diff ReadFids_script.R @ 7:122df1bf0a8c draft default tip

planemo upload for repository https://github.com/workflow4metabolomics/nmr_preprocessing commit 3d328007dd7716848ec2eeb6c2a472f27eeb2995

--- a/ReadFids_script.R	Wed May 13 03:52:09 2020 -0400
+++ b/ReadFids_script.R	Fri Jul 11 08:33:38 2025 +0000
@@ -2,403 +2,380 @@
 #
 #   Read FIDs in Bruker format
 #
-#
 ################################################################################################
-
-
-# vec2mat ==============================================================================
+# Function to convert a vector to a matrix
 vec2mat <- function(vec) {
-  return(matrix(vec, nrow = 1, dimnames = list(c(1), names(vec)))) 
+    return(matrix(vec, nrow = 1, dimnames = list(c(1), names(vec))))
 }
 
-
-# ReadFid ==============================================================================
+# Function to read FID files
+# Format of files depends on version of software used for spectra aquisition (< or > 4.x.x)
+# Read 1D FID using Bruker XWinNMR and TopSpin format.  It is inspired of the
+# matNMR matlab library which deals with 2D FID and also other formats
+# Based on aquisition parameters stored in the acqus file
 ReadFid <- function(path) {
-  
-  # Read 1D FID using Bruker XWinNMR and TopSpin format.  It is inspired of the
-  # matNMR matlab library which deals with 2D FID and also other formats
-  # Read also the parameters in the acqus file
-  
-  paramFile <- file.path(path, "acqus")
-  # BYTEORDA: 0 -> Little Endian 1 -> Big Endian
-  params <- readParams(paramFile, c("TD", "BYTORDA", "DIGMOD", "DECIM", "DSPFVS", 
-                                    "SW_h", "SW", "O1"))
-  
-  if (params[["DSPFVS"]] >= 20) {
+    paramFile <- file.path(path, "acqus")
+    params <- readParams(paramFile, c("TD", "BYTORDA", "DIGMOD", "DECIM", "DSPFVS", "SW_h", "SW", "O1"))
+
+    # Version of the TopSpin software
+    line <- readLines(paramFile)[1]
+    version <- as.numeric(substr(strsplit(line, " ")[[1]][5], 1, 1))
+    neo <- (version > 3)
+
     # The group delay first order phase correction is given directly from version 20
-    grpdly <- readParams(paramFile, c("GRPDLY"))
-    params[["GRPDLY"]] <- grpdly[["GRPDLY"]]
-  }
-  TD <- params[["TD"]]
-  endianness <- if (params$BYTORDA) 
-    "big" else "little"
-  if (TD%%2 != 0) {
-    stop(paste("Only even numbers are allowed for size in TD because it is complex 
-               data with the real and imaginary part for each element.", 
-               "The TD value is in the", paramFile, "file"))
-  }
-  
-  # Interpret params Dwell Time, time between 2 data points in the FID
-  params[["DT"]] <- 1/(2 * params[["SW_h"]])
-  
-  # Read fid
-  fidFile <- file.path(path, "fid")
-  fidOnDisk <- readBin(fidFile, what = "int", n = TD, size = 4L, endian = endianness)
-  
-  # Real size that is on disk (it should be equal to TD2, except for TopSpin/Bruker
-  # (which is our case) according to matNMR as just discussed
-  TDOnDisk <- length(fidOnDisk)
-  if (TDOnDisk < TD) {
-    warning("Size is smaller than expected, the rest is filled with zero so the size is the same for every fid")
-    fidGoodSize <- sapply(vector("list", length = TD), function(x) 0)
-    fidGoodSize[1:TDOnDisk] <- fidOnDisk
-    
-  } else if (TDOnDisk > TD) {
-    warning("Size is bigger than expected, the rest ignored so the size is the same for every fid")
-    fidGoodSize <- fidOnDisk(1:TD)
-    
-  } else {
-    fidGoodSize <- fidOnDisk
-  }
-  
-  fidRePart <- fidGoodSize[seq(from = 1, to = TD, by = 2)]
-  fidImPart <- fidGoodSize[seq(from = 2, to = TD, by = 2)]
-  fid <- complex(real = fidRePart, imaginary = fidImPart)
-  
-  return(list(fid = fid, params = params))
-}
+    if (params[["DSPFVS"]] >= 20) {
+        grpdly <- readParams(paramFile, c("GRPDLY"))
+        params[["GRPDLY"]] <- grpdly[["GRPDLY"]]
+    }
+    TD <- params[["TD"]]
+
+    endianness <- if (params$BYTORDA) {
+        "big"
+    } else {
+        "little"
+    }
+    if (TD %% 2 != 0) {
+        stop(paste("Only even numbers are allowed for size in TD because it is complex data with the real and imaginary part for each element.", "The TD value is in the", paramFile, "file"))
+    }
+
+    # Interpret params Dwell Time, time between 2 data points in the FID
+    params[["DT"]] <- 1 / (2 * params[["SW_h"]])
 
