# Rescaling our separate TUBE and PETAL datasets and combining them into virtual3D flowers

# Script written by Sara van de Keker (saravandekerke@outlook.com)
# Written December 2018, updated October 2019
# Accompanying van de Kerke et al. 2019 Capturing variation in floral shape; a virtual3D based morphospace for Pelargonium, PeerJ preprint

library(geomorph) # using version 3.1.2

### set working directory 
setwd("C:/Users/kerke001/OneDrive - WageningenUR/morphometrics_chapter_2-3/Rscripten/")



###########################################
############### Data import ###############
# Import the datalist created in the earlier steps (S2 Supplementary R script AnalyseFlowerShapaSeparate) 

load("pelargonium_TUBE_rotated.Rda")
load("pelargonium_PETAL_rotated.Rda")


# Import variablesfile containing only individuals for those species that show overlap between the datasets  
classifier.variables.PETAL.matching <- read.csv("Variables_PETAL_overlap.csv")
classifier.variables.TUBE.matching <- read.csv("variables_TUBE_overlap.csv", header=TRUE)



###########################################
############### Resize data ###############
# To be able to properly combine our datasets, we need as biologically accurate data as possible. We therefor want our individuals to
# be of accurate size, but not rotated in any way (because we then might accidentaly combine corolla with a tilted nectar tube). 
# The GPA coordinates for the separate datasets we have created earlier are already rotated in the GPA, but the size aspect has been
# removed as well in that step. We thus bring back the size component into the GPA coordinate data. For this, we use the Csize 
# extracted in the GPA analysis. We then end up with rotated individuals that are not corrected for size. 

### TUBE
# Bring back size by multiplying GPA coordinates with Csize
pelargonium.TUBE.rotated.initial <- 
  simplify2array(lapply(1:dim(pelargonium.TUBE.rotated$Coords.raw)[3], 
                        function(j) pelargonium.TUBE.rotated$Coords.GPA$coords[,,j] * 
                          pelargonium.TUBE.rotated$Coords.GPA$Csize[[j]]) )                                        

# Replace dimnames with correct ones
dimnames(pelargonium.TUBE.rotated.initial)[[3]] <- 
  dimnames(pelargonium.TUBE.rotated$Coords.raw)[[3]]

# Plot all specimens to check
plotAllSpecimens(pelargonium.TUBE.rotated$Coords.GPA$coords)
plotAllSpecimens(pelargonium.TUBE.rotated.initial)



### PETAL
# Bring back size by multiplying GPA coordinates with Csize
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]]) )                                        

# Replace dimnames with correct ones
dimnames(pelargonium.PETAL.rotated.initial)[[3]] <- 
  dimnames(pelargonium.PETAL.rotated$Coords.raw)[[3]]

# Plot all specimens to check
plotAllSpecimens(pelargonium.PETAL.rotated$Coords.GPA$coords)
plotAllSpecimens(pelargonium.PETAL.rotated.initial)



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

### Preparation

# First, we place the individuals from the variables files that are present in both datasets (complete overlap) in a new object.
petal<-classifier.variables.PETAL.matching
TUBE<-classifier.variables.TUBE.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]
TUBE$ID[TUBE$ID %in% dimnames(pelargonium.TUBE.rotated.initial)[[3]]==F]

# We make empty vectors to store the variable data in. 
individuals<-numeric(0)
nr.TUBE<-numeric(0)
nr.petal<-numeric(0)
sumnr.TUBE<-0
sumnr.petal<-0

# We determine for each species the nr. and sumnr.
# nr. holds the counted number of individuals in the original variable dataset
# sumnr. gives the rownumber each next species starts at in the variable dataset
for (i in 1:68) {
  nr.s<- sum(TUBE[,3]==i)
  nr.p<- sum(petal[,3]==i)
  
  nr.TUBE<- c(nr.TUBE,nr.s)
  nr.petal<-c(nr.petal, nr.p)
  sumnr.TUBE<-c(sumnr.TUBE, sum(nr.TUBE))
  sumnr.petal<-c(sumnr.petal, sum(nr.petal))
}

# Check
nr.TUBE
nr.petal

sumnr.TUBE
sumnr.petal



##### Start of the bootstrapping loop #####
# We randomly combine 6 specimens from the PETAL and TUBE datasets (the bootstrapping) per species and repeat that loop 20 times.

for (j in 1:20) {
  
  # We make empty vectors to store the bootstrap sampling data in.  
  bootTUBE.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) {
    bootTUBE.number <-c(bootTUBE.number,sumnr.TUBE[i]+sample(nr.TUBE[i],6, replace = TRUE))
    bootpetal.number <-c(bootpetal.number,sumnr.petal[i]+sample(nr.petal[i],6, replace = TRUE))
  }
  
  # Check content of bootTUBE and bootpetal
  bootTUBE.number
  bootpetal.number
  
  # Make sure the output of bootTUBE is in the same order as the variable data
  bootTUBE<-as.vector(TUBE$ID)[bootTUBE.number]
  bootpetal<-as.vector(petal$ID)[bootpetal.number]
  
  # Check
  head(bootTUBE)
  head(bootpetal)
  
  # Combine bootTUBE and bootpetal
  
  bbbb<-cbind(bootTUBE,bootpetal)
  
  # Connect numbers of bootTUBE/bootpetal to the individual numbers of the variable dataset and write to a file
  boots <- pelargonium.TUBE.rotated.initial[,,bootTUBE]
  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(bootTUBE)){
   
    # First, we extract the individual names
    
    a <- dimnames(boots)[[3]][i]
    b <- dimnames(bootp)[[3]][i]
    
    ### For the TUBE 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 TUBE 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 TUBE 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 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) usning 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
    
    
    # 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)
    
    
    ### Now the anchor point of the petal dataset (defined in refx2) is used to define the new x coordinate of the TUBE dataset. 
    # This x was previously defined to be zero, but has been 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 .txt file using 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 bootTUBE to species names of the variable dataset and write to a .csv file
  individuals.raw <- as.vector(TUBE$Species[bootTUBE.number])
  individuals <- cbind(j, i, bootTUBE, bootpetal, individuals.raw)
  write.table(individuals, file="individuals68species3.csv", append=TRUE, sep = ",", row.names = TRUE, col.names = FALSE)
  
  
  
}

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



# Import the .csv file created above containing individual species info of the VIRTUAL3D dataset and rename columns
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 ############### 
# We combine the 'raw' virtual3D coordinate data, the species names from the variables file, and information regarding the individuals
# of the original TUBE and PETAL datasets used in a list for further analyses
pelargonium.VIRTUAL3D.initial <- list(pelargonium.VIRTUAL3D.original, 
                                     classifier.individuals$species,
                                     classifier.individuals$combination)
names(pelargonium.VIRTUAL3D.initial) <- c("Coords",   "Species", "Combination")


save(pelargonium.VIRTUAL3D.initial, file = "Pelargonium_VIRTUAL3D_initial.Rda")