#########################################################################################################################################################
#
# Import, process and export in a standardized Newick format various language classifications (currently WALS, Ethnologue, Glottolog and AUTOTYP)
# and add branch length using a variety or methods (none, constant, proportional, grafen, nnls, ga and nj; see script "TreeHelperFunctions.r" for info)
#
# Copyright (C) 2013-2017  Dan Dediu
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
#
##########################################################################################################################################################

######################################################################################
#
# Global parametres controlling the various processing, import and export options:
#    tweak as desired (some require lots of computater resources and time)
# More parametres controlling the branch length methods are defined below.
#
######################################################################################
MATCH_CODES       = TRUE; # compute the equivalences between the ISO, WALS, AUTOTYP and GLOTTOLOG codes and save them to file?
PREOPTIMIZE_DISTS = TRUE; # pre-optimize the distance matrices for fast loading when required
COMPUTE_GEO_DISTS = TRUE; # compute the geographic distances between languages?

TRANSFORM_TREES   = TRUE; # transform the trees to the Newick notation (no branch length) and save them?

EXPORT_NEXUS      = TRUE;  # export the trees to a NEXUS file?
EXPORT_NEXUS_TRANSLATE_BLOCK = TRUE; # when exporting NEXUS files, generate a TRANSLATE block (useful when using programs such as BayesTraits that have issues parsing complicated taxa names)?
EXPORT_CSV        = TRUE;  # export the trees to a CSV file?

COMPUTE_BRLEN     = TRUE; # apply the various branch length methods to the Newick no branch length trees?

CPU_CORES = 4; # the number of parallel CPU cores to use overall

quotes="'"; # the quotes to use for the node (language and family) names

#########################
#
#    Tree functions
#
#########################
source( "./FamilyTrees.R" );


######################################################################################
#
# Data needed to map the ISO - WALS - Glottolog codes
#
######################################################################################

if( MATCH_CODES )
{
  # The Ethnologue data contains the ISO and language names:
  ethn.data <- read.table( "../input/ethnologue/LanguageCodes.tab", header=TRUE, sep="\t", quote="", stringsAsFactors=FALSE );
  ethn.data <- ethn.data[ , c("LangID", "Name") ];
  ethn.data$Name <- sapply( as.character(ethn.data$Name), function( s ){ normalize.language.name(s); } );
  
  # The wals data contains the WALS, ISO and Glottolog codes for a subset of languages:
  wals.data <- read.table( "../input/wals/languages_values.csv", header=TRUE, sep=",", quote="\"", stringsAsFactors=FALSE );
  wals.data <- wals.data[ , c("wals_code", "iso_code", "glottocode", "Name", "latitude", "longitude") ];
  wals.data$Name <- sapply( as.character(wals.data$Name), function( s ){ normalize.language.name(s); } );
  
  # The Glottolog data contains the ISO and Glottolog codes for a much larger set of languages:
  glott.data <- read.csv("../input/glottolog/languages_and_dialects_geo.csv");
  glott.data <- glott.data[, c("name", "glottocode", "isocodes", "latitude", "longitude")]; names(glott.data) <- c("Name", "Glottocode", "ISO", "latitude", "longitude");
  glott.data$Name <- sapply( as.character(glott.data$Name), function( s ){ normalize.language.name(s); } );
  glott.data$ISO[ glott.data$ISO == "" ] <- NA; 
  
  # The autotyp data contains the ISO, Glottolog and AUTOTYP (LID) mappings:
  autotyp.data <- read.table( "../input/autotyp/autotyp-trees.csv", header=TRUE, sep="\t", quote="", stringsAsFactors=FALSE );
  autotyp.data$language <- sapply( as.character(autotyp.data$language), function( s ){ normalize.language.name(s); } );
  
  # The mappings:
  code.mappings <- merge( glott.data, wals.data, by.x=c("Glottocode","ISO"), by.y=c("glottocode","iso_code"), all=TRUE, suffixes=c(".glot",".wals") );
  code.mappings <- merge( code.mappings, autotyp.data, by.x=c("Glottocode","ISO"), by.y=c("glottolog_LID.2014","ISO639.3.2013"), all=TRUE, suffixes=c(".glot",".auto") );
  code.mappings <- merge( code.mappings, ethn.data, by.x=c("ISO"), by.y=c("LangID"), all=TRUE, suffixes=c("",".ethn") );
  # Combine and harmonize the geographic corrdinates:
  # Fill in the missing coordinates in Glottolog from WALS:
  s <- (is.na(code.mappings$latitude.glot) | is.na(code.mappings$longitude.glot));
  code.mappings[s,c("latitude.glot","longitude.glot")] <- code.mappings[s,c("latitude.wals","longitude.wals")];
  # If the Glottolog and WALS coordinates are very different, display the differences and keep the WALS:
  s <- (!is.na(code.mappings$latitude.glot) & !is.na(code.mappings$longitude.glot) & !is.na(code.mappings$latitude.wals) & !is.na(code.mappings$longitude.wals) & 
        (abs(code.mappings$latitude.glot - code.mappings$latitude.wals) >= 1.0 | abs(code.mappings$latitude.glot - code.mappings$latitude.wals) >= 1.0));
  if( sum(s,na.rm=TRUE) > 0 )
  {
    cat( "There are disagreements larger than 1 degree between Glottolog and WALS for ", nrow(unique(code.mappings[s,c("ISO","wals_code","Glottocode")])), " languages (WALS codes): ", paste( unique(code.mappings$wals_code[s]), collapse=",", sep="" ), "; keeping the Glottolog coordinates!\n" );
    #code.mappings[s,c("latitude.glot","longitude.glot")] <- code.mappings[s,c("latitude.wals","longitude.wals")];
  }
  # See how many languages with ISO codes but no geographic coordinates there are:
  code.mappings$ISO[ is.na(code.mappings$ISO) ] <- "";
  s <- ((is.na(code.mappings$latitude.glot) | is.na(code.mappings$longitude.glot)) & code.mappings$ISO != "");
  if( sum(s,na.rm=TRUE) )
  {
    missing.geo <- unique(code.mappings$ISO[s]);
    cat( "There are ", length(missing.geo), " languages with ISO codes but no geographic information: I don't know where to get the coordinates from :(\n" );
  }
  # Data cleaning:
  code.mappings <- code.mappings[ , c( "ISO", "wals_code", "LID", "Glottocode", "Name", "Name.wals", "language", "Name.glot", "latitude.glot", "longitude.glot" ) ]; 
  names(code.mappings) <- c( "ISO", "WALS", "AUTOTYP", "Glottolog", "Name.ethn", "Name.wals", "Name.autotyp", "Name.glottlog", "Latitude", "Longitude" ); 
  code.mappings[is.na(code.mappings)] <- "";
  code.mappings$ISO           <- as.character(code.mappings$ISO);
  code.mappings$WALS          <- as.character(code.mappings$WALS);
  code.mappings$AUTOTYP       <- as.character(code.mappings$AUTOTYP);
  code.mappings$Glottolog     <- as.character(code.mappings$Glottolog);
  code.mappings$Name.ethn     <- as.character(code.mappings$Name.ethn);
  code.mappings$Name.wals     <- as.character(code.mappings$Name.wals);
  code.mappings$Name.autotyp  <- as.character(code.mappings$Name.autotyp);
  code.mappings$Name.glottlog <- as.character(code.mappings$Name.glottlog);
  code.mappings$Latitude      <- as.numeric(code.mappings$Latitude);
  code.mappings$Longitude     <- as.numeric(code.mappings$Longitude);
  code.mappings <- code.mappings[ order( code.mappings$ISO, code.mappings$WALS, code.mappings$AUTOTYP, code.mappings$Glottolog ), ];
  
  # Add a Universal Unique Language Identifier using the conventions introduced here:
  code.mappings <- cbind( code.mappings, "UULID"=paste("[i-",code.mappings$ISO,"][w-",code.mappings$WALS,"][a-",code.mappings$AUTOTYP,"][g-",code.mappings$Glottolog,"]",sep="") );
    
  # Write them to file:
  write.table( code.mappings, "../output/code_mappings_iso_wals_autotyp_glottolog.csv", sep="\t", quote=FALSE, row.names=FALSE );
  code.mappings$AUTOTYP <- as.character(code.mappings$AUTOTYP); code.mappings$AUTOTYP[ is.na(code.mappings$AUTOTYP) ] <- ""; # make sure the AUTOTYP code is treated as character with missing data as ""
} else
{
  code.mappings <- read.table( "../output/code_mappings_iso_wals_autotyp_glottolog.csv", sep="\t", quote="", header=TRUE, stringsAsFactors=FALSE );
  code.mappings$AUTOTYP <- as.character(code.mappings$AUTOTYP); code.mappings$AUTOTYP[ is.na(code.mappings$AUTOTYP) ] <- ""; # make sure the AUTOTYP code is treated as character with missing data as ""
}


