# This code cannot be reproduced without access to the underlying datasets
# It also relies on packages, functions, & variables that are called in SM.Rmd

#contains ugly fix due to simple ctc & lena ctc not controlling for length

remove_single_rec_kids = TRUE

for (data_set in data_sets){ #data_set="aclew";data_set="lena"
  if(data_set=="aclew"){
    mydat <- read.csv(paste0('../input/el1000-metrics/output/', data_set,'_metrics.csv'))
    #ugly fix for simple CTC & lena CTC not controlling for length
    mydat$simple_CTC_ph = mydat$simple_CTC/(mydat$duration_vtc/(60*60*1000))
  } else{
    mydat <- read.csv(paste0('../input/el1000-metrics/output/', data_set,'_metrics_avass.csv'))
    mydat=mydat[,colnames(mydat)!="X"] #for some reason, this version has row names - remove that col to avoid having issues downastream
    #ugly fix for simple CTC & lena CTC not controlling for length
    mydat$lena_CTC_ph=mydat$lena_CTC/(mydat$duration_its/(60*60*1000))
    mydat$lena_CVC_ph=mydat$lena_CVC/(mydat$duration_its/(60*60*1000))
  }
  
  # Remove Cougar non-normatives
  # if needed, in the terminal, do `datalad get input/el1000-metrics/EL1000/cougar/metadata/children.csv`
  cougar <- read.csv('../input/el1000-metrics/EL1000/cougar/metadata/children.csv')
  cougar_normative <- cougar[cougar$normative == "Y", ]
  mydat <- mydat[mydat$experiment != 'cougar' | mydat$child_id %in% cougar_normative$child_id, ]
  
  if(data_set=="aclew"){
    mydat2 <- read.csv(paste0('../input/laac-metrics/output/', data_set,'_metrics.csv'))
    #ugly fix for simple CTC & lena CTC not controlling for length
    mydat2$simple_CTC_ph = mydat2$simple_CTC/(mydat2$duration_vtc/(60*60*1000))
  } else{
    mydat2 <- read.csv(paste0('../input/laac-metrics/output/', data_set,'_metrics_avass.csv'))
    mydat2=mydat2[,colnames(mydat2)!="X"] #for some reason, this version has row names - remove that col to avoid having issues downastream
    #ugly fix for simple CTC & lena CTC not controlling for length
    mydat2$lena_CTC_ph=mydat2$lena_CTC/(mydat2$duration_its/(60*60*1000))
    mydat2$lena_CVC_ph=mydat2$lena_CVC/(mydat2$duration_its/(60*60*1000))
  }
  # Remove FauseyElse
  # if needed, in the terminal, do `datalad get input/laac-metrics/datasets/fausey-trio/metadata/recordings.csv`
  fausey_trio <- read.csv('../input/laac-metrics/datasets/fausey-trio/metadata/recordings.csv')
  fausey_trio_full <- fausey_trio[fausey_trio$Trio_Subset == "Trio_Full", ]
  fausey_trio_full$session_id <- paste0(fausey_trio_full$HomeBank_ID, "/", fausey_trio_full$fileName)
  mydat2 <- mydat2[mydat2$experiment=="quechua" | mydat2$session_id %in% fausey_trio_full$session_id,]
  #note that since we are only taking quechua & fausey-trio, then the fact that png2019 & tsimane2017 are repeated across laac & el1000 is not a problem
  
  #columns are not in the same order across el1000-metrics & laac-metrics
  mydat2=mydat2[,colnames(mydat)]
  
  
  mydat=rbind(mydat,mydat2)
  
  #ugly fix removing columns that we do not want to include
  mydat=mydat[,!(colnames(mydat)%in%c("lena_CTC","simple_CTC","lena_CVC","voc_chi"))]
  
  #remove no overlap metrics
  mydat = mydat[,grep("noov",colnames(mydat),invert=T)]
  mydat = mydat[,grep("no_overlap",colnames(mydat),invert=T)]
  
  #remove not included corpora
  mydat <- mydat[is.element(mydat$experiment, corpora),]
  
  #make sure child_id & session_id is unique across corpora
  mydat$child_id=paste(mydat$experiment,mydat$child_id)
  mydat$session_id=paste(mydat$child_id,mydat$session_id)
  
  #table(mydat$session_id)[table(mydat$session_id)>1] #no repeated sessions
  
  if(remove_single_rec_kids){
    #remove short recs
    if(data_set=='aclew') dur_col="duration_vtc" else dur_col="duration_its"
    mydat=mydat[mydat[,dur_col]/3.6e+6>=4,]
    
    #remove kids with only one rec
    rec_per_child = setNames(aggregate(data = mydat, session_id ~ child_id, function(session_id) length(unique(session_id))), c('child_id', 'n'))
    mydat=mydat[mydat$child_id %in% rec_per_child$child_id[rec_per_child$n>1],]
  }
  
  # Save data
  write.csv(mydat,paste0('../data_output/', data_set,'_base_data_set.csv'),row.names = F)
  print(paste0('Save to ', paste0('../data_output/', data_set,'_base_data_set.csv')))
}


for (data_set in data_sets){
  mydat <- read.csv(paste0('../data_output/', data_set,'_base_data_set.csv'))
  mydat$age_s=scale(mydat$age)
  mydat$age_s=(mydat$age - mean(mydat$age , na.rm=T))/sd(mydat$age , na.rm=T)

  metrics <- colnames(mydat)[!is.element(colnames(mydat), no.scale.columns)]
  
  #remove outliers
  for(metric in metrics) mydat[abs((mydat[,metric]-mean(mydat[,metric], na.rm=T))/sd(mydat[,metric], na.rm=T)) > 2.5 & !is.na(abs((mydat[,metric]-mean(mydat[,metric], na.rm=T))/sd(mydat[,metric], na.rm=T))), metric]<-NA 
  #NA values that are beyond 2.5 SD from mean
  
    
  write.csv(mydat,paste0('../data_output/', data_set,'_metrics.csv'),row.names = F) #all variables are unscaled, except for age
  print(paste0('Save to ', paste0('../data_output/', data_set,'_metrics.csv')))
  
  for(metric in metrics){ #metric="pc_mal_ph"
   # Scale
    mydat[, metric] <- (mydat[, metric] - mean(mydat[, metric], na.rm=T)) /sd(mydat[, metric], na.rm=T)
  }
  
  # Save data
  write.csv(mydat,paste0('../data_output/', data_set,'_metrics_scaled.csv'),row.names = F)
  print(paste0('Save to ', paste0('../data_output/', data_set,'_metrics_scaled.csv')))
}