-
-
+    # Read fid depend on version of aquisition software
+    fidFile <- file.path(path, "fid")
+    if (neo) {
+        fidOnDisk <- readBin(fidFile, what = "double", n = TD, size = NA_integer_, signed = TRUE, endian = .Platform$endian)
+    } else {
+        fidOnDisk <- readBin(fidFile, what = "int", n = TD, size = 4L, endian = endianness)
+    }
 
-# getDirsContainingFid ==============================================================================
-getDirsContainingFid <- function(path) {
-  subdirs <- dir(path, full.names = TRUE)
-  if (length(subdirs) > 0) {
-    cond <- sapply(subdirs, function(x)  {
-      content <- dir(x)
-      # subdirs must contain fid, acqu and acqus files
-      return("fid" %in% content && "acqu" %in% content && "acqus" %in% content)
-    })
-    subdirs <- subdirs[cond]
-  }
-  return(subdirs)
+    # Real size that is on disk (it should be equal to TD2, except for TopSpin/Bruker
+    # (which is our case) according to matNMR as just discussed
+    TDOnDisk <- length(fidOnDisk)
+    if (TDOnDisk < TD) {
+        warning("Size is smaller than expected, the rest is filled with zero so the size is the same for every fid")
+        fidGoodSize <- sapply(vector("list", length = TD), function(x) 0)
+        fidGoodSize[1:TDOnDisk] <- fidOnDisk
+    } else if (TDOnDisk > TD) {
+        warning("Size is bigger than expected, the rest ignored so the size is the same for every fid")
+        fidGoodSize <- fidOnDisk(1:TD)
+    } else {
+        fidGoodSize <- fidOnDisk
+    }
+
+    fidRePart <- fidGoodSize[seq(from = 1, to = TD, by = 2)]
+    fidImPart <- fidGoodSize[seq(from = 2, to = TD, by = 2)]
+    fid <- complex(real = fidRePart, imaginary = fidImPart)
+
+    return(list(fid = fid, params = params))
 }
 