######################################################################################
#
# Load (and possibly compute) various distance matrices to be used for computing
# the branch lengths.
#
######################################################################################

# Pre-optimize and cache the distance matrices for on-demand loading only when really needed:
if( PREOPTIMIZE_DISTS )
{
  # Make sure the ../output/preoptimized-distances directory exists:
  dir.create("../output/preoptimized-distances", showWarnings=FALSE);
  
  # ASJP16 (the row and column names are the iso codes):
  load( "../input/distances/ASJP/asjp16-dists.RData" ); # asjp16.dm
  d <- asjp16.dm; save(d, file="../output/preoptimized-distances/asjp16-dm.RData", compress="xz"); rm(asjp16.dm,d); # cache it pre-optimized for on-demand loading
  
  # WALS (the row and column names are the wals codes); gowler and manhattan are quite similar; using both NAs and imputation of NAs with the mode (per variable):
  load( "../input/distances/WALS/wals-gower-dm.RData" ); wals.gower.dm <- as.matrix(wals.gower.dm); # wals.gower.dm (with missing data)
  d <- wals.gower.dm; save(d, file="../output/preoptimized-distances/wals-gower-dm.RData", compress="xz"); rm(wals.gower.dm,d); # cache it pre-optimized for on-demand loading
  
  load( "../input/distances/WALS/wals-gower-mode-dm.RData" ); wals.gower.mode.dm <- as.matrix(wals.gower.mode.dm); # wals.gower.mode.dm (with mode imputation for missing data)
  d <- wals.gower.mode.dm; save(d, file="../output/preoptimized-distances/wals-gower-mode-dm.RData", compress="xz"); rm(wals.gower.mode.dm,d); # cache it pre-optimized for on-demand loading
  
  load( "../input/distances/WALS/wals-euclidean-dm.RData" ); wals.euclidean.dm <- as.matrix(wals.euclidean.dm); # wals.euclidean.dm (with missing data)
  d <- wals.euclidean.dm; save(d, file="../output/preoptimized-distances/wals-euclidean-dm.RData", compress="xz"); rm(wals.euclidean.dm,d); # cache it pre-optimized for on-demand loading
  
  load( "../input/distances/WALS/wals-euclidean-mode-dm.RData" ); wals.euclidean.mode.dm <- as.matrix(wals.euclidean.mode.dm); # wals.euclidean.mode.dm (with mode imputation for missing data)
  d <- wals.euclidean.mode.dm; save(d, file="../output/preoptimized-distances/wals-euclidean-mode-dm.RData", compress="xz"); rm(wals.euclidean.mode.dm,d); # cache it pre-optimized for on-demand loading
  
  # AUTOTYP (the row and column names are the AUTOTYP codes (LID)):
  load("../input/distances/AUTOTYP/autotyp-dist.RData"); # autotyp.dm (with missing data, using ony features with single values per language, courtesy of Balthasar Bickel)
  d <- autotyp.dm; save(d, file="../output/preoptimized-distances/autotyp-dm.RData", compress="xz"); rm(autotyp.dm,d); # cache it pre-optimized for on-demand loading
  
  # Geographic distances (the row and column names are the Glottolog codes):
  if( COMPUTE_GEO_DISTS )
  {
    # Compute the great-circle geographic distances between all languages with geographic coordinates
    library(geosphere);
    # The languages with geographic coordinates:
    s <- (!is.na(code.mappings$Latitude) & !is.na(code.mappings$Longitude)); 
    # of which those that don't have a specific type of code:
    s.ISO       <- s & (code.mappings$ISO == "");
    s.WALS      <- s & (code.mappings$WALS == "");
    s.AUTOTYP   <- s & (code.mappings$AUTOTYP == "");
    s.Glottolog <- s & (code.mappings$Glottolog == "");
    # Decide which code to use so that as many as possible of the languages with geographic info are covered by the coding scheme:
    code.to.use <- c("ISO","WALS","AUTOTYP","Glottolog")[ which.min( c( sum(s.ISO,na.rm=TRUE), sum(s.WALS,na.rm=TRUE), sum(s.AUTOTYP,na.rm=TRUE), sum(s.Glottolog,na.rm=TRUE) ) ) ];
    s.use <- s & !(get( paste( "s.", code.to.use, sep="" ) )); # obtain (s & (code != ""))
    tmp <- code.mappings[s.use,c(code.to.use,"Longitude","Latitude")];
    cat( "For the geographic distances (ellipsoid) the row and column names use the ", code.to.use, " codes.\n" );
    # Some codes (with geographic coordinates) are duplicated: keep the frst entries only:
    tmp <- tmp[!duplicated(tmp[,1]),];
    
    ## Plot these languages:
    #library(maps);
    #X11(10,10);
    #map("world");
    #points( tmp$Longitude, tmp$Latitude, cex=0.5, col=grey(0.5) );
    
    # Compute the distances:
    m <- as.matrix( tmp[,c("Longitude","Latitude")] ); rownames(m) <- tmp[,1];
    system.time( geo.dm <- geosphere::distm( m, fun=distMeeus ) ); # much faster than distVincentyEllipsoid but still a very decent estimate of geographic distance
    rownames(geo.dm) <- colnames(geo.dm) <- rownames(m);
    
    rm(m); rm(tmp);
    
    save( geo.dm, file="../input/distances/geo-great-circle-ellipsoid-dists.RData", compress="xz" );
  }
  load( file="../input/distances/geo-great-circle-ellipsoid-dists.RData" ); # geo.dm 
  # Normalize it between 0 and 1 by dividing it with the maximum distance (which is very close to the maximum possible distance on Earth of about 21,000km):
  geo.dm <- geo.dm / max(as.numeric(geo.dm),na.rm=TRUE);
  d <- geo.dm; save(d, file="../output/preoptimized-distances/geo-dm.RData", compress="xz"); rm(geo.dm,d); # cache it pre-optimized for on-demand loading
  
  
  # The MG2015 [Maurits, L. & Griffiths, T.L. (2015) PNAS 111 (37):13576--13581] distances
  # there's one distance matrix for each of the four classifications (WALS, Ethnologue, Glottolog and AUTOTYP)
  # and each has a rwo & column names the full language names and code -> we need to replace them by the appropriate codes
  .change.row.col.names.to.codes <- function(m, code=c("iso","wals","autotyp","glottolog"))
  {
    if( is.null(m) || !inherits(m,"matrix") || nrow(m) != ncol(m) || any(rownames(m) != colnames(m)) )
    { 
      stop("Problems with the distance matrix!\n"); 
      return (NULL); 
    } else
    {
      # Change the row and column names to the WALS codes:
      new.names <- vapply(rownames(m), function(s)
        { 
          tmp <- extract.name.and.codes(s); 
          if( is.null(tmp[[code]]) || length(tmp[[code]]) != 1 || tmp[[code]]=="" )
          { 
            warning(paste0("Language ",s," must have a single ", code, " code\n")); 
            return (s); 
          } else 
          {
            return(tmp[[code]]); 
          }
        }, character(1));
      return (new.names);
    }
  }
  
  load("../input/distances/MG2015/MG2015-wals-alpha=0.69.RData"); rownames(MG2015.WALS) <- colnames(MG2015.WALS) <-.change.row.col.names.to.codes(MG2015.WALS, "wals");  # MG2015.WALS
  d <- MG2015.WALS; save(d, file="../output/preoptimized-distances/mg2015-wals-dm.RData", compress="xz"); rm(MG2015.WALS,d); # cache it pre-optimized for on-demand loading
  
  load("../input/distances/MG2015/MG2015-ethnologue-alpha=0.69.RData"); rownames(MG2015.Ethnologue) <- colnames(MG2015.Ethnologue) <-.change.row.col.names.to.codes(MG2015.Ethnologue, "iso");  # MG2015.Ethnologue
  d <- MG2015.Ethnologue; save(d, file="../output/preoptimized-distances/mg2015-ethnologue-dm.RData", compress="xz"); rm(MG2015.Ethnologue,d); # cache it pre-optimized for on-demand loading
  
  load("../input/distances/MG2015/MG2015-glottolog-alpha=0.69.RData"); rownames(MG2015.Glottolog) <- colnames(MG2015.Glottolog) <-.change.row.col.names.to.codes(MG2015.Glottolog, "glottolog");  # MG2015.Glottolog
  d <- MG2015.Glottolog; save(d, file="../output/preoptimized-distances/mg2015-glottolog-dm.RData", compress="xz"); rm(MG2015.Glottolog,d); # cache it pre-optimized for on-demand loading
  
  load("../input/distances/MG2015/MG2015-autotyp-alpha=0.69.RData"); rownames(MG2015.AUTOTYP) <- colnames(MG2015.AUTOTYP) <-.change.row.col.names.to.codes(MG2015.AUTOTYP, "autotyp");  # MG2015.AUTOTYP
  d <- MG2015.AUTOTYP; save(d, file="../output/preoptimized-distances/mg2015-autotyp-dm.RData", compress="xz"); rm(MG2015.AUTOTYP,d); # cache it pre-optimized for on-demand loading
  
  gc(); # call the garbage collector to make sure the space is freed after this memory-hungry step
}


