R/FMM_internal_restr.R

In alexARC26/FMM: Rhythmic Patterns Modeling by FMM Models

################################################################################
# Auxiliary functions for the fit of restricted FMM models
# Functions:
#   iterateOmegaGrid:   iterate fitting on omegasIter grid (FMM_unit_restr models
#                       with fixed omega)
#   backfittingRestr:   backfitting algorithm with fixed omegas.
#   fitFMM_unit_restr:  to fit monocomponent FMM models with fixed omega.
#   stepOmega:          to optimize omega.
#   step2FMM_restr:     second step of FMM fitting process with fixed omega
#   angularMean:        to compute the angular mean.
################################################################################

################################################################################
# Internal function: iterate fitting on omegasIter grid (FMM_unit_restr models
# with fixed omega)
# Arguments:
#   vData: data to be fitted an FMM model.
#   omegasIter: omegas grid (depending on original omega grid and on
#               restrictions)
#   betaRestrictions: beta's constraint vector.
#   timePoints: one single period time points.
#   lengthAlphaGrid: precision of the grid of alpha and omega parameters.
#   alphaGrid: grid of alpha parameter.
#   numReps: number of times the alpha-omega grid search is repeated.
#   nback: number of FMM components to fit.
#   maxiter: maximum number of iterations for the backfitting algorithm.
#   parallelize: TRUE to use parallelized procedure to fit FMM model
# Returns a list with objects of class FMM.
################################################################################
iterateOmegaGrid <- function(vData, omegasIter, betaRestrictions, timePoints = seqTimes(length(vData)),
                             lengthAlphaGrid = 48, alphaGrid = seq(0,2*pi,length.out = lengthAlphaGrid),
                             numReps = numReps, nback = nback, maxiter = maxiter, stopFunction = alwaysFalse,
                             parallelize = FALSE){
  if(parallelize){
    # cores for parallelized procedure registered: 'foreach' is used
    objectFMMList <- foreach::foreach(omegas = iterators::iter(omegasIter, by="row")) %dopar% {
      backfittingRestr(vData = vData, timePoints = timePoints, omegas = omegas,
                       lengthAlphaGrid = lengthAlphaGrid, alphaGrid = alphaGrid, numReps = numReps,
                       nback = nback, maxiter = maxiter, betaRestrictions = betaRestrictions,
                       stopFunction = stopFunction)
    }
  }else{
    # 'for' loop (non-parallelized version)
    objectFMMList <- list()
    for(i in 1:length(omegasIter[,1])){
      omegas <- omegasIter[i,]
      objectFMMList[[i]] <-  backfittingRestr(vData = vData, timePoints = timePoints,
                                              omegas = omegas, lengthAlphaGrid = lengthAlphaGrid,
                                              alphaGrid = alphaGrid, numReps = numReps, nback = nback,
                                              maxiter = maxiter, betaRestrictions = betaRestrictions,
                                              stopFunction = stopFunction)
    }
  }
  return(objectFMMList)

}