-
-
-
-
-
-# beginTreatment ==============================================================================
-
-beginTreatment <- function(name, Signal_data = NULL, Signal_info = NULL, 
-                           force.real = FALSE) {
-  
-  cat("Begin", name, "\n")
-  
-  
-  # Formatting the Signal_data and Signal_info -----------------------
-  
-  vec <- is.vector(Signal_data)
-  if (vec) {
-    Signal_data <- vec2mat(Signal_data)
-  }
-  if (is.vector(Signal_info)) {
-    Signal_info <- vec2mat(Signal_info)
-  }
-  if (!is.null(Signal_data)) {
-    if (!is.matrix(Signal_data)) {
-      stop("Signal_data is not a matrix.")
+# Function to obtain path to files
+getDirsContainingFid <- function(path) {
+    subdirs <- dir(path, full.names = TRUE)
+    if (length(subdirs) > 0) {
+        cond <- sapply(subdirs, function(x) {
+            content <- dir(x)
+            # subdirs must contain fid, acqu and acqus files
+            return("fid" %in% content && "acqu" %in% content && "acqus" %in% content)
+        })
+        subdirs <- subdirs[cond]
     }
-    if (!is.complex(Signal_data) && !is.numeric(Signal_data)) {
-      stop("Signal_data contains non-numerical values.")
-    }
-  }
-  if (!is.null(Signal_info) && !is.matrix(Signal_info)) {
-    stop("Signal_info is not a matrix.")
-  }
-  
-  
-  Original_data <- Signal_data
-  
-  # Extract the real part of the spectrum ---------------------------
-  
-  if (force.real) {
-    if (is.complex(Signal_data)) {
-      Signal_data <- Re(Signal_data)
-    } else {
-      # The signal is numeric Im(Signal_data) is zero anyway so let's avoid
-      # using complex(real=...,imaginary=0) which would give a complex signal
-      # in endTreatment()
-      force.real <- FALSE
-    }
-  }
-  
-  
-  # Return the formatted data and metadata entries --------------------
-  
-  return(list(start = proc.time(), vec = vec, force.real = force.real, 
-              Original_data = Original_data, Signal_data = Signal_data, Signal_info = Signal_info))
+    return(subdirs)
 }
 
-
-# endTreatment ==============================================================================
+# Function to xxx
+beginTreatment <- function(name, Signal_data = NULL, Signal_info = NULL, force.real = FALSE) {
+    cat("Begin", name, "\n")
 
-endTreatment <- function(name, begin_info, Signal_data) {
-  end_time = proc.time() # record it as soon as possible
-  start_time = begin_info[["start"]]
-  delta_time = end_time - start_time
-  delta = delta_time[]
-  cat("End", name, "\n")
-  cat("It lasted",
-      round(delta["user.self"], 3), "s user time,",
-      round(delta["sys.self"] , 3), "s system time and",
-      round(delta["elapsed"]  , 3), "s elapsed time.\n")
-  if (begin_info[["force.real"]]) {
-    # The imaginary part is left untouched
-    i <- complex(real=0, imaginary=1)
-    Signal_data = Signal_data + i * Im(begin_info[["Original_data"]])
-  }
-  if (begin_info[["vec"]]) {
-    Signal_data = Signal_data[1,]
-  }
-  return(Signal_data)
+    # Formatting the Signal_data and Signal_info
+    vec <- is.vector(Signal_data)
+    if (vec) {
+        Signal_data <- vec2mat(Signal_data)
+    }
+    if (is.vector(Signal_info)) {
+        Signal_info <- vec2mat(Signal_info)
+    }
+    if (!is.null(Signal_data)) {
+        if (!is.matrix(Signal_data)) {
+            stop("Signal_data is not a matrix.")
+        }
+        if (is.vector(Signal_info)) {
+            Signal_info <- vec2mat(Signal_info)
+        }
+    }
+    if (!is.null(Signal_info) && !is.matrix(Signal_info)) {
+        stop("Signal_info is not a matrix.")
+    }
+
+    Original_data <- Signal_data
+
+    # Extract the real part of the spectrum
+    if (force.real) {
+        if (is.complex(Signal_data)) {
+            Signal_data <- Re(Signal_data)
+        } else {
+            # The signal is numeric Im(Signal_data) is zero anyway so let's avoid
+            # using complex(real=...,imaginary=0) which would give a complex signal
+            # in endTreatment()
+            force.real <- FALSE
+        }
+    }
+
+    # Return the formatted data and metadata entries
+    return(list(start = proc.time(), vec = vec, force.real = force.real, Original_data = Original_data, Signal_data = Signal_data, Signal_info = Signal_info))
 }
 
-
-# checkArg ==============================================================================
+# Function to get information on reading end
+endTreatment <- function(name, begin_info, Signal_data) {
+    end_time <- proc.time() # record it as soon as possible
+    start_time <- begin_info[["start"]]
+    delta_time <- end_time - start_time
+    delta <- delta_time[]
+    cat("End", name, "\n")
+    cat("It lasted", round(delta["user.self"], 3), "s user time,", round(delta["sys.self"], 3), "s system time and", round(delta["elapsed"], 3), "s elapsed time.\n")
+    if (begin_info[["force.real"]]) {
+        # The imaginary part is left untouched
+        i <- complex(real = 0, imaginary = 1)
+        Signal_data <- Signal_data + i * Im(begin_info[["Original_data"]])
+    }
+    if (begin_info[["vec"]]) {
+        Signal_data <- Signal_data[1, ]
+    }
+    return(Signal_data)
+}
 
-checkArg <- function(arg, checks, can.be.null=FALSE) {
-  check.list <- list(bool=c(is.logical, "a boolean"),
-                     int =c(function(x){x%%1==0}, "an integer"),
-                     num =c(is.numeric, "a numeric"),
-                     str =c(is.character, "a string"),
-                     pos =c(function(x){x>0}, "positive"),
-                     pos0=c(function(x){x>=0}, "positive or zero"),
-                     l1 =c(function(x){length(x)==1}, "of length 1")
-  )
-  if (is.null(arg)) {
-    if (!can.be.null) {
-      stop(deparse(substitute(arg)), " is null.")
+# Function to check arguments
+checkArg <- function(arg, checks, can.be.null = FALSE) {
+    check.list <- list(bool = c(is.logical, "a boolean"), int = c(function(x) {
+        x %% 1 == 0
+    }, "an integer"), num = c(is.numeric, "a numeric"), str = c(is.character, "a string"), pos = c(function(x) {
+        x > 0
+    }, "positive"), pos0 = c(function(x) {
+        x >= 0
+    }, "positive or zero"), l1 = c(function(x) {
+        length(x) == 1
+    }, "of length 1"))
+    if (is.null(arg)) {
+        if (!can.be.null) {
+            stop(deparse(substitute(arg)), " is null.")
+        }
     }
-  } else {
-    if (is.matrix(arg)) {
-      stop(deparse(substitute(arg)), " is not scalar.")
-    }
-    for (c in checks) {
-      if (!check.list[[c]][[1]](arg)) {
-        stop(deparse(substitute(arg)), " is not ", check.list[[c]][[2]], ".")
-      }
-    }
-  }
 }
 
 
-# getArg ==============================================================================
-
-getArg <- function(arg, info, argname, can.be.absent=FALSE) {
-  if (is.null(arg)) {
-    start <- paste("impossible to get argument", argname, "it was not given directly and");
-    if (!is.matrix(info)) {
-      stop(paste(start, "the info matrix was not given"))
+# Function to get arguments
+getArg <- function(arg, info, argname, can.be.absent = FALSE) {
+    if (is.null(arg)) {
+        start <- paste("impossible to get argument", argname, "it was not given directly and")
+        if (!is.matrix(info)) {
+            stop(paste(start, "the info matrix was not given"))
+        }
+        if (!(argname %in% colnames(info))) {
+            if (can.be.absent) {
+                return(NULL)
+            } else {
+                stop(paste(start, "is not in the info matrix"))
+            }
+        }
+        if (nrow(info) < 1) {
+            stop(paste(start, "the info matrix has no row"))
+        }
+        arg <- info[1, argname]
+        if (is.na(arg)) {
+            stop(paste(start, "it is NA in the info matrix"))
+        }
     }
-    if (!(argname %in% colnames(info))) {
-      if (can.be.absent) {
-        return(NULL)
-      } else {
-        stop(paste(start, "is not in the info matrix"))
-      }
-    }
-    if (nrow(info) < 1) {
-      stop(paste(start, "the info matrix has no row"))
-    }
-    arg <- info[1,argname]
-    if (is.na(arg)) {
-      stop(paste(start, "it is NA in the info matrix"))
-    }
-  }
-  return(arg)
+    return(arg)
 }
 
-
-
-# getTitle ==============================================================================
-
-# Get the name of the signal from the title file or fromt the name of the subdirectory
-# Get the name of the signal from the title file or fromt the name of the subdirectory
-
+# Function to get title and name samples
+# Get the name of the signal from the title file or from the name of the subdirectory
 getTitle <- function(path, l, subdirs) {
-  title <- NULL
-  title_file <- file.path(file.path(file.path(path, "pdata"), "1"), "title")
-  if (file.exists(title_file)) {
-    lines <- readLines(title_file, warn = FALSE)
-    if (length(lines) >= 1)  {
-      first_line <- gsub("^\\s+|\\s+$", "", lines[l])
-      if (nchar(first_line) >= 1)  {
-        title <- first_line
-      } else {
-        warning(paste("The", l ,"line of the title file is blank for directory ", 
-                      path, "and the (sub)dirs names are used instead"))
-      }
+    title <- NULL
+    title_file <- file.path(file.path(file.path(path, "pdata"), "1"), "title")
+    if (file.exists(title_file)) {
+        lines <- readLines(title_file, warn = FALSE)
+        if (length(lines) >= 1) {
+            first_line <- gsub("^\\s+|\\s+$", "", lines[l])
+            if (nchar(first_line) >= 1) {
+                title <- first_line
+            } else {
+                warning(paste("The", l, "line of the title file is blank for directory ", path, "and the (sub)dirs names are used instead"))
+            }
+        } else {
+            warning(paste("Title file doesn't exists for directory ", path, "\n the (sub)dirs names are  used instead"))
+        }
     } else {
-      warning(paste("The title file is empty for directory ", path, "and the (sub)dirs names are used instead"))
+        warning(paste("Title file doesn't exists for directory ", path, "\n the (sub)dirs names are  used instead"))
     }
-  } else {
-    warning(paste("Title file doesn't exists for directory ", path, "\n the (sub)dirs names are  used instead"))
-  }
-  if (is.null(title)) {
-    if(subdirs) {
-      separator <- .Platform$file.sep
-      path_elem <- strsplit(path,separator)[[1]]
-      title <- paste(path_elem[length(path_elem)-1], path_elem[length(path_elem)], sep = "_")
-    } else{title <- basename(path)} 
-  }
-  return(title)
+    if (is.null(title)) {
+        if (subdirs) {
+            separator <- .Platform$file.sep
+            path_elem <- strsplit(path, separator)[[1]]
+            title <- paste(path_elem[length(path_elem) - 1], path_elem[length(path_elem)], sep = "_")
+        } else {
+            title <- basename(path)
+        }
+    }
+    return(title)
 }
 
-
-
-
-# readParams ==============================================================================
-# Read parameter values for Fid_info in the ReadFids function
-
+# Function to read parameter values for Fid_info in the ReadFids function
 readParams <- function(file, paramsName) {
-  
-  isDigit <- function(c) {
-    return(suppressWarnings(!is.na(as.numeric(c))))
-  }
-  lines <- readLines(file)
-  params <- sapply(paramsName, function(x) NULL)
-  
-  for (paramName in paramsName)  {
-    # Find the line with the parameter I add a '$' '=' in the pattern so that for
-    # example 'TD0' is not found where I look for 'TD' and LOCSW and WBSW when I look
-    # for 'SW'
-    pattern <- paste("\\$", paramName, "=", sep = "")
-    occurences <- grep(pattern, lines)
-    if (length(occurences) == 0L)  {
-      stop(paste(file, "has no field", pattern))
+    isDigit <- function(c) {
+        return(suppressWarnings(!is.na(as.numeric(c))))
     }
-    if (length(occurences) > 1L) {
-      warning(paste(file, "has more that one field", pattern, " I take the first one"))
+    lines <- readLines(file)
+    params <- sapply(paramsName, function(x) NULL)
+
+    for (paramName in paramsName) {
+        # Find the line with the parameter I add a '$' '=' in the pattern so that for
+        # example 'TD0' is not found where I look for 'TD' and LOCSW and WBSW when I look
+        # for 'SW'
+        pattern <- paste("\\$", paramName, "=", sep = "")
+        occurences <- grep(pattern, lines)
+        if (length(occurences) == 0L) {
+            stop(paste(file, "has no field", pattern))
+        }
+        lines <- readLines(file)
+        params <- sapply(paramsName, function(x) NULL)
+
+        for (paramName in paramsName) {
+            # Find the line with the parameter I add a '$' '=' in the pattern so that for
+            # example 'TD0' is not found where I look for 'TD' and LOCSW and WBSW when I look
+            # for 'SW'
+            pattern <- paste("\\$", paramName, "=", sep = "")
+            occurences <- grep(pattern, lines)
+            if (length(occurences) == 0L) {
+                stop(paste(file, "has no field", pattern))
+            }
+            if (length(occurences) > 1L) {
+                warning(paste(file, "has more that one field", pattern, " I take the first one"))
+            }
+            line <- lines[occurences[1]]
+
+            # Cut beginning and end of the line '##$TD= 65536' -> '65536'
+            igual <- as.numeric(regexpr("=", line))
+
+            first <- igual
+            while (first <= nchar(line) & !isDigit(substr(line, first, first))) {
+                first <- first + 1
+            }
+            last <- nchar(line)
+            while (last > 0 & !isDigit(substr(line, last, last))) {
+                last <- last - 1
+            }
+            params[paramName] <- as.numeric(substr(line, first, last))
+        }
+        line <- lines[occurences[1]]
+
+        # Cut beginning and end of the line '##$TD= 65536' -> '65536'
+        igual <- as.numeric(regexpr("=", line))
+
+        first <- igual
+        while (first <= nchar(line) & !isDigit(substr(line, first, first))) {
+            first <- first + 1
+        }
+        last <- nchar(line)
+        while (last > 0 & !isDigit(substr(line, last, last))) {
+            last <- last - 1
+        }
+        params[paramName] <- as.numeric(substr(line, first, last))
     }
-    line <- lines[occurences[1]]
-    
-    # Cut beginning and end of the line '##$TD= 65536' -> '65536'
-    igual = as.numeric(regexpr("=", line))
-    
-    first <- igual
-    while (first <= nchar(line) & !isDigit(substr(line, first, first))) {
-      first <- first + 1
-    }
-    last <- nchar(line)
-    while (last > 0 & !isDigit(substr(line, last, last)))  {
-      last <- last - 1
-    }
-    params[paramName] <- as.numeric(substr(line, first, last))
-  }
-  return(params)
+    return(params)
 }
 
-
-
-# ReadFids ==============================================================================
-
+# Function to read all fid's in the directory
 ReadFids <- function(path, l = 1, subdirs = FALSE, dirs.names = FALSE) {
-  
-  # Data initialisation and checks ----------------------------------------------
-  begin_info <- beginTreatment("ReadFids")
-  checkArg(path, c("str"))
-  checkArg(l, c("pos"))
-  if (file.exists(path) == FALSE) {
-    stop(paste("Invalid path:", path))
-  }
-  
-  
-  # Extract the FIDs and their info ----------------------------------------------
-  
-  if (subdirs == FALSE) {
-    fidDirs <- getDirsContainingFid(path)
-    n <- length(fidDirs)
-    if (n == 0L)  {
-      stop(paste("No valid fid in", path))
-    }
-    if (dirs.names) {
-      separator <- .Platform$file.sep
-      path_elem <- strsplit(fidDirs,separator)
-      fidNames <- sapply(path_elem, function(x) x[[length(path_elem[[1]])]])
-    }else {fidNames <- sapply(X = fidDirs, FUN = getTitle, l = l, subdirs = subdirs,  USE.NAMES = F)}
-    
-    for (i in 1:n)  {
-      fidList <- ReadFid(fidDirs[i])
-      fid <- fidList[["fid"]]
-      info <- fidList[["params"]]
-      m <- length(fid)
-      if (i == 1)  {
-        Fid_data <- matrix(nrow = n, ncol = m, dimnames = list(fidNames, 
-                                                               info[["DT"]] * (0:(m - 1))))
-        Fid_info <- matrix(nrow = n, ncol = length(info), dimnames = list(fidNames, 
-                                                                          names(info)))
-      }
-      Fid_data[i, ] <- fid
-      Fid_info[i, ] <- unlist(info)
+    # Data initialisation and checks
+    begin_info <- beginTreatment("ReadFids")
+    checkArg(path, c("str"))
+    checkArg(l, c("pos"))
+    if (file.exists(path) == FALSE) {
+        stop(paste("Invalid path:", path))
     }
-    
-  } else  {
-    maindirs <- dir(path, full.names = TRUE) # subdirectories
-    Fid_data <- numeric()
-    Fid_info <- numeric()
-    
-    fidDirs <- c()
-    for (j in maindirs) {
-      fd <- getDirsContainingFid(j) # recoved FIDs from subdirectories
-      n <- length(fd)
-      if (n > 0L)  {
-        fidDirs <- c(fidDirs, fd)
-      } else {warning(paste("No valid fid in",j ))}
+
+    # Extract the FIDs and their info
+    if (subdirs == FALSE) {
+        fidDirs <- getDirsContainingFid(path)
+        n <- length(fidDirs)
+        if (n == 0L) {
+            stop(paste("No valid fid in", path))
+        }
+        if (dirs.names) {
+            separator <- .Platform$file.sep
+            path_elem <- strsplit(fidDirs, separator)
+            fidNames <- sapply(path_elem, function(x) x[[length(path_elem[[1]])]])
+        } else {
+            fidNames <- sapply(X = fidDirs, FUN = getTitle, l = l, subdirs = subdirs, USE.NAMES = F)
+        }
+
+        for (i in 1:n) {
+            fidList <- ReadFid(fidDirs[i])
+            fid <- fidList[["fid"]]
+            info <- fidList[["params"]]
+            m <- length(fid)
+            if (i == 1) {
+                Fid_data <- matrix(nrow = n, ncol = m, dimnames = list(fidNames, info[["DT"]] * (0:(m - 1))))
+                Fid_info <- matrix(nrow = n, ncol = length(info), dimnames = list(fidNames, names(info)))
+            }
+            Fid_data[i, ] <- fid
+            Fid_info[i, ] <- unlist(info)
+        }
+    } else {
+        maindirs <- dir(path, full.names = TRUE) # subdirectories
+        Fid_data <- numeric()
+        Fid_info <- numeric()
+
+        fidDirs <- c()
+        for (j in maindirs) {
+            fd <- getDirsContainingFid(j) # recoved FIDs from subdirectories
+            n <- length(fd)
+            if (n > 0L) {
+                fidDirs <- c(fidDirs, fd)
+            } else {
+                warning(paste("No valid fid in", j))
+            }
+        }
+
+        if (dirs.names == TRUE) {
+            if (length(fidDirs) != length(dir(path))) { # at least one subdir contains more than 1 FID
+                separator <- .Platform$file.sep
+                path_elem <- strsplit(fidDirs, separator)
+                fidNames <- sapply(path_elem, function(x) paste(x[[length(path_elem[[1]]) - 1]], x[[length(path_elem[[1]])]], sep = "_"))
+            } else {
+                fidNames <- dir(path)
+            }
+        } else {
+            fidNames <- sapply(X = fidDirs, FUN = getTitle, l = l, subdirs = subdirs, USE.NAMES = F)
+        }
+
+        for (i in 1:length(fidNames)) {
+            fidList <- ReadFid(fidDirs[i])
+            fid <- fidList[["fid"]]
+            info <- fidList[["params"]]
+            m <- length(fid)
+            if (i == 1) {
+                Fid_data <- matrix(nrow = length(fidNames), ncol = m, dimnames = list(fidNames, info[["DT"]] * (0:(m - 1))))
+                Fid_info <- matrix(nrow = length(fidNames), ncol = length(info), dimnames = list(fidNames, names(info)))
+            }
+            print(paste("i=", i, "Fid_data=", ncol(Fid_data), "fid=", length(fid)))
+
+            Fid_data[i, ] <- fid
+            Fid_info[i, ] <- unlist(info)
+        }
     }
-    
-    if (dirs.names==TRUE) {
-      if (length(fidDirs)!= length(dir(path))) { # at least one subdir contains more than 1 FID
-        separator <- .Platform$file.sep
-        path_elem <- strsplit(fidDirs,separator)
-        fidNames <- sapply(path_elem, function(x) paste(x[[length(path_elem[[1]])-1]],
-                                                        x[[length(path_elem[[1]])]], sep = "_"))
-      }else {fidNames <- dir(path)}
-      
-    } else {fidNames <- sapply(X = fidDirs, FUN = getTitle, l = l, subdirs = subdirs, USE.NAMES = F)}
-    
-    for (i in 1:length(fidNames))  {
-      fidList <- ReadFid(fidDirs[i])
-      fid <- fidList[["fid"]]
-      info <- fidList[["params"]]
-      m <- length(fid)
-      if (i == 1)  {
-        Fid_data <- matrix(nrow = length(fidNames), ncol = m, dimnames = list(fidNames, 
-                                                                              info[["DT"]] * (0:(m - 1))))
-        Fid_info <- matrix(nrow = length(fidNames), ncol = length(info), dimnames = list(fidNames, 
-                                                                                         names(info)))
-      }
-      Fid_data[i, ] <- fid
-      Fid_info[i, ] <- unlist(info)
+
+    # Check for non-unique IDs
+    NonnuniqueIds <- sum(duplicated(row.names(Fid_data)))
+    cat("dim Fid_data: ", dim(Fid_data), "\n")
+    cat("IDs: ", rownames(Fid_data), "\n")
+    cat("non-unique IDs?", NonnuniqueIds, "\n")
+    if (NonnuniqueIds > 0) {
+        warning("There are duplicated IDs: ", Fid_data[duplicated(Fid_data)])
     }
-    
-    
-  }
-  
-  # Check for non-unique IDs ----------------------------------------------
-  NonnuniqueIds <- sum(duplicated(row.names(Fid_data)))
-  cat("dim Fid_data: ", dim(Fid_data), "\n")
-  cat("IDs: ", rownames(Fid_data), "\n")
-  cat("non-unique IDs?", NonnuniqueIds, "\n")
-  if (NonnuniqueIds > 0) {
-    warning("There are duplicated IDs: ", Fid_data[duplicated(Fid_data)])
-  }
-  
-  
-  # Return the results ----------------------------------------------
-  return(list(Fid_data = endTreatment("ReadFids", begin_info, Fid_data), Fid_info = Fid_info))
-  
+
+    # Return the results
+    return(list(Fid_data = endTreatment("ReadFids", begin_info, Fid_data), Fid_info = Fid_info))
 }
-