######################################################################################
#
# Global parametres controlling the branch length methods to apply
#
######################################################################################

# the branch length methods to be used; can be "all", "none", or any subset of {"constant", "proportional", "grafen", "nnls", "ga", "nj"}
CLASSIFICATIONS = c("wals", "glottolog"); #c("wals", "ethnologue", "glottolog", "autotyp");  --> update just WALS and Glottolog trees, but don't touch the others!
METHODS   = c("constant", "proportional", "grafen", "nnls", "nj"); # , "ga"); --> GA has issues: don't update it for now!  
CONSTANT  = 1.0; # the positive constant required by some methods
DISTS.CODES = read.table(text="Distance              ShortName  File                                                           Code      \n
                               asjp16                asjp       ../output/preoptimized-distances/asjp16-dm.RData               iso       \n
                               wals(gower)           w:g        ../output/preoptimized-distances/wals-gower-dm.RData           wals      \n
                               wals(gower,mode)      w:gm       ../output/preoptimized-distances/wals-gower-mode-dm.RData      wals      \n
                               wals(euclidean)       w:e        ../output/preoptimized-distances/wals-euclidean-dm.RData       wals      \n
                               wals(euclidean,mode)  w:em       ../output/preoptimized-distances/wals-euclidean-mode-dm.RData  wals      \n
                               autotyp               atyp       ../output/preoptimized-distances/autotyp-dm.RData              autotyp   \n
                               mg2015(wals)          mg:w       ../output/preoptimized-distances/mg2015-wals-dm.RData          wals      \n
                               mg2015(ethnologue)    mg:e       ../output/preoptimized-distances/mg2015-ethnologue-dm.RData    iso       \n
                               mg2015(glottolog)     mg:g       ../output/preoptimized-distances/mg2015-glottolog-dm.RData     glottolog \n
                               mg2015(autotyp)       mg:a       ../output/preoptimized-distances/mg2015-autotyp-dm.RData       autotyp   \n
                               geo                   geo        ../output/preoptimized-distances/geo-dm.RData                  glottolog \n",
                         sep="", quote="", header=TRUE, stringsAsFactors=FALSE); # the distances required by some methods