#######################################################################################
# Internal function: backfitting algorithm with fixed omegas.
# Arguments:
#   vData: data to be fitted an FMM model.
#   omegas: fixed omegas to fit FMM models.
#   betaRestrictions: beta's constraint vector.
#   timePoints: one single period time points.
#   alphaGrid: grid of alpha parameter.
#   numReps: number of times the alpha-omega grid search is repeated.
#   nback: number of FMM components to fit.
#   maxiter: maximum number of iterations for the backfitting algorithm.
#   parallelize: TRUE to use parallelized procedure to fit FMM model
#   stopFunction: Function to check the convergence criterion for the
#                 backfitting algorithm
# Returns an object of class FMM.
#######################################################################################
backfittingRestr <- function(vData, omegas, nback, betaRestrictions,
                             timePoints = seqTimes(length(vData)), lengthAlphaGrid = 48,
                             alphaGrid = seq(0,2*pi,length.out = lengthAlphaGrid),
                             numReps = 3, maxiter = nback, stopFunction = alwaysFalse){
  nObs <- length(vData)
  omegas <- as.numeric(omegas)
  fittedValuesPerComponent <- matrix(0, ncol = nback, nrow = nObs)
  fittedFMMPerComponent <- list()
  prevFittedFMMvalues <- NULL

  # Backfitting algorithm: iteration
  for(i in 1:maxiter){
    # Backfitting algorithm: component
    for(j in 1:nback){
      # data for component j: difference between vData and all other components fitted values
      backFittingData <- vData - apply(as.matrix(fittedValuesPerComponent[,-j]), 1, sum)
      # component j fitting using fitFMM_unit_restr function
      fittedFMMPerComponent[[j]] <- fitFMM_unit_restr(backFittingData, omegas[j], timePoints = timePoints,
                                                      lengthAlphaGrid = lengthAlphaGrid, 
                                                      #alphaGrid = alphaGrid[[j]],
                                                      numReps = numReps)
      fittedValuesPerComponent[,j] <- getFittedValues(fittedFMMPerComponent[[j]])
    }

    # Check stop criterion
    # Fitted values as sum of all components
    fittedFMMvalues <- apply(fittedValuesPerComponent, 1, sum)

    if(!is.null(prevFittedFMMvalues)){

      if(PV(vData, prevFittedFMMvalues) > PV(vData, fittedFMMvalues)){
        fittedFMMPerComponent <- previousFittedFMMPerComponent
        fittedFMMvalues <- prevFittedFMMvalues
        break
      }
      if(stopFunction(vData, fittedFMMvalues, prevFittedFMMvalues)){
        break
      }
    }
    prevFittedFMMvalues <- fittedFMMvalues
    previousFittedFMMPerComponent <- fittedFMMPerComponent
  }

  # alpha, beta y omega estimates
  alpha <- unlist(lapply(fittedFMMPerComponent, getAlpha))
  beta <- unlist(lapply(fittedFMMPerComponent, getBeta))
  omega <- unlist(lapply(fittedFMMPerComponent, getOmega))

  # beta restrictions: calculate angular mean of beta parameters
  # the nearest betas are chosen
  restBeta <- beta
  markedBetas <- rep(0, length(betaRestrictions))
  betaIndexVector <- 1:length(betaRestrictions)
  for(indRes in unique(betaRestrictions)){
    numComponents <- sum(betaRestrictions == indRes)
    betaIndex <- betaIndexVector[markedBetas == 0][1]
    distance <- abs(restBeta - restBeta[betaIndex])
    distance[markedBetas == 1] <- Inf
    nearestBetasIndex <- order(distance)[1:numComponents]
    markedBetas[nearestBetasIndex] <- 1
    restBeta[nearestBetasIndex] <- angularMean(restBeta[nearestBetasIndex])%%(2*pi)
  }

  # A and M estimates are recalculated by linear regression
  cosPhi <- calculateCosPhi(alpha = alpha, beta = restBeta, omega = omega, timePoints = timePoints)
  regresion <- lm(vData ~ cosPhi)

  M <- as.vector(coefficients(regresion)[1])
  A <- as.vector(coefficients(regresion)[-1])
  fittedFMMvalues <- predict(regresion)

  # Residual sum of squares
  SSE <- ifelse(sum(A < 0, na.rm = TRUE) > 0, Inf, sum((fittedFMMvalues-vData)^2))

  # Returns an object of class FMM
  return(FMM(
    M = M,
    A = A,
    alpha = alpha,
    beta = restBeta,
    omega = omega,
    timePoints = timePoints,
    summarizedData = vData,
    fittedValues = fittedFMMvalues,
    SSE = SSE,
    R2 = PVj(vData, timePoints, alpha, beta, omega),
    nIter = i
  ))
}

