R/MainFunctions.R

In FixedPoint: Algorithms for Finding Fixed Point Vectors of Functions

#' A function for finding the fixed point of a contraction mapping
#'
#' This function takes in a function and an initial guess for the fixed point of that function. It then uses a fixed point method to find the
#' fixed point of that function.
#' @export
#' @param Function This is the function for which a fixed point is sought. This function must take and return a vector of the same dimension.
#' @param Inputs This can be either a vector of values that is an initial guess for a fixed point or it can be an N x A matrix of previous inputs
#' for which corresponding outputs are available. In this case N is the dimensionality of the fixed point vector you are seeking (Hence each
#' column is a matrix that is input to the "Function") and A is the number of previous Inputs/Outputs that are being provided to the fixed point.
#' Where a matrix is input, a corresponding outputs must be provided or the last column of the outputs matrix is taken as a startpoint guess and
#' the rest of the inputs and output matrices are discarded.
#' @param Outputs (Optional) This is a matrix of the Function values for each column of the input. It must be provided so that column k of the
#' outputs matrix is equal to Function(Column k of inputs matrix).
#' @param Method This is the fixed point method to be used. It can be "Anderson", "Simple", "Aitken", "Newton", "MPE", "RRE", "VEA" or "SEA". See
#' vignette and references to see explanations of these methods.
#' @param ConvergenceMetric This is a function that takes in a vector of residuals from one iterate of the function (defined as f(x) - x for
#' vector x and function f) and returns a scalar. This scalar should be low when convergence is close to being achieved. By default this is
#' the maximum residual by absolute value (the sup norm in the space of residuals).
#' @param ConvergenceMetricThreshold This is the threshold for terminating the algorithm. The algorithm will terminate when the scalar that
#' ConvergenceMetric returns is less than ConvergenceMetricThreshold. This can be set to a negative number in which case the algorithm will
#' run until MaxIter is hit or an error occurs (Note that an error is likely in trying to use any method other than "Simple" when a fixed point
#' is already found).
#' @param MaxIter This is the maximum number of iterates that will be undertaken.
#' @param MaxM This is the maximum number of saved iterates that are used in the Anderson algorithm. It has no effect if another method is chosen.
#' Note that the number of previous iterates that will actually be used is the minimum of MaxIter, the dimensionality of the function's vector and
#' the number of inputs that have been tried to previously (the width of the Outputs matrix at each given stage of the algorithm). If
#' PrintReports = TRUE, the number of previous iterates actually used is reported as the algorithm is running.
#' @param ExtrapolationPeriod This is the number of simple iterates to perform before extrapolating. This is used for the MPE, RRE, VEA and SEA
#' methods and has no effect if another method is chosen. Where an epsilon algorithm is used this should be one plus a multiple of two, ie (3,5,7,etc).
#' @param Dampening This is the dampening parameter. By default it is 1 which means no dampening takes place. It can also be less than 1
#' (indicating dampening) or more than 1 (indicating extrapolation).
#' @param ReplaceInvalids This determines how NAs and Infs generated by the extrapolation method are handled. If it is "NoAction" then if
#' the extrapolation algorithm generates an NA or Inf then no action is taken. This will likely result in these NAs and Infs being fed into
#' the function and an error resulting. If it is "ReplaceVector" then the entire extrapolated vector will be replaced by a vector of a simple
#' iterate. If it is "ReplaceElements" then the elements containing an NA or Inf will be replaced by the corresponding elements from a simple iterate.
#' @param PrintReports This is a boolean describing whether to print ongoing ConvergenceMetric values for each iterate.
#' @param ReportingSigFig This is the number of significant figures that will be used in printing the convergence values to the console
#' (only if PrintReports is TRUE).
#' @param ConditionNumberThreshold This is a threshold for what condition number is acceptable for solving the least squares problem for the
#' Anderson Method. If the condition number is larger than this threshold then fewer previous iterates will be used in solving the problem.
#' This has no effect unless the "Anderson" method is used.
#' @param Plot This determines whether a plot should be drawn for every iterate. It can be "NoPlot", "ConvergenceFig" or "ChangePerIterate".
#' By default it is "NoPlot" and no plot is drawn. If it is "ConvergenceFig" then a plot is shown with iterates on the x axis and convergence
#' (as defined by the ConvergenceMetric) is on the y axis. If it is "ChangePerIterate" then there is the index of the array value on the x axis
#' and the value of the array value on the y axis. The previous iterate is also shown so the change per iterate can be visualised.
#' @param ConvergenceFigLags This only affects anything if Plot == "ConvergenceFig". This gives how many previous iterates should be shown on
#' the x axis. By default it is 5. To see them all set it to a high number.
#' @param ChangePerIteratexaxis This only affects anything if Plot == "ChangePerIterate". Sometimes there is a more appropriate xaxis value
#' to use than (the default) value index for this figure. For instance in the consumption smoothing problem in the vignette every value is a
#' value function value at a given budget level. In this case the budget levels could be used for this xaxis.
#' @param DropOldIterates This drops old iterates that will not be used in determining the next fixedpoint guess under the current setting.
#' @return A list containing the FixedPoint, the Inputs and corresponding Outputs, and convergence values (which are computed under the
#' "ConvergenceMetric"). The list will also include a "Finish" statement describing why it has finished. This is often going to be due to
#' either MaxIter or ConvergenceMetricThreshold being  reached. It may also terminate due to an error in generating a new input guess or using
#' the function with that guess. If this occurs the function will terminate early and the "Finish" statement will describe the issue.
#' In this event there will also be additional objects returned in the list "NewInputVector" and possibly "NewOutputVector" for use in debugging
#' the issue.
#' @examples
#' # For the simplest possible example we can seek the fixed point of the cos function with a scalar.
#' Inputs = 0.3
#' Function = function(x){ cos(x) }
#' A = FixedPoint(Function, Inputs, Method = "Aitken", Dampening = 0.5)
#' B = FixedPoint(Function, Inputs, Method = "Anderson", Dampening = 1.0)
#'
#' # For this next one the ConvergenceMetricThreshold is negative so the algorithm
#' # will keep running until MaxIter is met.
#' C = FixedPoint(Function, Inputs, Method = "Simple", MaxIter = 4, ConvergenceMetricThreshold = -1)
#' # This is not yet close to the fixed point but we can continue solving for this fixed point
#' # by inputting the previous inputs and outputs to the function. We can also switch methods
#' # and will switch below to the Newton Method.
#' D = FixedPoint(Function, C$Inputs, C$Outputs, Method = "Newton")
#'
#' # We can also find a 4 dimensional fixed point vector of this function.
#' Inputs = c(0.3, 98, 0, pi)
#' E = FixedPoint(Function, Inputs, Method = "Anderson")
#' F = FixedPoint(Function, Inputs, Method = "Anderson", MaxM = 4, ReportingSigFig = 13)
FixedPoint = function(Function, Inputs, Outputs = c(), Method = c("Anderson", "Simple", "Aitken", "Newton", "MPE", "RRE", "VEA", "SEA") ,
                      ConvergenceMetric  = function(Residuals){max(abs(Residuals))} , ConvergenceMetricThreshold = 1e-10, MaxIter = 1e3, MaxM = 10,
                      ExtrapolationPeriod = 7, Dampening = 1, ReplaceInvalids = c("ReplaceElements","ReplaceVector","NoAction"), PrintReports = FALSE,
                      ReportingSigFig = 5, ConditionNumberThreshold = 1e3, Plot = c("NoPlot", "ConvergenceFig", "ChangePerIterate"),
                      ConvergenceFigLags = 5, ChangePerIteratexaxis = c(), DropOldIterates = FALSE){
  # Checking inputs to the function
  Method = Method[1]
  Plot   = Plot[1]
  if (ConditionNumberThreshold < 1){stop("ConditionNumberThreshold must be at least 1.")}
  if (is.vector(Inputs)){
    Inputs  = matrix(Inputs, ncol = 1)
    Outputs = NULL
  }
  if (is.null(Outputs)){
    if (dim(Inputs)[2] > 1.5){
      warning("If you do not give outputs to the function then you can only give one vector of inputs. As you have input a matrix of input values
              everything but the last column has been discarded.\n")
      Inputs = Inputs[,dim(Inputs)[2]]
      }
  } else {
    if (sum(dim(Inputs) != dim(Outputs)) > 0.5) {
      warning("If you input a matrix of outputs as well as a matrix of inputs then inputs and outputs must be the same shape. As they differ in
               this case the last column of the inputs matrix has been taken as the starting point and everything else discarded.\n")
      Inputs  = Inputs[,dim(Inputs)[2]]
      Outputs = NULL
    }
  }
  SimpleStartIndex = dim(Inputs)[2]

  # Do an initial run if no runs have been done:
  if (is.null(Outputs)){
    Inputs          = matrix(Inputs, ncol = 1)
    NewOutputVector = try(Function(Inputs[,1]))
    # In the event of an error we return informative information on the error.
    if (inherits(NewOutputVector, "try-error")){return(list(Inputs = Inputs, Outputs = NA, Convergence = NA, FixedPoint = NA,
                                                           Finish = "Could not execute function with input vector", fpevals = 1))}
    if (sum(is.na(NewOutputVector))){return(list(Inputs = Inputs, Outputs = NA, Convergence = NA, FixedPoint = NA,
                                                 Finish = "New output vector contains NAs", NewInputVector = Inputs,
                                                 NewOutputVector = NewOutputVector, fpevals = 1))}
    # If no error we put the results into the Outputs matrix.
    Outputs = matrix(Function(Inputs[,1]), ncol = 1)
  } else {
    # This ensures that MaxIter refers to max iter excluding any previous passed in results
    MaxIter = MaxIter + dim(Outputs)[2]
    # This is to take into account previously passed in simple iterations (without jumps).
    SimpleStartIndex = SimpleStartIndex- (dim(PutTogetherIteratesWithoutJumps(Inputs, Outputs))[2] )
  }

  # First running through the last column of Inputs to test if we already have a fixed point.
  Resid = Outputs - Inputs
  iter = dim(Resid)[2]
  ConvergenceVector = sapply(1:iter, function(x) ConvergenceMetric(Resid[,x]) )
  if (ConvergenceVector[iter] < ConvergenceMetricThreshold){
    if (PrintReports){cat("The last column of Inputs matrix is already a fixed point under input convergence metric and convergence threshold")}
    return(list(FixedPoint = Outputs[,iter] , Inputs = Inputs, Outputs = Outputs, Convergence = ConvergenceVector, fpevals = iter))
  }
  # Printing a report for initial convergence
  Convergence = ConvergenceVector[iter]
  if (PrintReports){cat(paste0(format(" ", width = 49, justify = "right"), "Method: ", format(Method, width = 8, justify = "right")   ,
                               ". Iteration: ", format(iter, digits =  5, width = 5,scientific = FALSE),
                               ". Convergence: ",formatC(Convergence,  big.interval = 1, small.interval = ReportingSigFig, format = "e"),"\n"))}
  if (Plot == "ConvergenceFig"){ConvergenceFig(Inputs, Outputs,  Input_Convergence = ConvergenceVector, FromIterate = 1)}
  if (Plot == "ChangePerIterate"){ChangePerIterate(Inputs, Outputs, ConvergenceVector, FromIterate = dim(Inputs)[2],
                                                   ToIterate = dim(Inputs)[2], xaxis = ChangePerIteratexaxis, secondhold = -1)}
  iter = iter + 1

  while ((Convergence > ConvergenceMetricThreshold) & (iter <= MaxIter)){
    # Generating new input guess
    NewInputVector = try({FixedPointNewInput(Inputs = Inputs, Outputs = Outputs, Method = Method, MaxM = MaxM,
                                              SimpleStartIndex = SimpleStartIndex, ExtrapolationPeriod = ExtrapolationPeriod,
                                              Dampening = Dampening, ConditionNumberThreshold = ConditionNumberThreshold,
                                              PrintReports = PrintReports, ReplaceInvalids = ReplaceInvalids)})
    if (inherits(NewInputVector, "try-error")){return(list(Inputs = Inputs, Outputs = Outputs, Convergence = ConvergenceVector,
                                                          FixedPoint = NA, Finish = "Could not generate new input vector", fpevals = iter))}
    if (sum(is.na(NewInputVector))>0.5){return(list(Inputs = Inputs, Outputs = Outputs, Convergence = ConvergenceVector,
                                                FixedPoint = NA, Finish = "New input vector contains NAs", NewInputVector = NewInputVector, fpevals = iter))}
    if (sum(is.infinite(NewInputVector))>0.5){return(list(Inputs = Inputs, Outputs = Outputs, Convergence = ConvergenceVector,
                                                    FixedPoint = NA, Finish = "New input vector contains Infs", NewInputVector = NewInputVector, fpevals = iter))}
    if (Method != "Anderson" & PrintReports){cat(paste0(format(" ", width = 49, justify = "right")))}
    # Putting this input through the function to get a new output guess.
    NewOutputVector = try({Function(NewInputVector)})
    if (inherits(NewOutputVector, "try-error")){return(list(Inputs = Inputs, Outputs = Outputs, Convergence = ConvergenceVector,
                                                           FixedPoint = NA, Finish = "Could not execute function with generated vector",
                                                           NewInputVector = NewInputVector, fpevals = iter))}
    if (sum(is.na(NewOutputVector))){return(list(Inputs = Inputs, Outputs = Outputs, Convergence = ConvergenceVector,
                                                 FixedPoint = NA, Finish = "New output vector contains NAs", NewInputVector = NewInputVector,
                                                 NewOutputVector = NewOutputVector, fpevals = iter))}
    # Adding input and output to saved matrices.
    width = length(ConvergenceVector)
    adjustedData = adjustData(Inputs, Outputs, Resid, ConvergenceVector, NewInputVector, NewOutputVector, SimpleStartIndex, ConvergenceMetric, width, iter, MaxM, ExtrapolationPeriod, Method, DropOldIterates)
    Inputs = adjustedData$Inputs
    Outputs = adjustedData$Outputs
    Resid = adjustedData$Resids
    ConvergenceVector = adjustedData$ConvergenceVector
    SimpleStartIndex = adjustedData$SimpleStartIndex
    newWidth = length(ConvergenceVector)
    # Checking and recording convergence
    Convergence = ConvergenceVector[newWidth]
    # Output of report and going to next iteration.
    if (PrintReports){cat(paste0("Method: ", format(Method, width = 8, justify = "right")   , ". Iteration: ",
                                 format(iter, digits =  5, width = 5,scientific = FALSE), ". Convergence: ",
                                 formatC(Convergence,  big.interval = 1, small.interval = ReportingSigFig, format = "e"), "\n"))}
    if (Plot == "ConvergenceFig"){ConvergenceFig(Inputs, Outputs,  Input_Convergence =  ConvergenceVector,
                                                 FromIterate = max(1, dim(Inputs)[2]- ConvergenceFigLags))}
    if (Plot == "ChangePerIterate"){ChangePerIterate(Inputs, Outputs, ConvergenceVector, secondhold = -1,
                                                     FromIterate = dim(Inputs)[2], ToIterate = dim(Inputs)[2],
                                                     xaxis = ChangePerIteratexaxis)}
    iter  = iter + 1
  }
  # Extracting FixedPoint
  FixedPoint = Outputs[,dim(Outputs)[2]]
  # Returning result
  if (Convergence < ConvergenceMetricThreshold){Finish = "Reached Convergence Threshold"} else {Finish = "Reached MaxIter"}
  return(list(Inputs = Inputs, Outputs = Outputs, Convergence = ConvergenceVector, FixedPoint = FixedPoint, Finish = Finish, fpevals = iter-1))
}

adjustData = function(oldInputs, oldOutputs, oldResids, oldConvergences, newInput, newOutput, SimpleStartIndex, ConvergenceMetric, width, iter, MaxM, ExtrapolationPeriod, Method, DropOldIterates){
  newResid = newOutput - newInput
  newConvergence = ConvergenceMetric(newResid)
  expandedWidth = width+1
  Inputs   = matrix(c(oldInputs , newInput)           , ncol = expandedWidth, byrow = FALSE)
  Outputs  = matrix(c(oldOutputs, newOutput)          , ncol = expandedWidth, byrow = FALSE)
  Resid    = matrix(c(oldResids , newResid - newInput), ncol = expandedWidth, byrow = FALSE)
  ConvergenceVector =  c(oldConvergences, newConvergence)
  if (!(DropOldIterates)){
    return(list(Inputs = Inputs, Outputs = Outputs, Resids = Resid, ConvergenceVector = ConvergenceVector, SimpleStartIndex = SimpleStartIndex))
  }
  from = 1
  if (Method == "Simple"){
    from = expandedWidth
  } else if (Method == "Anderson"){
    from = max(1,expandedWidth - (MaxM+1))
  } else if (Method %in% c("Aitken", "Newton")){
    from = max(1, expandedWidth - 4)
  } else if (Method %in% c("MPE", "RRE", "VEA", "SEA")){
    from = max(1, expandedWidth - (ExtrapolationPeriod+1))
  }
  newWidth = expandedWidth - from + 1
  iteratesLost = (newWidth - width)-1
  newSimpleStartIndex = SimpleStartIndex
  if (Method == "Simple"){
    newSimpleStartIndex = 1
  } else if (Method == "Anderson"){
    newSimpleStartIndex = newWidth
  } else if (Method %in% c("Aitken", "Newton")){
    startIndices = c(SimpleStartIndex - iteratesLost, SimpleStartIndex + 3)
    newSimpleStartIndex = min(startIndices[startIndices > 0])
  } else if (Method %in% c("MPE", "RRE", "VEA", "SEA")){
    startIndices = c(SimpleStartIndex - iteratesLost, SimpleStartIndex + ExtrapolationPeriod)
    newSimpleStartIndex = min(startIndices[startIndices > 0])
  }

  return(list(Inputs = matrix(Inputs[,from:expandedWidth], ncol = newWidth), Outputs           = matrix(Outputs[,from:expandedWidth]         , ncol = newWidth),
              Resids = matrix(Resid[,from:expandedWidth] , ncol = newWidth), ConvergenceVector = matrix(ConvergenceVector[from:expandedWidth], ncol = newWidth),
              SimpleStartIndex = newSimpleStartIndex))
}

#' FixedPointNewInput
#' This function takes the previous inputs and outputs from the FixedPoint function and determines what vector to try next in seeking
#' a fixed point.
#' @export
#' @param Inputs This is an N x A matrix of previous inputs for which corresponding outputs are available. In this case N is the
#' dimensionality of the fixed point vector that is being sought (Hence each column is a matrix that is input to the "Function") and A
#' is the number of previous Inputs/Outputs that are being provided to the fixed point.
#' @param Outputs This is a matrix of Function values for the each column of the "Inputs" matrix.
#' @param Method This is the fixed point method to be used. It can be "Anderson", "Simple", "Aitken", "Newton", "MPE", "RRE", "VEA", "SEA".
#' @param MaxM This is the number of saved iterates that are used in the Anderson algorithm. This has no role if another method is used.
#' @param SimpleStartIndex This is the index for what column in the input/output matrices did the algorithm start doing simple iterates
#' without jumps. This is used for all methods except the simple and Anderson methods where it has no effect.
#' @param ExtrapolationPeriod This is the number of simple iterates to perform before extrapolating. This is used for the MPE, RRE, VEA and SEA
#' methods and has no effect if another method is chosen. Where an epsilon algorithm is used this should be one plus a multiple of two, ie (3,5,7,etc).
#' @param Dampening This is the dampening parameter. By default it is 1 which means no dampening takes place. It can also be less than 1
#' (indicating dampening) or more than 1 (indicating extrapolation).
#' @param ConditionNumberThreshold This is what threshold should be chosen to drop previous iterates if the matrix is ill conditioned.
#' Only used in Anderson acceleration.
#' @param PrintReports This is a boolean describing whether to print ongoing ConvergenceMetric values for each iterate.
#' @param ReplaceInvalids This determines how NAs and Infs generated by the extrapolation method are handled. If it is "NoAction" then if
#' the extrapolation algorithm generates an NA or Inf then no action is taken. This will likely result in these NAs and Infs being fed into
#' the function and an error resulting. If it is "ReplaceVector" then the entire extrapolated vector will be replaced by a vector of a simple
#' iterate. If it is "ReplaceElements" then the elements containing an NA or Inf will be replaced by the corresponding elements from a simple iterate.
#' @return A vector containing the next guess in seeking a fixed point.
#' @examples
#' FPFunction = function(x){c(0.5*sqrt(abs(x[1] + x[2])), 1.5*x[1] + 0.5*x[2])}
#' A = FixedPoint( Function = FPFunction, Inputs = c(0.3,900), MaxIter = 6, Method = "Simple")
#' NewGuessAnderson = FixedPointNewInput(A$Inputs, A$Outputs, Method = "Anderson")
#' NewGuessVEA = FixedPointNewInput(A$Inputs, A$Outputs, Method = "VEA")
#' NewGuessMPE = FixedPointNewInput(A$Inputs, A$Outputs, Method = "MPE")
#' NewGuessAitken = FixedPointNewInput(A$Inputs, A$Outputs, Method = "Aitken")
FixedPointNewInput = function(Inputs, Outputs, Method = "Anderson", MaxM = 10, SimpleStartIndex = 1,
                               ExtrapolationPeriod = 7, Dampening = 1, ConditionNumberThreshold = 1e3,
                              PrintReports = FALSE, ReplaceInvalids = c("ReplaceElements","ReplaceVector","NoAction") ) {
  CompletedIters = dim(Outputs)[2]
  ReplaceInvalids = ReplaceInvalids[1]

  if (Method == "Simple"){
    return( as.vector(Outputs[,CompletedIters] ))
  } else if (Method == "Anderson"){
    if (CompletedIters < 1.5){
      if (PrintReports){cat(paste0(format(" ", width = 32, justify = "right"), "  Using",  format(3, width = 3, justify = "right")," lags. "))}
      return(as.vector( Outputs[,CompletedIters] ))
    }
    # Extract variables for anderson acceleration algorithm.
    VectorLength   = dim(Outputs)[1]
    M = min(MaxM-1,CompletedIters-1,VectorLength)
    Coeffs = NA
    Truncated_Outputs = matrix(Outputs[, seq(CompletedIters-M, CompletedIters,1)], ncol = M+1)
    Truncated_Inputs  = matrix(Inputs[ , seq(CompletedIters-M, CompletedIters,1)], ncol = M+1)
    Resid             = Truncated_Outputs - Truncated_Inputs
    DeltaOutputs      = matrix(Truncated_Outputs[,2:(M+1)] - Truncated_Outputs[,1:M], ncol = M)
    DeltaResids       = matrix(Resid[  ,2:(M+1)] - Resid[  ,1:M], ncol = M)
    LastResid       = Resid[,M+1]
    LastOutput      = Truncated_Outputs[,M+1]
    # We nest the fitting of coefficients in a while loop. This allows us to drop columns and try again if we get NAs from a lack of numerical stability.
    while (sum(is.na(Coeffs))>0.5){
      # For numerical stability reasons we monitor the condition number
      ConditionNumber = rcond(DeltaResids)
      if (ConditionNumber > ConditionNumberThreshold){
        M = M-1
        DeltaOutputs= matrix(DeltaOutputs[, 2:(M+1)], ncol = M)
        DeltaResids = matrix(DeltaResids[,2:(M+1)], ncol = M)
        next
      }
      # We fit a linear model to get coefficients.
      Fit             = stats::lm.fit(x = DeltaResids, y = LastResid, tol = 1e-13)
      Coeffs          = Fit$coefficients
      # If we have NAs then we drop a previous iterate and will try again in the next loop.
      if (sum(is.na(Coeffs))>0.5){
      M = M-1
      if (M < 1.5){
        # This happens occasionally in test cases where the iteration is very close to a fixed point.
        if (PrintReports){cat(paste0(format(" ", width = 32, justify = "right"), "  Using",
                                     format(3, width = 3, justify = "right")," lags. "))}
        return(LastOutput)
      }
      DeltaOutputs= matrix(DeltaOutputs[, 2:(M+1)], ncol = M)
      DeltaResids = matrix(DeltaResids[,2:(M+1)], ncol = M)
      }
    }
    # Print information on the use of the anderson method and return the new guess.
    if (PrintReports){cat(paste0("Condition number is ", format(format(ConditionNumber, digits = 5,scientific = TRUE), width = 12, justify = "right")
                               ,". Used:",  format(M+1, width = 3, justify = "right"), " lags. "))}
      NewGuess = LastOutput - Coeffs  %*% t(DeltaOutputs)
      NewGuess = as.vector(NewGuess)
  } else if (Method == "Aitken"){
    if ((CompletedIters + SimpleStartIndex) %% 3 == 0){
      # If we are in 3rd, 6th, 9th, 12th iterate from when we started Acceleration then we want to do a jumped Iterate,
      # First we extract the guess that started this run of 3 iterates (x), the Function applied to it (fx) and the function applied to that (ffx)
      x = Inputs[,CompletedIters -1]
      fx = Outputs[,CompletedIters -1]
      ffx = Outputs[,CompletedIters ]
      # Now using the appropriate formula to make a new guess. Note that if a vector is input here it is used elementwise.
      NewGuess = x - ((fx - x)^2/(ffx - 2*fx + x))
    } else {
      # We just do a simple iterate. We do an attempt with the latest iterate.
      return( as.vector(Outputs[,CompletedIters] ))
      }
  } else if (Method == "Newton"){
    if ((CompletedIters + SimpleStartIndex) %% 2 == 1 & CompletedIters > 1){
      # If we are in 3rd, 6th, 9th, 12th iterate from when we started Acceleration then we want to do a jumped Iterate,
      # First we extract the guess that started this run of 3 iterates (xk1), the Function applied to it (fxk1, xk) and the function applied to that (fxk)
      xk1 = Inputs[,CompletedIters -1]
      fxk1 = Outputs[,CompletedIters -1]
      gxk1 = fxk1 - xk1
      xk  = Inputs[,CompletedIters ]
      fxk = Outputs[,CompletedIters ]
      gxk = fxk - xk
      # Now using the appropriate formula to make a new guess. Note that if a vector is input here it is used elementwise.
      derivative = (gxk-gxk1)/(xk - xk1)
      NewGuess = xk - (gxk /derivative)
    } else {
      # We just do a simple iterate. We do an attempt with the latest iterate.
      return( as.vector(Outputs[,CompletedIters] ))
    }
  } else if (Method == "MPE" | Method == "RRE"){
    SimpleIteratesMatrix = PutTogetherIteratesWithoutJumps(Inputs, Outputs)
    if ( dim(SimpleIteratesMatrix)[2] %% ExtrapolationPeriod == 0  ){
      NewGuess = PolynomialExtrapolation(SimpleIteratesMatrix, Method = Method)
    } else {
      # We just do a simple iterate.
      return(as.vector(Outputs[,CompletedIters]))
    }
  } else if (Method == "VEA" | Method == "SEA"){
    SimpleIteratesMatrix = PutTogetherIteratesWithoutJumps(Inputs, Outputs)
    if ( dim(SimpleIteratesMatrix)[2] %% ExtrapolationPeriod == 0  ){
      NewGuess = EpsilonExtrapolation(SimpleIteratesMatrix, Method = Method)
    } else {
      # We just do a simple iterate.
      return(as.vector(Outputs[,CompletedIters] ))
    }
  }
  # Here we use the ReplaceInvalids input to remove NAs and Infs.
  if (ReplaceInvalids == "ReplaceElements"){
    if (sum(is.infinite(NewGuess)) + sum(is.na(NewGuess)) > 0.5){
      NewGuess[(is.infinite(NewGuess) | is.na(NewGuess))] = Outputs[(is.infinite(NewGuess) | is.na(NewGuess)) ,CompletedIters]
  } else if (ReplaceInvalids == "ReplaceVector"){
    if (sum(is.infinite(NewGuess)) + sum(is.na(NewGuess)) > 0.5){
      NewGuess = Outputs[,CompletedIters]
     }
  }
  }
  # Here we apply any dampening and return the new guess.
  return(as.vector(Dampening*NewGuess + (1-Dampening)*Outputs[,CompletedIters] ))
}


#' PolynomialExtrapolation
#' This function performs Minimal Polynomial extrapolation (MPE) or Reduced Rank Extrapolation (RRE) given a matrix of previous iterates of the
#' function.
#' @export
#' @import MASS
#' @param Iterates A matrix of inputs and outputs excluding jumps. Can be pieced together from Inputs and Outputs matrices of the FixedPoint
#' function using the PutTogetherIteratesWithoutJumps function.
#' @param Method The method for polynomial extrapolation. Should be either "MPE" for minimal polynomial extrapolation or "RRE" for reduced
#' rank extrapolation.
#' @return A vector containing the extrapolated vector.
#' @examples
#' FPFunction = function(x){c(0.5*sqrt(abs(x[1] + x[2])), 1.5*x[1] + 0.5*x[2])}
#' A = FixedPoint( Function = FPFunction, Inputs = c(0.3,900), MaxIter = 6, Method = "Simple")
#' Iterates = PutTogetherIteratesWithoutJumps(A$Inputs, A$Outputs)
#' PolynomialExtrapolation(Iterates, "MPE")
#' B = FixedPoint( function(x){cos(x)}, Inputs = 1, MaxIter = 5, Method = "Simple")
#' Iterates = PutTogetherIteratesWithoutJumps(B$Inputs, B$Outputs)
#' PolynomialExtrapolation(Iterates, "RRE")
PolynomialExtrapolation = function(Iterates, Method = c("MPE", "RRE")){
  Method = Method[1]
  if (!(Method %in% c("MPE", "RRE"))){stop("Invalid method input. PolynomialExtrapolation function can only take Method as MPE or RRE.")}
  if (Method == "MPE"){
    TotalColumnsOfIterates = dim(Iterates)[2]
    Differences          = matrix(Iterates[,2:(TotalColumnsOfIterates-1)] - Iterates[,1:(TotalColumnsOfIterates-2)],
                                   ncol = TotalColumnsOfIterates-2)
    LastDifference       = matrix(  (Iterates[,TotalColumnsOfIterates]-Iterates[,(TotalColumnsOfIterates-1)])     , ncol = 1)
    InverseDifferences   = MASS::ginv(Differences)
    cVector              = -InverseDifferences %*% LastDifference
    cVector              = rbind(cVector,1)
    s                    = (Iterates[,2:TotalColumnsOfIterates]%*%cVector)/sum(cVector)
    return(s)
  }
  if (Method == "RRE"){
    TotalColumnsOfIterates = dim(Iterates)[2]
    FirstColumn          = matrix(Iterates[,1], ncol = 1)
    Differences          = matrix(Iterates[,2:(TotalColumnsOfIterates)] - Iterates[,1:(TotalColumnsOfIterates-1)],
                                   ncol = TotalColumnsOfIterates-1)
    SecondDifferences    = matrix(Differences[,2:(TotalColumnsOfIterates-1)] - Differences[,1:(TotalColumnsOfIterates-2)],
                                   ncol = TotalColumnsOfIterates-2)
    FirstDifference      = matrix(Differences[,1], ncol = 1)
    Differences          = matrix(Differences[, 1:(TotalColumnsOfIterates-2)], ncol =  TotalColumnsOfIterates-2)
    InverseSecondDifferences = MASS::ginv(SecondDifferences)
    s = FirstColumn - Differences %*% InverseSecondDifferences %*% FirstDifference
    return(s)
  }
}


#' EpsilonExtrapolation
#' This function takes a matrix with previous iterates and extrapolates the limit of the sequence.
#' @export
#' @import MASS
#' @param Iterates A matrix representing different iterates with one iterate per column. Can be pieced together from Inputs and Outputs matrices
#' of the FixedPoint function using the PutTogetherIteratesWithoutJumps function.
#' @param Method Method for epsilon extrapolation. Should be either "VEA" for the vector extrapolation algorithm or "SEA" for the scalar epsilon
#' algorithm.
#' @return A vector with the extrapolated vector.
#' @examples
#' FPFunction = function(x){c(0.5*sqrt(abs(x[1] + x[2])), 1.5*x[1] + 0.5*x[2])}
#' A = FixedPoint( Function = FPFunction, Inputs = c(0.3,900), MaxIter = 6, Method = "Simple")
#' Iterates = PutTogetherIteratesWithoutJumps(A$Inputs, A$Outputs)
#' EpsilonExtrapolation(Iterates, "VEA")
#' B = FixedPoint( function(x){cos(x)}, Inputs = 1, MaxIter = 5, Method = "Simple")
#' Iterates = PutTogetherIteratesWithoutJumps(B$Inputs, B$Outputs)
#' EpsilonExtrapolation(Iterates, "SEA")
EpsilonExtrapolation = function(Iterates, Method = c("VEA", "SEA")){
    if (dim(Iterates)[2] %% 2 == 0){Iterates = matrix(Iterates[,2:dim(Iterates)[2]], ncol = dim(Iterates)[2]-1 )}
    # The function cannot do anything to a one dimensional input so will return input unchanged.
    if (dim(Iterates)[2] == 1){return(Iterates)}
    Method = Method[1]
    if (!(Method %in% c("VEA", "SEA"))){stop("Invalid method input. EpsilonExtrapolation function can only take Method as VEA or SEA")}

    Matrix = Iterates
    RowsOfMatrix    = dim(Matrix)[1]
    TotalColumnsOfMatrix = dim(Matrix)[2]
    PreviousMatrix = matrix(rep(0, RowsOfMatrix*(TotalColumnsOfMatrix-1)), nrow = RowsOfMatrix )

    for (MatrixColumns in TotalColumnsOfMatrix:(2)){
      NewMatrix = PreviousMatrix + EpsilonExtrapolationVectorOfInverses(matrix(Matrix[,2:MatrixColumns] - Matrix[,1:(MatrixColumns-1)],
                                                                                nrow = RowsOfMatrix), Method = Method)
      PreviousMatrix = Matrix[,2:(MatrixColumns-1)]
      Matrix = NewMatrix
    }

    # The function can get NAs from the inversion (ie if differenceMatrix contains a zero for SEA then there is division by zero).
    # To avert this we try with 2 less columns.
    if (any(is.na(Matrix))){Matrix = EpsilonExtrapolation( matrix(Iterates[,3:dim(Iterates)[2]], ncol = dim(Iterates)[2]-2 )  ,Method)    }
  return(matrix(Matrix, nrow = RowsOfMatrix))
}

#' EpsilonExtrapolationVectorOfInverses
#' This is a helper function for EpsilonExtrapolation
#' @param DifferenceMatrix The matrix of the differences in elements to be inverted.
#' @param Method "SEA" or "VEA".
#' @return A vector of the result of inverting each (column) vector in a mmatrix.
EpsilonExtrapolationVectorOfInverses = function(DifferenceMatrix, Method){
  if (dim(DifferenceMatrix)[1] == 1 | Method == "SEA"){return(1/DifferenceMatrix)
    } else {
  return(apply(DifferenceMatrix, 2, function(x) MASS::ginv(x)))
    }
}

#' PutTogetherIteratesWithoutJumps
#' This function takes the previous inputs and outputs and assembles a matrix with both excluding jumps.
#' @export
#' @param Inputs This is an N x A matrix of previous inputs for which corresponding outputs are available. In this case N is the
#' dimensionality of the fixed point vector that is being sought (and each column is a matrix that is input to the "Function") and A is the
#' number of previous Inputs/Outputs that are being provided to the fixed point.
#' @param Outputs This is a matrix of "Function" values for each column of the "Inputs" matrix.
#' @param AgreementThreshold A parameter for determining when a column in Inputs and a column in Outputs match. They are deemed to match if
#' the sum of the absolute values of the difference in the columns is less than AgreementThreshold.
#' @return A matrix of inputs and outputs excluding jumps.
#' @examples
#'A = FixedPoint( function(x){c(0.5*sqrt(abs(x[1] + x[2])), 1.5*x[1] + 0.5*x[2])},
#' Inputs = c(0.3,900), MaxIter = 5, Method = "Simple")
#'A = FixedPoint( function(x){c(0.5*sqrt(abs(x[1] + x[2])), 1.5*x[1] + 0.5*x[2])},
#' Inputs = A$Inputs, Outputs = A$Outputs, MaxIter = 5, Method = "Aitken")
#'A = FixedPoint( function(x){c(0.5*sqrt(abs(x[1] + x[2])), 1.5*x[1] + 0.5*x[2])},
#' Inputs = A$Inputs, Outputs = A$Outputs, MaxIter = 5, Method = "Simple")
#'CombinedSequence = PutTogetherIteratesWithoutJumps(A$Inputs, A$Outputs)
PutTogetherIteratesWithoutJumps = function(Inputs, Outputs, AgreementThreshold = 10 * .Machine$double.eps){
  # Input checking
  if (!(is.matrix(Inputs))){stop("Inputs is not a matrix.")}
  if (!(is.matrix(Outputs))){stop("Outputs is not a matrix.")}
  if (any(dim(Inputs) != dim(Outputs))){stop("Inputs and Outputs matrices are not comformable.")}
  Dimensions = dim(Inputs)
  # We get a matrix of the differences between outputs and the next inputs. If they were all simple iterates these should match.
  if (Dimensions[2] == 1){return(matrix(cbind(matrix(Inputs, ncol = 1), matrix(Outputs, ncol = 1)), ncol =2))}
  Difference = matrix(Inputs[,2:(Dimensions[2])] - Outputs[,1:(Dimensions[2]-1)], nrow = Dimensions[1])
  # We add up differences to get which iterates were simple iterates.
  Sum_Of_Differences = apply(Difference, 2, function(x) {sum(abs(x))})
  Agreements = Sum_Of_Differences < AgreementThreshold
  if (all(Agreements)){return(cbind(matrix( Inputs[,1:(Dimensions[2])], nrow=Dimensions[1]), matrix(Outputs[,Dimensions[2]],nrow=Dimensions[1])))}
  # If not all iterates were simple we just take the latter iterates and return them.
  LocationsOfBreaks = which(!Agreements)
  LastBreak = LocationsOfBreaks[length(LocationsOfBreaks)]
  return(matrix(cbind(matrix(Inputs[,(LastBreak+1):(Dimensions[2])], ncol = Dimensions[2] - LastBreak),
                      matrix(Outputs[, Dimensions[2]], ncol = 1)), ncol = Dimensions[2] - LastBreak + 1))
}