######################
#
# WALS trees
#
######################

if( TRANSFORM_TREES )
{
  # WALS languages info:
  wals.data <- read.table( "../input/wals/languages_values.csv", header=TRUE, sep=",", quote="\"", stringsAsFactors=FALSE );
  wals.data <- wals.data[, c("wals_code", "iso_code", "glottocode", "Name", "latitude", "longitude", "genus", "family" ) ];
  wals.data <- wals.data[ !is.na(wals.data$iso_code) & wals.data$iso_code != "", ]; # keep only the actual languages
  
  # Replace troublesome characters in language and family names:
  wals.data$Name <- sapply( as.character(wals.data$Name), function( s ){ normalize.language.name(s); } );
  wals.data$family <- sapply( as.character(wals.data$family), function( s ){ normalize.language.name(s); } );
  wals.data$genus <- sapply( as.character(wals.data$genus), function( s ){ normalize.language.name(s); } );
  
  # Create the WALS trees:
  roots <- NULL;
  for( i in 1:nrow(wals.data) )
  {
    # Process this entry:
    if( !is.na(wals.data$family[i]) && !(wals.data$family[i] %in% c("other","Unclassified","Creole","Mixed language")) )
    {
      autotypcode <- unique( as.character( code.mappings$AUTOTYP[ code.mappings$WALS == as.character(wals.data$wals_code[i]) ] ) );
      if( length(autotypcode) > 1 )
      {
        # Get rid of empty codes:
        autotypcode <- autotypcode[ autotypcode != "" ];
        if( length(autotypcode) == 0 )
        {
          autotypcode = ""; # all were ""
        } else
        {
          autotypcode <- paste( autotypcode, collapse="-", sep="" ); # in principle it can happen to have more than one AUTOTYP code
        }
      }
      roots <- add.language.to.languageclassification( roots, 
                                                       path=c(create.name.and.codes(wals.data$family[i], iso="", wals="", autotyp="", glottocode="", quotes=quotes),  
                                                              create.name.and.codes(wals.data$genus[i],  iso="", wals="", autotyp="", glottocode="", quotes=quotes),
                                                              create.name.and.codes(wals.data$Name[i], 
                                                                                    iso=       as.character(wals.data$iso_code[i]), 
                                                                                    wals=      as.character(wals.data$wals_code[i]), 
                                                                                    autotyp=   autotypcode, 
                                                                                    glottocode=as.character(wals.data$glottocode[i]), 
                                                                                    quotes) ),
                                                       brlens=NULL );
    }
  }
  attr(roots, "classif.name") <- "wals";
  # Fix the family and language names (get rid of the codes):
  for( i in 1:length(roots) )
  {
    roots[[i]] <- .fix.names(roots[[i]], quotes=quotes);
    attr(roots[[i]], "tree.name") <- str_trim(strsplit(strsplit(get.name(roots[[i]]),quotes,fixed=TRUE)[[1]][2],"[",fixed=TRUE)[[1]][1]);
  }
  # Save them to file:
  export.languageclassification( roots, dir.name="../output", classification="wals", export.nexus=EXPORT_NEXUS, nexus.translate.block=EXPORT_NEXUS_TRANSLATE_BLOCK, export.csv=EXPORT_CSV );
  
  rm(wals.data);
}


######################
#
# Ethnologue trees
#
######################