#######################################################################################
# Internal function: to fit monocomponent FMM models with fixed omega.
# Arguments:
#   vData: data to be fitted an FMM model.
#   omega: fixed value of the omega parameter.
#   timePoints: one single period time points.
#   lengthAlphaGrid: precision of the grid of alpha parameter.
#   alphaGrid: grid of alpha parameter.
#   numReps: number of times the alpha-omega grid search is repeated.
# Returns an object of class FMM.
#######################################################################################
fitFMM_unit_restr<-function(vData, omega, timePoints = seqTimes(length(vData)),
                            lengthAlphaGrid = 48, alphaGrid = seq(0, 2*pi, length.out = lengthAlphaGrid),
                            numReps = 3){
  usedApply = getApply(FALSE)[[1]]

  ## Step 1: initial values of M, A, alpha, beta and omega
  # alpha and omega are fixed and cosinor model is used to calculate the
  # rest of the parameters.
  # step1FMM function is used to make this estimate
  grid <- expand.grid(alphaGrid, omega)
  step1 <- usedApply(FUN = step1FMM, X = grid, vData = vData, timePoints = timePoints)
  colnames(step1) <- c("M","A","alpha","beta","omega","RSS")

  # We find the optimal initial parameters,
  # minimizing Residual Sum of Squared with several stability conditions.
  # We use bestStep1 internal function
  bestPar <- bestStep1(vData,step1)

  # When the value of the fixed omega is too extreme, making the fit impossible,
  # a null fitted model is returned.
  if(is.null(bestPar)){
    outMobius <- FMM(
      M = 0,
      A = 0,
      alpha = 0,
      beta = 0,
      omega = omega,
      timePoints = timePoints,
      summarizedData = vData,
      fittedValues = rep(0,length(vData)),
      SSE = sum(vData^2),
      R2 = PV(vData, rep(0,length(vData))),
      nIter = 0
    )
    return(outMobius)
  }

  ## Step 2: Nelder-Mead optimization. 'step2FMM_restr' function is used.
  if(!is.infinite(step2FMM_restr(bestPar[1:4], vData = vData, timePoints = timePoints, omega = omega))){
    nelderMead <- optim(par = bestPar[1:4], fn = step2FMM_restr, vData = vData,
                        timePoints = timePoints, omega = omega)
    parFinal <- c(nelderMead$par, omega)
  } else {
    parFinal <- c(bestPar[1:4],omega)
  }

  names(parFinal) <- c("M", "A", "alpha", "beta", "omega")
  fittedFMMvalues <- parFinal["M"] + parFinal["A"]*
    cos(parFinal["beta"] + 2*atan(parFinal["omega"]*tan((timePoints-parFinal["alpha"])/2)))
  SSE <- sum((fittedFMMvalues-vData)^2)

  # alpha and beta between 0 and 2pi
  parFinal[3] <- parFinal[3]%%(2*pi)
  parFinal[4] <- parFinal[4]%%(2*pi)

  # the grid search is repeated numReps
  numReps <- numReps - 1
  while(numReps > 0){
    # new grid for alpha between 0 and 2pi
    nAlphaGrid <- length(alphaGrid)
    amplitudeAlphaGrid <- 1.5*mean(diff(alphaGrid))
    alphaGrid <- seq(parFinal[3]-amplitudeAlphaGrid,parFinal[3]+amplitudeAlphaGrid,
                     length.out = nAlphaGrid)
    alphaGrid <- alphaGrid%%(2*pi)

    ## Step 1: initial parameters
    grid <- as.matrix(expand.grid(alphaGrid,omega))
    step1 <- usedApply(FUN = step1FMM, X = grid, vData = vData, timePoints = timePoints)
    colnames(step1) <- c("M","A","alpha","beta","omega","RSS")
    previousBestPar <- bestPar
    bestPar <- bestStep1(vData,step1)

    # None satisfies the conditions
    if(is.null(bestPar)){
      bestPar <- previousBestPar
      numReps <- 0
    }

    ## Step 2: Nelder-Mead optimization
    if(!is.infinite(step2FMM_restr(bestPar[1:4],vData = vData, timePoints = timePoints, omega = omega))){
      nelderMead <- optim(par = bestPar[1:4], fn = step2FMM_restr, vData = vData,
                          timePoints = timePoints, omega = omega)
      parFinal <- c(nelderMead$par,omega)
    } else {
      parFinal <- c(bestPar[1:4],omega)
    }

    # alpha and beta between 0 and 2pi
    parFinal[3] <- parFinal[3]%%(2*pi)
    parFinal[4] <- parFinal[4]%%(2*pi)

    numReps <- numReps - 1
  }
  names(parFinal) <- c("M","A","alpha","beta","omega")

  # Returns an object of class FMM
  fittedFMMvalues <- parFinal["M"] + parFinal["A"]*
    cos(parFinal["beta"] + 2*atan(parFinal["omega"]*tan((timePoints-parFinal["alpha"])/2)))
  SSE <- sum((fittedFMMvalues-vData)^2)

  return(FMM(
    M = parFinal["M"],
    A = parFinal["A"],
    alpha = parFinal[3],
    beta = parFinal[4],
    omega = parFinal[5],
    timePoints = timePoints,
    summarizedData = vData,
    fittedValues = fittedFMMvalues,
    SSE = SSE,
    R2 = 0,
    nIter = 0
  ))
}