if( TRANSFORM_TREES && FALSE ) # Ethnologue has changed too much and it is not such an important language classification system anymore, so don't touch it.
{
  # The Enthologue trees: we need to download and parse them from the Ethnologue website:
  fams.file <- paste("../input/ethnologue/fams.html",sep="");
  fams.URL <- "http://www.ethnologue.com/browse/families";
  if( !file.exists(fams.file) || file.info(fams.file)$size < 1 )
  {
    options(timeout=180); # Ethnologue is a slow site for some reason...
    options(HTTPUserAgent="Mozilla/5.0 (X11; Linux x86_64; rv:17.0) Gecko/20130917 Firefox/17.0"); # emulate Firefox
    if( download.file( fams.URL, fams.file, quiet=T ) != 0 )
    {
      stop( "Error retrieving the Ethnologue main families page !\n" );
    }
  }
  # Load the file:
  HTML.content <- readLines( fams.file, warn=FALSE ); # read the file:
  # Find the links to the families themselves:
  fams.links <- HTML.content[ grep( "/subgroups/", HTML.content, fixed=TRUE ) ];
  # Extract the links:
  fams.links <- sapply( fams.links, function(s){ 
                                                 tmp <- strsplit( s, "<a href=\"/subgroups/", fixed=TRUE )[[1]];
                                                 if( length(tmp) < 2 )
                                                 {
                                                   cat( "Error retrieving the URL for ", s, "\n" );
                                                   return (NULL);
                                                 }
                                                 tmp <- strsplit( tmp[2], "\">", fixed=TRUE )[[1]];
                                                 if( length(tmp) < 1 )
                                                 {
                                                   cat( "Error retrieving the URL for ", s, "\n" );
                                                   return (NULL);
                                                 }
                                                 return (tmp[1]);
                                               } );
  
  # Load and process one language family
  fetch.ethnologue.family <- function( family )
  {
    if( is.na( family ) )
    {
      return( NULL );
    }
    else
    {
      # The base URL for families:
      fam.base.URL <- "http://www.ethnologue.com/subgroups/";
      fam.URL <- paste( fam.base.URL, family, sep="" );
      fam.file <- paste("../input/ethnologue/fam-",family,".html",sep="");
      if( !file.exists(fam.file) || file.info(fam.file)$size < 1 )
      {
        options(timeout=180); # Ethnologue is a slow site for some reason...
        options(HTTPUserAgent="Mozilla/5.0 (X11; Linux x86_64; rv:17.0) Gecko/20130917 Firefox/17.0"); # emulate Firefox
        if( download.file( fam.URL, fam.file, quiet=T ) != 0 )
        {
          cat( "Error retrieving family ", family, " !\n" );
          return( NULL );
        }
      }
      
      # Use the downloaded file:
      HTML.content <- readLines( fam.file, warn=F ); # read the file:
      
      # Get the family official name:
      tag.string <- "<h1 class=\"title\" id=\"page-title\">";
      matches <- grep( tag.string, HTML.content, fixed=TRUE ); 
      if( length(matches) != 1 )
      {
        cat( "Error locating family name for ", family, "\n" )
        return ( NULL );
      } else
      {
        # Get the family name:
        family.name <- strsplit( strsplit( HTML.content[ matches ], tag.string, fixed=TRUE )[[1]][2], "</h1>", fixed=TRUE )[[1]][1];
      }
      
      # The tree in preorder with special symbols ( and ) representing an increse/decrease in level
      # the names are "node" for an internal node, "leaf" for a leaf, and "special" for ( and ):
      tree <- c("node"=create.name.and.codes( family.name, iso="", wals="", autotyp="", glottocode="", quotes=quotes));
      
      # Preprocess the HTML text to make sure each of the important tags appear only once per line (i.e., insert newlines after each </a> and </ul>):
      HTML.content.new <- gsub("</a>", "</a>\n", HTML.content[(matches+1):length(HTML.content)], fixed=TRUE);
      HTML.content.new <- gsub("</ul>", "</ul>\n", HTML.content.new, fixed=TRUE);
      HTML.content.new <- unlist(strsplit(HTML.content.new, "\n", fixed=TRUE));
      
      # Build the tree by adding the subfamilies and languages as they appear in the file:
      for( i in 1:length(HTML.content.new) )
      {
        if( length( grep( "<a href=\"/subgroups/", HTML.content.new[i], fixed=TRUE ) ) == 1 )
        {
          # This seems to be a subgroup: extract its name:
          subgroup <- strsplit( strsplit( strsplit( HTML.content.new[i], "<a href=\"/subgroups/", fixed=TRUE )[[1]][2], "\">", fixed=TRUE )[[1]][2], "</a>", fixed=TRUE )[[1]][1];
          subgroup <- str_trim( strsplit( subgroup, "\\([[:digit:]]+\\)" )[[1]][1] );
          if( is.null(subgroup) || subgroup == "" )
          {
            cat( "Error extracting subgroup from ", HTML.content.new[i], "\n" );
            return ( NULL );
          }
          # Add this to the tree:
          tree <- c(tree, "node"=create.name.and.codes(subgroup, iso="", wals="", autotyp="", glottocode="", quotes=quotes)); 
        } else if( length( grep( "<span class=\"field-content\">", HTML.content.new[i], fixed=TRUE ) ) == 1 && length( grep( "<a href=\"/language/", HTML.content.new[i], fixed=TRUE ) ) == 1 )
        {
          # Seems to be language: extract it and its iso:
          lg.name <- strsplit( HTML.content.new[i], "<span class=\"field-content\">", fixed=TRUE )[[1]][2];
          lg.name <- str_trim( strsplit( lg.name, "<a href=\"/language/", fixed=TRUE )[[1]][1] );
          lg.iso <- strsplit( HTML.content.new[i], "<a href=\"/language/", fixed=TRUE )[[1]][2];
          lg.iso <- strsplit( lg.iso, "\">[", fixed=TRUE )[[1]][2];
          lg.iso <- strsplit( lg.iso, "]</a>", fixed=TRUE )[[1]][1];
          if( is.null(lg.name) || lg.name == "" ||  is.null(lg.iso) || lg.iso == "" )
          {
            cat( "Error extracting language from ", HTML.content.new[i], "\n" );
            return ( NULL );
          }
          # Get the other codes as well:
          wals <- autotyp <- glottocode <- "";
          s <- code.mappings[ code.mappings$ISO == lg.iso, ];
          if( !is.null(s) && nrow(s) > 0 )
          {
            wals       <- unique(s$WALS);
            autotyp    <- unique(s$AUTOTYP);
            glottocode <- unique(s$Glottolog);
          }
          # Add this language to the tree:
          tree <- c(tree, "leaf"=create.name.and.codes(lg.name, iso=lg.iso, wals=wals, autotyp=autotyp, glottocode=glottocode, quotes=quotes));
        } else if( length( grep( "<ul>", HTML.content.new[i], fixed=TRUE ) ) == 1 )
        {
          # Subgroup starts:
          tree <- c(tree, "special"="(");
        } else if( length( grep( "</ul>", HTML.content.new[i], fixed=TRUE ) ) == 1 )
        {
          # Subgroup ends:
          tree <- c(tree, "special"=")");
        } else if( length( grep( "</section>", HTML.content.new[i], fixed=TRUE ) ) == 1 )
        {
          # The end of the tree section
          break;
        }
      }
      
      # Build the tree reusing as much code as possible: simply add the full path for each language:
      root <- NULL;
      cur.path <- tree[1];
      for( i in 2:length(tree) )
      {
        #cat(paste0("Tree[",i,"]=",tree[i],": "));
        if( names(tree)[i] == "special" && tree[i] == "(" )
        {
          # Simply ignore, we're going up one level
        } else if( names(tree)[i] == "special" && tree[i] == ")" )
        {
          if( i < length(tree) && names(tree)[i+1] == "special" && tree[i+1] == "(" ) next; # artifact of the way languages and sub-families are mixed in the HTML file
          # Down one level, discard the last node in the path:
          cur.path <- cur.path[ -length(cur.path) ];
          #cat("cur.path=",paste0(cur.path,collapse=", "));
        } else if( names(tree)[i] == "leaf" )
        {
          # Add the full path to the tree:
          root <- add.tree.path(root, c(cur.path, tree[i]));
          #cat(paste0("add.tree.path(",paste0(c(cur.path, tree[i]),collapse=", "),")"));
        } else if( names(tree)[i] == "node" )
        {
          # Add the node to the path:
          cur.path <- c(cur.path, tree[i]);
          #cat("cur.path=",paste0(cur.path,collapse=", "));
        }
        #cat("\n");
      }
      attr(root, "tree.name") <- family.name;
      
      return( root );
    }
  }
  fetch.ethnologue.family <- cmpfun(fetch.ethnologue.family);
  
  # Fetch the families:
  roots <- list();
  for( i in 1:length(fams.links) )
  {
    cat( "Processing Ethnologue family ", fams.links[i], "...\n" );
    root <- fetch.ethnologue.family( fams.links[i] );
    if( is.null(root) )
    {
      cat( "  >>>  Error fetching Ethnologue family ", fams.links[i], "...\n" );
      next;
    }
    roots[[ length(roots) + 1 ]] <- .fix.names( root, quotes=quotes );
  }
  # Save them to file:
  attr(roots, "classif.name") <- "Ethnologue";
  class(roots) <- append("languageclassification", class(roots));
  export.languageclassification( roots, dir.name="../output", classification="ethnologue", export.nexus=EXPORT_NEXUS, nexus.translate.block=EXPORT_NEXUS_TRANSLATE_BLOCK, export.csv=EXPORT_CSV );
  
  rm(HTML.content); rm(fams.links);
}


######################
#
# Glottolog trees
#
######################

if( TRANSFORM_TREES )
{
  # The Glottolog trees: parse them and export them in the Family\tTree format, also updating the convention on the language codes:
  glott.fams <- readLines( "../input/glottolog/tree_glottolog_newick.txt", warn=FALSE );
  mapping <- code.mappings[,c("ISO","WALS","AUTOTYP","Glottolog")]
  roots <- NULL;
  for( i in 1:length(glott.fams) )
  {
    # Pre-process the tree by replacing any "\\'" within langauge names by "`", and taking care of names with ":" by replacing it with "-":
    gfam <- gsub( "\\'", "`", str_trim(glott.fams[i]), fixed=TRUE );
    tmp <- strsplit( gfam, "'", fixed=TRUE )[[1]];
    k <- grep( "[^:]+:", tmp );
    if( length(k) > 0 )
    {
      # ":" in names!
      for( j in k ) tmp[j] <- gsub( ":", "-", tmp[j], fixed=TRUE );
      gfam <- paste( tmp, collapse="'", sep="" );
    }
    
    # Make sure the Newick is enclosed between '(' and ');':
    if( substr(gfam,1,2) != "(" || substr(gfam,nchar(gfam)-1,nchar(gfam)) != ");" )
    {
      if( substr(gfam,nchar(gfam),nchar(gfam)) == ";" ) gfam <- paste0("(",substr(gfam,1,nchar(gfam)-1),");") else gfam <- paste0("(",gfam,");");
    }
    
    # Parse the tree:
    root <- familytree(gfam);
    # Set the family name (the root's only child):
    attr(root, "tree.name") <- str_trim(strsplit(strsplit(.get.node.name(root,root$edge[root$edge[,1] == find.root(root), 2][1]),"[",fixed=TRUE)[[1]][1],quotes,fixed=TRUE)[[1]][2]);
    # Convert the Glottolog convention to our conventions in what concerns the codes:
    root <- .convert.glottolog.convention(root, mapping);
    # Erase the branch length info:
    root <- compute.brlen(root, "none")$tree;
    
    # Fix the names and add it to the collection:
    roots[[ length(roots) + 1 ]] <- .fix.names(root, quotes);
  }
  # Save them to file:
  attr(roots, "classif.name") <- "Glottolog";
  class(roots) <- append("languageclassification", class(roots));
  export.languageclassification( roots, dir.name="../output", classification="glottolog", export.nexus=EXPORT_NEXUS, nexus.translate.block=EXPORT_NEXUS_TRANSLATE_BLOCK, export.csv=EXPORT_CSV );
  
  rm(glott.fams); rm(mapping);
}


######################
#
# AUTOTYP trees
#
######################

if( TRANSFORM_TREES && FALSE ) # I did not attempt to update the AUTOTYP trees so leave it untouched
{
  # AUTOTYP data:
  autotyp.data <- read.table( "../input/autotyp/autotyp-trees.csv", header=TRUE, sep="\t", quote="", stringsAsFactors=FALSE );
  
  # Replace troublesome characters in language, and branch names:
  autotyp.data$language <- sapply( as.character(autotyp.data$language), function( s ){ if( !is.na(s) ) { normalize.language.name(s); } else { NA; } } );
  autotyp.data$lsbranch <- sapply( as.character(autotyp.data$lsbranch), function( s ){ if( !is.na(s) ) { normalize.language.name(s); } else { NA; } } );
  autotyp.data$ssbranch <- sapply( as.character(autotyp.data$ssbranch), function( s ){ if( !is.na(s) ) { normalize.language.name(s); } else { NA; } } );
  autotyp.data$sbranch  <- sapply( as.character(autotyp.data$sbranch),  function( s ){ if( !is.na(s) ) { normalize.language.name(s); } else { NA; } } );
  autotyp.data$mbranch  <- sapply( as.character(autotyp.data$mbranch),  function( s ){ if( !is.na(s) ) { normalize.language.name(s); } else { NA; } } );
  autotyp.data$stock    <- sapply( as.character(autotyp.data$stock),    function( s ){ if( !is.na(s) ) { normalize.language.name(s); } else { NA; } } );
  
  # Replace the NAs in Glottolog and ISO codes by "":
  autotyp.data$glottolog_LID.2014[ is.na(autotyp.data$glottolog_LID.2014) ] <- "";
  autotyp.data$ISO639.3.2013[ is.na(autotyp.data$ISO639.3.2013) ] <- "";
  
  # Create the AUTOTYP trees:
  roots <- NULL;
  for( i in 1:nrow(autotyp.data) )
  {
    # Process this entry:
    if( !is.na(autotyp.data$stock[i]) )
    {
      # Build the path; the stock is always given:
      path <- create.name.and.codes(autotyp.data$stock[i], iso="", wals="", autotyp="", glottocode="", quotes=quotes); 
      # For the mbranch, sbranch, ssbranch & lsbranch things are a bit more complex in the sense that sometimes an intermediate level is skipped (NA) but deeper levels are given: add them with later tree-wide disambiguation:
      defined.levels <- !is.na(autotyp.data[ i, c( "mbranch", "sbranch", "ssbranch", "lsbranch" ) ]);
      if( any(defined.levels) )
      {
        # There is such a level defined
        deepest.defined.level <- max(which( !is.na(autotyp.data[ i, c( "mbranch", "sbranch", "ssbranch", "lsbranch" ) ]) ));
        mbranch.column <- which(names(autotyp.data)=="mbranch")-1;
        for( j in 1:deepest.defined.level )
        {
          if( is.na(autotyp.data[ i, j+mbranch.column ]) )
          {
            path <- c(path, create.name.and.codes(paste0("Group",j), iso="", wals="", autotyp="", glottocode="", quotes=quotes));
          } else
          {
            path <- c(path, create.name.and.codes(autotyp.data[ i, j+mbranch.column ], iso="", wals="", autotyp="", glottocode="", quotes=quotes));
          }
        }
      }
      # Add the language, also finding the other matching code(s):
      s <- which(code.mappings$AUTOTYP == as.character(autotyp.data$LID[i]));
      path <- c(path, create.name.and.codes(autotyp.data$language[i], 
                                            iso=as.character(autotyp.data$ISO639.3.2013[i]), 
                                            wals=ifelse( length(s) == 0, "", paste(code.mappings$WALS[s],collapse="-",sep="") ), 
                                            autotyp=as.character(autotyp.data$LID[i]),
                                            glottocode=as.character(autotyp.data$glottolog_LID.2014[i])));
      
      roots <- add.language.to.languageclassification( roots, 
                                                       path=path,
                                                       brlens=NULL );
    }
  }
  attr(roots, "classif.name") <- "autotyp";
  # Fix the family and language names (get rid of the codes):
  for( i in 1:length(roots) )
  {
    roots[[i]] <- .fix.names(roots[[i]], quotes=quotes);
    attr(roots[[i]], "tree.name") <- str_trim(strsplit(strsplit(get.name(roots[[i]]),quotes,fixed=TRUE)[[1]][2],"[",fixed=TRUE)[[1]][1]);
  }
  # Save them to file:
  export.languageclassification( roots, dir.name="../output", classification="autotyp", export.nexus=EXPORT_NEXUS, nexus.translate.block=EXPORT_NEXUS_TRANSLATE_BLOCK, export.csv=EXPORT_CSV );
  
  rm(autotyp.data);
}