################################################################################
# Internal function: to optimize omega in the extra optimization step of omega,
# within fitFMM_restr function.
# Arguments:
#   uniqueOmegas: grid of omega parameters.
#   indOmegas: omegas' constraint vector.
#   objFMM: FMM object to refine the omega fitting.
#   omegaMax: max value for omega.
################################################################################
stepOmega <- function(uniqueOmegas, indOmegas, objFMM, omegaMax){

  timePoints <- getTimePoints(objFMM)
  M <- getM(objFMM)
  A <- getA(objFMM)
  alpha <- getAlpha(objFMM)
  beta <- getBeta(objFMM)
  omega <- uniqueOmegas[indOmegas]
  vData <- getSummarizedData(objFMM)

  # FMM fitting and RSS
  fittedValues <- A%*%t(calculateCosPhi(alpha = alpha, beta = beta, omega = omega, timePoints = timePoints)) + M
  RSS <- sum((fittedValues - vData)^2)

  # If the integrity conditions are valid, it returns RSS
  # else it returns infinite.
  rest1 <- all(uniqueOmegas <= omegaMax)
  rest2 <- all(uniqueOmegas >= 0)
  if(rest1 & rest2){
    return(RSS)
  }else{
    return(Inf)
  }
}

################################################################################
# Internal function: second step of FMM fitting process with fixed omega.
# Arguments:
#   param: M, A, alpha, beta initial parameter estimations.
#   vData: data to be fitted an FMM model.
#   timePoints: one single period time points.
#   omega: fixed value of omega.
################################################################################
step2FMM_restr <- function(parameters, vData, timePoints, omega){

  nObs <- length(timePoints)
  # FMM model and residual sum of squares
  modelFMM <- parameters[1] + parameters[2] *
    cos(parameters[4]+2*atan2(omega*sin((timePoints - parameters[3])/2),
                              cos((timePoints - parameters[3])/2)))
  residualSS <- sum((modelFMM - vData)^2)/nObs
  sigma <- sqrt(residualSS*nObs/(nObs - 5))

  # When amplitude condition is valid, it returns RSS
  # else it returns infinite.
  amplitudeUpperBound <- parameters[1] + parameters[2]
  amplitudeLowerBound <- parameters[1] - parameters[2]
  rest1 <- amplitudeUpperBound <= max(vData) + 1.96*sigma
  rest2 <- amplitudeLowerBound >= min(vData) - 1.96*sigma

  # Other integrity conditions that must be met
  rest3 <- parameters[2] > 0 #A > 0
  if(rest1 & rest2 & rest3){
    return(residualSS)
  }else{
    return(Inf)
  }
}

################################################################################
# Internal function: to compute the angular mean.
# Arguments:
#   angles: input vector of angles.
################################################################################
angularMean <- function(angles){
  return(atan2(sum(sin(angles)), sum(cos(angles))))
}