#################################
#
# Compute the branch lengths
#
#################################

if( COMPUTE_BRLEN )
{
  # Apply various methods for computing branch length to a classification:
  compute.brlength.trees <- function( dir.name,                   # the directory where the original trees are and where the resulting trees will be saved
                                      classification,             # the classification: there must be a trees file in ./output/classification/classification-newick.csv
                                      export.nexus=TRUE,          # export the resulting trees to a NEXUS file?
                                      nexus.translate.block=TRUE, # when exporting to NEXUS, also generate a translate block?
                                      export.csv=TRUE,            # export the resulting trees to a CSV file?
                                      methods="all",              # the branch length methods to be used; can be "all", "none", or any subset of {"constant", "proportional", "grafen", "nnls", "ga", "nj"}
                                      constant=1.0,               # the positive constant required by the methods "constant" and "proportional"
                                      dists.codes=NULL,           # the data.frame of distance matrices to be used by "nnls", "ga" and "nj", and the codes used as the distance matrices' row and column names
                                      replace.NA.brlen.with=NA,   # some methods produce NA branch length (if the actual brlen cannot be estimated from the data): replace it by another numeric value?
                                      restore.collapsed.singles=TRUE, # some standard methods need to collapse singles, but we can restore them
                                      remove.all.missing.distances=TRUE, # remove languages with missing distance info?
                                      parallel.mc.cores=4,
                                      quotes="'"
  )
  {
    start.time <- Sys.time(); cat(paste0("Staring processing classification '",classification,"' @ ",start.time,"\n"));
    
    # Helper function: given a distances matrix with row and column names as language codes and which code to use, generate the corresponding submatrix with row and column full language names:
    dists.submatrix <- function( roots, dist.mat, code=c("iso","wals","autotyp","glottolog") )
    {
      # Sanity checks:
      if( is.null(roots) || length(roots) == 0 || is.null(dist.mat) || nrow(dist.mat) == 0 || ncol(dist.mat) == 0 )
      {
        return (NULL);
      }
      
      # Adapt the distance methods for these trees (keep only the languages in the intersection and generate the full names):
      lgs <- data.frame( "full"=unlist( lapply( roots, function(root){ extract.languages( root ); } ) ), "name"=NA, "iso"=NA, "wals"=NA, "autotyp"=NA, "glottolog"=NA, stringsAsFactors=FALSE );
      for( i in 1:nrow(lgs) )
      {
        tmp <- extract.name.and.codes( lgs$full[i], quotes );
        lgs$name[i]      <- normalize.language.name(tmp$name);
        lgs$iso[i]       <- paste0(tmp$iso,collapse="-");
        lgs$wals[i]      <- paste0(tmp$wals,collapse="-");
        lgs$autotyp[i]   <- paste0(tmp$autotyp,collapse="-");
        lgs$glottolog[i] <- paste0(tmp$glottolog,collapse="-");
      }
      
      # Expand the required code in the long format, keeping only the actually defined ones:
      codes <- lapply( 1:nrow(lgs), function(i){ tmp <- ifelse( lgs[i,code] == "", "", strsplit(lgs[i,code],"-",fixed=TRUE)[[1]] ); data.frame( "name"=lgs$full[i], "code"=tmp ); } );
      codes <- do.call( rbind, codes ); 
      codes <- codes[ codes$code != "", ];
      
      # Overlap santiy check:
      if( length( intersect( codes$code, rownames(dist.mat) ) ) == 0 )
      {
        return (NULL);
      }
      
      # The shared codes:
      shared.codes <- codes[ codes$code %in% rownames(dist.mat), ];
      if( is.null(shared.codes) || length(shared.codes) == 0 )
      {
        return (NULL);
      }
      # Select the submatrix corresponding to these codes
      dm <- dist.mat[ as.character(shared.codes$code), as.character(shared.codes$code), drop=FALSE ];
      if( is.null(dm) || !is.matrix(dm) || nrow(dm) == 0 || ncol(dm) == 0 )
      {
        return (NULL);
      }
      shared.codes <- shared.codes[ shared.codes$code %in% rownames(dm), , drop=FALSE ];
      if( is.null(shared.codes) || nrow(shared.codes) == 0 )
      {
        return (NULL);
      }
      full.language.names <- as.character(shared.codes$name);
      
      # For duplicated languages, keep the code that minimized the matrix-wide amount of missing data:
      duplgs <- duplicated( shared.codes$name ); 
      if( sum(duplgs,na.rm=TRUE) > 0 )
      {
        # The duplicated language names:
        duplgnames <- unique(shared.codes$name[ duplgs ]);
        # The tried combinations, each composed of an iso code per duplicated language and the resulting number of missing data in the distances matrix:
        lgcombs <- list();
        for( i in 1:length(duplgnames) )
        {
          lgcombs[[i]] <- shared.codes$code[ shared.codes$name==duplgnames[i] ];
        }
        lgcombs <- expand.grid( lgcombs );
        names(lgcombs) <- sapply( 1:length(duplgnames), function(i){ duplgnames[i] } );
        lgcombs <- cbind( lgcombs, "NoNAs"=NA );
        # The non-duplicated language codes:
        nondupcodes <- as.character(shared.codes$code[ !(shared.codes$name %in% duplgnames) ]);
        # Compute the number of NAs for each combination:
        for( i in 1:nrow(lgcombs) )
        {
          # The distances matrix:
          combcodes <- c( nondupcodes, as.character(unlist(lgcombs[i,1:ncol(lgcombs)-1])) );
          d <- dm[ combcodes, combcodes ];
          lgcombs$NoNAs[i] <- sum( is.na(d) );
        }
        if( sum( !is.na(lgcombs$NoNAs) ) > 0 )
        {
          # Pick one that minimizes the number of NAs:
          combcodes <- c( nondupcodes, as.character(unlist(lgcombs[which.min(lgcombs$NoNAs),1:ncol(lgcombs)-1])) );
          dm <- dm[ combcodes, combcodes ];
          # And the corresponding language names:
          full.language.names <- c( as.character(shared.codes$name[ !(shared.codes$name %in% duplgnames) ]), as.character(duplgnames) );
        }
      } 
      # Make the row and column names the full language names:
      rownames(dm) <- colnames(dm) <- full.language.names;
      
      # Return this submatrix:
      return (dm);
    }
    
    # Sanity checks:
    if( constant <= 0 )
    {
      cat( "The constant must be positive\n" );
      return (FALSE);
    }
    
    # Expand the methods to be used:
    if( length(methods) == 1 && methods == "none" )
    {
      # Nothing to do:
      return (TRUE);
    } else if( length(methods) == 1 && methods == "all" )
    {
      # All methods
      methods = c("constant", "proportional", "grafen", "nj", "nnls", "ga");
    } else if( !all( methods %in% c("constant", "proportional", "grafen", "nj", "nnls", "ga"), na.rm=TRUE ) )
    {
      cat( "Unknown methods\n" );
      return (FALSE);
    }   
    
    # Read the roots from file:
    cat( paste( "Reading the trees for classification '", classification, "'...\n", sep="") );
    roots <- languageclassification( classif.name=classification, csv.file=paste(dir.name,"/",classification,"/",classification,"-newick.csv",sep=""), csv.name.column="Family", csv.tree.column="Tree" );
    if( is.null(roots) || length(roots) == 0 || !inherits(roots,"languageclassification") )
    {
      return (FALSE);
    }
    
    # Compute the branch lengths and export them as separate files, one per method (by distance or constant):
    for( method in methods )
    {
      # Decide which constant and/or distances to use depending on the method (and store that in the d.c.k data.frame):
      if( method == "grafen" )
      {
        d.c.k <- data.frame("Distance"=NA, "ShortName"=NA, "File"=NA, "Code"=NA, "k"=NA);
      } else if( method %in% c("constant", "proportional") )
      {
        d.c.k <- data.frame("Distance"=NA, "ShortName"=NA, "File"=NA, "Code"=NA, "k"=constant);
      } else if( method %in% c("nnls", "ga", "nj") )
      {
        d.c.k <- cbind(dists.codes, "k"=NA);
      } else
      {
        # Unknown method; how did I get here?
        next;
      }
      
      # Apply the method:
      for( i in 1:nrow(d.c.k) )
      {
        # Print progress message:
        cat( paste( "\t\tProcesssing the '", classification, "' trees with method '", as.character(method), "'", 
                    ifelse( is.na(d.c.k$k[i]), "", paste( " with constant=", sprintf("%.02f",d.c.k$k[i]), sep="" ) ), 
                    ifelse( is.na(d.c.k$Distance[i]), "", paste( " with distance matrix '", d.c.k$Distance[i], "'", sep="" ) ),
                    "...\n", sep="" ) );
        
        # Load and adapt the distance matrix:
        if( !is.na(d.c.k$Distance[i]) )
        {
          load(d.c.k$File[i]); if( is.null(d) || !inherits(d, "matrix") || nrow(d) != ncol(d) || nrow(d) < 1 ){ cat("Distance ", d.c.k$Distance[i], " is malformed\n"); return (FALSE); }
          d <- dists.submatrix(roots, d,  d.c.k$Code[i]); 
          if( is.null(d) || !inherits(d, "matrix") || nrow(d) != ncol(d) || nrow(d) < 1 ){ cat("Distance ", d.c.k$Distance[i], " is empty for classification ", classification, "\n"); return (FALSE); }
        } else
        {
          d <- NULL;
        }
        
        # Apply the method:
        export.languageclassification( roots,
                                       dir.name=dir.name,
                                       classification=classification, 
                                       export.nexus=export.nexus,
                                       nexus.translate.block=nexus.translate.block,
                                       export.csv=export.csv,
                                       method=as.character(method), 
                                       constant=d.c.k$k[i], 
                                       distmatrix=d,
                                       replace.NA.brlen.with=replace.NA.brlen.with, 
                                       restore.collapsed.singles=restore.collapsed.singles,
                                       remove.all.missing.distances=remove.all.missing.distances,
                                       filename.postfix=paste( "-", method, 
                                                               ifelse( is.na(d.c.k$k[i]), "", paste( "=", sprintf("%.02f",d.c.k$k[i]), sep="") ), 
                                                               ifelse( is.na(d.c.k$Distance[i]), "", paste( "+", d.c.k$Distance[i], sep="") ), 
                                                               sep="" ), 
                                       quotes=quotes, 
                                       parallel.mc.cores=1 ); #parallel.mc.cores );
        
        # Free up the space:
        rm(d);
      }
    }
    
    end.time <- Sys.time(); cat(paste0("Ended processing classification '",classification,"' @ ",end.time," (took ",sprintf("%.1f minutes",difftime(end.time,start.time,units="mins")),")","\n"));
    
    return (TRUE);
  }
  compute.brlength.trees <- cmpfun(compute.brlength.trees);

  # Compute the branch lengths for each classification:
  x <- mclapply( CLASSIFICATIONS, function( classification )
  {
    # Apply only the relevant mg2015 distances:
    if( length(grep("mg2015", DISTS.CODES$Distance, fixed=TRUE)) > 0 )
    {
      dists.codes <- DISTS.CODES[ -grep("mg2015", DISTS.CODES$Distance, fixed=TRUE), ]; 
      dists.codes <- rbind(dists.codes, DISTS.CODES[ DISTS.CODES$Distance == paste0("mg2015(",classification,")"), ]);
    } else
    {
      dists.codes <- DISTS.CODES;
    }
    # Compute the branch lengths:
    compute.brlength.trees( dir.name="../output", classification=classification, 
                            export.nexus=EXPORT_NEXUS, nexus.translate.block=EXPORT_NEXUS_TRANSLATE_BLOCK, export.csv=EXPORT_CSV,
                            methods=METHODS, constant=CONSTANT, dists.codes=dists.codes, replace.NA.brlen.with=NA, restore.collapsed.singles=TRUE,
                            parallel.mc.cores=CPU_CORES, quotes=quotes );
  }, mc.cores=1); #length(CLASSIFICATIONS) );
}





