R/modelTree.R

Defines functions reduce.to.canonical.tree order.tree order.and.return.string stretch.tree stretch.block find.block.terminal.node block.of.ops.is.stretched check.node.is.in.canonical.form create.model.tree.builtin.scaled create.model.tree.builtin create.model.tree.from.list generate.next.models delete.node insert.kernel.instance add.kernel summary.ModelTree print.ModelTree as.character.ModelTree node.to.string find.parent.node find.root.node create.model.tree prod.hess sum.hess prod.grad sum.grad kern.prod kern.sum generate.id

Documented in add.kernel as.character.ModelTree create.model.tree create.model.tree.builtin create.model.tree.builtin.scaled create.model.tree.from.list delete.node insert.kernel.instance print.ModelTree summary.ModelTree

ModelTree_class_name <- "ModelTree"

generate.id <- function(existing.ids, prefix=NULL) {
  id <- NA
  counter <- 0
  while(is.na(id) | !is.na(match(id, existing.ids))) {
    counter <- counter + 1
    id <- paste(prefix, counter, sep=".")
  }
  return(id)
}

kern.sum <- function(a, b) {
  return(a + b)
}

kern.prod <- function(a, b) {
  return(a * b)
}

sum.grad <- function(a.grad, b.grad, ...) {
  return(a.grad + b.grad)
}

prod.grad <- function(a, b, a.grad, b.grad) {
  a1 <- b1 <- b.grad
  for (i in 1:dim(b.grad)[3]) {
    a1[, , i] <- a
    b1[, , i] <- b
  }
  return(a1 * b.grad + b1 * a.grad)
}

sum.hess <- function(a.hess, b.hess, ...) {
  return(a.hess + b.hess)
}

prod.hess <- function(a, b, a.grad, b.grad, a.hess, b.hess) {
  a1 <- b1 <- b.hess
  for (i in 1:dim(b.grad)[3]) {
    for (j in 1:dim(b.grad)[3]) {
        a1[, , i, j] <- a
        b1[, , i, j] <- b
    }
  }

  a.grad1 <- a.grad2 <- b.grad1 <- b.grad2 <- b.hess
  for (i in 1:dim(b.grad)[3]) {
    a.grad1[, , i, ] <- a.grad2[, , , i] <- a.grad
    b.grad1[, , i, ] <- b.grad2[, , , i] <- b.grad
  }

  return(a.hess * b1 + a.grad1 * b.grad2 + a.grad2 * b.grad1 + a1 * b.hess)
}

#' Create an Empty Model Tree
#'
#' Creates an instance of a modelTree object with no associated kernels.
#'
#' @return a modelTree object
#' @export
#'
#' @examples
#' mt <- create.model.tree()
create.model.tree <- function() {
  model.tree <- list()
  class(model.tree) <- ModelTree_class_name

  model.tree$operation.functions = list("+"=kern.sum,
                                        "*"=kern.prod
  )

  model.tree$operation.grad.functions = list("+"=sum.grad,
                                             "*"=prod.grad
  )

  model.tree$operation.hess.functions = list("+"=sum.hess,
                                             "*"=prod.hess
  )

  model.tree$tree <- data.frame(leftDaughter=numeric(0),
                                rightDaughter=numeric(0),
                                kernelInstanceID=character(0),
                                operationID=character(0),
                                stringsAsFactors=FALSE
  )

  model.tree$kernel.instances <- data.frame(kernelInstanceID=character(0),
                                            kernelClassName=character(0),
                                            stringsAsFactors=FALSE
  )

  model.tree$all.hyper.params <- c()

  # These are accessed by string IDs to make sure there's no issue
  # with accidentally accessing the nth element when we want the element
  # whose ID is "n".

  # This gives the indices for the hyperparameters in the hyperparameters
  # vector for a specific instance of a kernel. Indexed by kernelInstanceID.
  model.tree$kernel.inst.hyper.param.indices <- list()

  # This holds the base kernels (of class Kernel) available to the tree.
  # Indexed by kernelClassName.
  model.tree$kernel.objects <- list()

  return(model.tree)
}

find.root.node <- function(model.tree) {
  num.nodes <- nrow(model.tree$tree)
  if (num.nodes == 0) {
    root.node <- 1
  } else {
    child.nodes <- na.omit(c(model.tree$tree$leftDaughter, model.tree$tree$rightDaughter))
    root.node <- which(!1:num.nodes %in% child.nodes)
  }
  if (length(root.node) == 0){
    stop("No root node found.")
  }
  if (length(root.node) > 1){
    stop("Multiple root nodes found.")
  }
  return(root.node)
}

find.parent.node <- function(model.tree, node) {
  parent.node <- c(which(node == model.tree$tree$leftDaughter),
                   which(node == model.tree$tree$rightDaughter)
  )
  if (length(parent.node) > 1) {
    stop("Multiple parents found.")
  }
  if (length(parent.node) == 0) {
    # Presumably this is the root node. Return NULL.
    parent.node <- NULL
  }
  return(parent.node)
}

node.to.string <- function(model.tree, node, show.node.ids=FALSE) {
  if (nrow(model.tree$tree) == 0) {
    ret <- NULL
  } else {
    op <- model.tree$tree[node, "operationID"]
    ld <- model.tree$tree[node, "leftDaughter"]
    rd <- model.tree$tree[node, "rightDaughter"]
    if (is.na(op)) {
      ret <- model.tree$tree[node, "kernelInstanceID"]
    } else {
      ret <- paste("(",
                   node.to.string(model.tree, ld, show.node.ids),
                   op,
                   node.to.string(model.tree, rd, show.node.ids),
                   ")",
                   sep=" "
      )
    }
  }
  if (show.node.ids) {
    ret <- paste("[", node, ": ", ret, " :", node, "]",
                 sep=""
                 )
  }
  return(ret)
}

#' Convert a Model Tree to a String
#'
#' Return or print a string representation of the combination of kernels represented by a given model tree
#'
#' @param model.tree a modelTree object
#' @param show.node.ids boolean indicating whether the string should be annotated with the associated node numbers.
#'
#' @return a string representing a modelTree object
#' @export
#'
#' @examples
#' mt <- create.model.tree.builtin()
#' mt <- insert.kernel.instance(mt, 1, "SE", NULL)
#' mt <- insert.kernel.instance(mt, 1, "SE", "+")
#'
#' as.character(mt) #"( SE.1 + SE.2 )"
#'
#' as.character(mt, show.node.ids=TRUE) #"[3: ( [1: SE.1 :1] + [2: SE.2 :2] ) :3]"
#'
#' print(mt) #"( SE.1 + SE.2 )"
#'
#' print(mt, show.node.ids=TRUE) #"[3: ( [1: SE.1 :1] + [2: SE.2 :2] ) :3]"
as.character.ModelTree <- function(model.tree, show.node.ids=FALSE) {
  root.node <- find.root.node(model.tree)
  return(node.to.string(model.tree, root.node, show.node.ids))
}

#' @rdname as.character.ModelTree
#' @export
print.ModelTree <- function(model.tree, show.node.ids=FALSE) {
  cat(as.character.ModelTree(model.tree, show.node.ids))
  cat("\n")
}

#' Print a Summary of a Model Tree
#'
#' @param model.tree ModelTree object to summarise
#'
#' @export
#'
summary.ModelTree <- function(model.tree) {
  cat("Operation functions:\n")
  for (op in names(model.tree$operation.functions)) {
    cat(paste('- "', op, '"\n', sep=""))
  }

  cat("\nKernel Classes:\n")
  for (kernel in names(model.tree$kernel.objects)) {
    cat(paste('- ', kernel, '\n', sep=""))
    print(model.tree$kernel.objects[[kernel]])
  }

  cat("\nKernel Instances:\n")
  for (i in seq_along(model.tree$kernel.instances[, 1])) {
    kinst_ID <- model.tree$kernel.instances$kernelInstanceID[i]
    kinst_class <- model.tree$kernel.instances$kernelClassName[i]
    cat(paste('- ', kinst_ID, '\n', sep=""))
    cat(paste('  - Class:', kinst_class, '\n', sep=""))
    cat(paste('  - Hyperparameters:\n', sep=""))
    hp_vec <- model.tree$kernel.inst.hyper.param.indices[[kinst_ID]]
    for (hp in seq_along(hp_vec)) {
      if (is.null(model.tree$all.hyper.params[hp_vec[hp]])) {
        hp_val <- "NULL"
      } else {
        hp_val <- model.tree$all.hyper.params[hp_vec[hp]]
      }
      cat(paste('    - ', model.tree$kernel.objects[[kinst_class]][["hyperparam_names"]][hp],
                ": ", hp_val, '\n', sep=""))
    }
  }

  cat("\nModel Tree Formula:\n")
  cat(paste("  ", as.character(model.tree), "\n", sep=""))

  cat("\nModel Tree Formula (with node number annotations):\n")
  cat(paste("  ", as.character(model.tree, show.node.ids=T), "\n", sep=""))
}

#' Add a Kernel Class to a Model Tree
#'
#' @param orig.model.tree a ModelTree object
#' @param kernel.class.name the name of the kernel class (arbitrary, must not
#' match the class name of any kernel already in the model tree).
#' @param kernel either an R function defining a kernel, or a string indicating a built-in kernel type
#' @param hyper.param.names a vector containing the names of the kernel hyperparameters
#' @param kernel.additional.params a list containing any additional kernel parameters
#' @param kernel.grad either an R function returning the gradient of \code{kernel}, or NULL. If \code{kernel}
#' is a string specifying a built-in kernel, \code{kernel.grad} should be NULL. If \code{kernel} is an
#' R function and \code{kernel.grad} is NULL, \code{\link[numDeriv]{grad}} will be used to derive a numerical
#' gradient.
#'
#' @return a ModelTree object
#' @export
#'
#' @examples
#' mt <- create.model.tree()
#' mt <- add.kernel(mt, "SE", "squaredExponential", c("l"))
add.kernel <- function(orig.model.tree,
                       kernel.class.name,
                       kernel) {
  model.tree <- orig.model.tree

  if (kernel.class.name %in% names(model.tree$kernel.objects)) {
    stop("Kernel with that name is already included")
  }

  if (grepl(pattern="#", x=kernel.class.name)) {
    stop("Kernel class name cannot contain '#'")
  }

  model.tree$kernel.objects[[kernel.class.name]] <- kernel


  return(model.tree)
}

#' Insert Kernel Instance into Model Tree
#'
#' @param orig.model.tree a ModelTree object
#' @param node the node at which to insert the kernel instance
#' @param kernel.class.name a string indicating the class of kernel to insert
#' @param operation.id a string indicating how to insert the kernel at \code{node}.
#' Must be "*" or "+", or NULL if the model tree is currently empty.
#' @param hyper.params the hyper parameters of the inserted kernel. Can be NULL to
#' fit hyper parameters at a later time.
#'
#' @return a ModelTree object
#' @export
#'
#' @examples
#' mt <- create.model.tree.builtin()
#'
#' mt <- insert.kernel.instance(mt, 1, "SE", NULL)
#' print(mt) # "SE.1"
#'
#' mt <- insert.kernel.instance(mt, 1, "SE", "*")
#' print(mt) # "( SE.1 * SE.2 )"
#'
#' mt <- insert.kernel.instance(mt, 3, "RQ", "+")
#' print(mt) # "( ( SE.1 * SE.2 ) + RQ.1 )"
#'
insert.kernel.instance <- function(orig.model.tree,
                                   node,
                                   kernel.class.name,
                                   operation.id,
                                   hyper.params=NULL) {
  model.tree <- orig.model.tree

  if (!kernel.class.name %in% names(model.tree$kernel.objects)) {
    stop("Unknown kernel specified.")
  }
  if (node > max(nrow(model.tree$tree), 1) | node <= 0 | floor(node) != node) {
    stop("Invalid node specified.")
  }

  # Check there's no operation if it's the first node
  if(nrow(model.tree$tree)==0 & node==1) {
    if (!is.null(operation.id)) {
      stop("Tree is empty - no operation should be specified for first node inserted.")
    }
  } else {
    # If it's not the first node we need a valid operation.
    if (is.null(operation.id)) {
      stop("No operation specified.")
    } else if (!operation.id %in% names(model.tree$operation.functions)) {
      stop("Invalid operation specified")
    }
  }

  # Check the hyperparameters look right
  hyperparam_names <- model.tree$kernel.objects[[kernel.class.name]]$hyperparam_names
  if (length(hyperparam_names) > 0) {
    if (!is.null(hyper.params)) {
      # If we expect some and have been given some, do the two match?
      if (!all(names(hyper.params) %in% hyperparam_names)) {
        stop("Supplied hyperparameters do not match hyperparameters expected for this kernel.")
      }
    }
  } else {
    # We're not expecting any hyperparams.
    if (!is.null(hyper.params)) {
      stop("Hyperparameters supplied when none are expected for this kernel.")
    }
  }


  inst.id <- generate.id(model.tree$kernel.instances[,1], kernel.class.name)

  # Add to the list of instances
  model.tree$kernel.instances[inst.id, ] <- c(kernelInstanceID=inst.id,
                                              kernelClassName=kernel.class.name
  )


  # Append or create hyper parameters. NAs indicate random start.
  # Are we expecting any hyperparameters?
  if (length(hyperparam_names) > 0) {
    # Have any been supplied?
    if (!is.null(hyper.params)) {
      # Re-order if necessary to match our stored list.
      hyper.params <- hyper.params[hyperparam_names]
    } else {
      # Make NA hyperparams to indicate random start
      hyper.params <- rep(NA, length(hyperparam_names))
    }
    names(hyper.params) <- paste(inst.id, hyperparam_names, sep="#")

    # Add them to the list of hyper.params
    start.index <- length(model.tree$all.hyper.params) + 1
    model.tree$all.hyper.params <- c(model.tree$all.hyper.params, hyper.params)
    end.index <- length(model.tree$all.hyper.params)

    model.tree$kernel.inst.hyper.param.indices[[inst.id]] <- start.index:end.index

  } else {
    # We're not expecting any hyperparams. Create an empty hyperparam indices list.
    model.tree$kernel.inst.hyper.param.indices[[inst.id]] <- c()
  }

  kernel.inst.node.entry <- c(leftDaughter=NA,
                              rightDaughter=NA,
                              kernelInstanceID=inst.id,
                              operationID=NA
  )

  # Find the root node of this model tree. environment() returns the containing model tree env.
  root.node <- find.root.node(model.tree)

  # Are we inserting the first node?
  if (nrow(model.tree$tree)==0 & node==1) {
    model.tree$tree[1, ] <- kernel.inst.node.entry
  } else {
    new.inst.node <- nrow(model.tree$tree) + 1
    new.op.node <- nrow(model.tree$tree) + 2


    new.op.node.entry <- c(leftDaughter=node,
                           rightDaughter=new.inst.node,
                           kernelInstanceID=NA,
                           operationID=operation.id
    )

    # If node is the root node we're done. Otherwise we need to fix the parent
    # to point at our new operation node.
    if (node!=root.node) {
      # Find the parent to this node so we can have it point at our new
      # operation node rather than the current node.
      parent.node <- find.parent.node(model.tree, node)

      node.is.left.daughter <- model.tree$tree[parent.node, "leftDaughter"] == node

      # point the parent at our new op node
      if (node.is.left.daughter) {
        model.tree$tree[parent.node, "leftDaughter"] <- new.op.node
      } else {
        # If it's not left then it must be right
        model.tree$tree[parent.node, "rightDaughter"] <- new.op.node
      }
    }

    # Insert our new nodes
    model.tree$tree[new.inst.node, ] <- kernel.inst.node.entry
    model.tree$tree[new.op.node, ] <- new.op.node.entry

  }
  return(model.tree)
}

#' Delete a Node in a Model Tree
#'
#' @param model.tree a ModelTree object
#' @param node the node to delete (must not be the root node)
#'
#' @return a copy of \code{model.tree} with \code{node} deleted.
#' @export
delete.node <- function(model.tree, node) {
  root.node <- find.root.node(model.tree)
  if (node == root.node) {
    stop("Can't delete the root node.")
  }

  new.mt <- model.tree

  # Repoint the grandparent node at this node's sibling
  parent.node <- find.parent.node(new.mt, node)

  if (new.mt$tree[parent.node, "leftDaughter"] == node) {
    sibling.node <- new.mt$tree[parent.node, "rightDaughter"]
    new.mt$tree[parent.node, "rightDaughter"] <- NA
  } else {
    sibling.node <- new.mt$tree[parent.node, "leftDaughter"]
    new.mt$tree[parent.node, "leftDaughter"] <- NA
  }

  if (!parent.node == root.node) {
    grandparent.node <- find.parent.node(new.mt, parent.node)

    if (new.mt$tree[grandparent.node, "leftDaughter"] == parent.node) {
      new.mt$tree[grandparent.node, "leftDaughter"] <- sibling.node
    } else {
      new.mt$tree[grandparent.node, "rightDaughter"] <- sibling.node
    }
  }

  # Now delete the orphaned nodes (by overwriting with NAs then deleting the rows)
  deleted.nodes <- numeric(0)
  node.queue <- c(parent.node)
  while (length(node.queue) > 0) {
    current.node <- node.queue[1]
    node.queue <- node.queue[-1]
    node.row <- new.mt$tree[current.node, ]
    if (!is.na(node.row$operationID)) {
      node.queue <- as.numeric(c(node.queue,
                      new.mt$tree[current.node, "leftDaughter"],
                      new.mt$tree[current.node, "rightDaughter"]))
      node.queue <- node.queue[!is.na(node.queue)]
    }
    new.mt$tree[current.node, ] <- NA

    deleted.nodes <- c(deleted.nodes, current.node)
  }

  node.indices <- seq(nrow(new.mt$tree))[-deleted.nodes]

  new.mt$tree <- new.mt$tree[-deleted.nodes, ]

  rownames(new.mt$tree) <- 1:nrow(new.mt$tree)

  # We need to update the node references in the tree
  for (i in 1:nrow(new.mt$tree)) {
    for (child in c("leftDaughter", "rightDaughter")) {
      if (!is.na(new.mt$tree[i, child])) {
        new.mt$tree[i, child] <- which(node.indices == new.mt$tree[i, child])
      }
    }
  }

  remaining.inst.ids <- unique(new.mt$tree[, "kernelInstanceID"])
  deleted.hyper.param.indices <- c()

  for (i in seq_along(new.mt$kernel.instances[, "kernelInstanceID"])) {
    inst.id <- new.mt$kernel.instances[i, "kernelInstanceID"]
    if (!(inst.id %in% remaining.inst.ids)) {
      new.mt$kernel.instances <- new.mt$kernel.instances[-i, ]

      deleted.hyper.param.indices <- c(deleted.hyper.param.indices,
                                       new.mt$kernel.inst.hyper.param.indices[[inst.id]])

      temp.index <- which(names(new.mt$kernel.inst.hyper.param.indices) == inst.id)
      new.mt$kernel.inst.hyper.param.indices <- new.mt$kernel.inst.hyper.param.indices[-temp.index]
    }
  }

  # If we deleted some hyper parameters we need to translate the old indices to
  # the new indices.
  if (length(deleted.hyper.param.indices) > 0) {
    hyper.param.indices <- seq_along(new.mt$all.hyper.params)[-deleted.hyper.param.indices]
    new.mt$all.hyper.params <- new.mt$all.hyper.params[-deleted.hyper.param.indices]
    for (i in seq_along(new.mt$kernel.inst.hyper.param.indices)) {
      for (j  in seq_along(new.mt$kernel.inst.hyper.param.indices[[i]])) {
        new.mt$kernel.inst.hyper.param.indices[[i]][j] <-
          which(hyper.param.indices == new.mt$kernel.inst.hyper.param.indices[[i]][j])
      }
    }
  }

  return(new.mt)
}

generate.next.models <- function(model.tree) {

  # Generate a list of models
  models <- list()
  i <- 1
  if (nrow(model.tree$tree) == 0) {
    for (kernel in names(model.tree$kernel.objects)) {
      models[[i]] <- insert.kernel.instance(orig.model.tree=model.tree,
                                            node=1,
                                            kernel.class.name=kernel,
                                            operation.id=NULL)
      i <- i + 1
    }
  } else {
    for (operation in names(model.tree$operation.functions)) {
      for (kernel in names(model.tree$kernel.objects)) {
        for (node in 1:nrow(model.tree$tree)) {
          models[[i]] <- insert.kernel.instance(orig.model.tree=model.tree,
                                                node=node,
                                                kernel.class.name=kernel,
                                                operation.id=operation)
          i <- i + 1
        }
      }
    }
    root.node <- find.root.node(model.tree)
    for (node in 1:nrow(model.tree$tree)) {
      if (node != root.node) {
        models[[i]] <- delete.node(model.tree, node)
        i <- i + 1
      }
    }
  }

  # Reduce to canonical form and dedupe

  model.strings <- character(length(models))
  for (i in seq_along(models)) {
    models[[i]] <- reduce.to.canonical.tree(models[[i]])
    model.strings[i] <- as.character(models[[i]])
  }

  duplicates <- which(duplicated(model.strings))

  if (length(duplicates) > 0) {
    models <- models[-duplicates]
  }

  return(models)
}

#' Create a Model Tree from a list of kernel objects
#'
#' Returns a model tree instance populated with a basis set of kernels.
#'
#' @param kernel_list a named list of kernel objects
#'
#' @return A modelTree object
#' @export
create.model.tree.from.list <- function(kernel_list) {

  mt <- create.model.tree()
  for (kernel_name in names(kernel_list)) {
    kernel <- kernel_list[[kernel_name]]
    mt <- add.kernel(mt,
                     kernel_name,
                     kernel
    )
  }

  return(mt)
}

#' Create a Model Tree Containing All Built-in Kernels
#'
#' Returns a model tree instance populated with a basis set of all of the built in kernels.
#'
#' @return A modelTree object
#' @export
#'
#' @examples
#' mt <- create.model.tree.builtin()
create.model.tree.builtin <- function() {
  builtin.kernels <- list_built_in_kernels()

  mt <- create.model.tree.from.list(builtin.kernels)

  return(mt)
}

#' Create a Model Tree Containing All Built-in Kernels (scaled)
#'
#' Returns a model tree instance populated with a basis set of all of the built in kernels.
#'
#' @return A modelTree object
#' @export
#'
#' @examples
#' mt <- create.model.tree.builtin.scaled()
create.model.tree.builtin.scaled <- function() {
  builtin.kernels <- list_built_in_kernels_scaled()

  mt <- create.model.tree.from.list(builtin.kernels)

  return(mt)
}


check.node.is.in.canonical.form <- function(model.tree, node) {
  is.canonical <- TRUE

  node.row <- model.tree$tree[node, ]

  # If this is a terminal (non-operation) node then it is canonical. If not...
  if (!is.na(node.row$operationID)) {
    # check left daughter operation does not match this node's
    if (!identical(model.tree$tree[node.row$leftDaughter, "operationID"], node.row$operationID)) {
      is.canonical <- FALSE
    } else {
      # check the left daughter of this node is smaller than the left daughter of it's right daughter,
      # if they are the same type of operation
      right.daughter.node.row <- model.tree$tree[node.row$rightDaughter, ]
      if (identical(right.daughter.node.row$operationID, node.row$operationID)) {
        right.daughters.left.daughter.string <- node.to.string(right.daughter.node.row$leftDaughter)
        left.daughter.string <- node.to.string(node.row$leftDaughter)
        if (right.daughters.left.daughter.string < left.daughter.string) {
          is.canonical <- FALSE
        }
      }
    }
  }
  return(is.canonical)
}

block.of.ops.is.stretched <- function(model.tree, node, operationID) {
  node.row <- model.tree$tree[node, ]

  if (is.na(operationID)) {
    return(TRUE)
  }

  if (identical(operationID, node.row$operationID)){
    left.daughter.op <- model.tree$tree[node.row$leftDaughter, "operationID"]
    if (identical(operationID, left.daughter.op)) {
      return(FALSE)
    } else {
      return(block.of.ops.is.stretched(model.tree, node.row$rightDaughter, operationID))
    }
  } else {
    return(TRUE)
  }
}

find.block.terminal.node <- function(model.tree, node, operationID) {
  tree <- model.tree$tree
  node.row <- tree[node,]
  if (!identical(node.row$operationID, operationID)) {
    stop("block root node is not of the specified operation")
  }
  current.node <- node
  next.node <- tree[current.node, "rightDaughter"]
  while (identical(tree[next.node, "operationID"], operationID)) {
    current.node <- next.node
    next.node <- tree[next.node, "rightDaughter"]
  }
  return(current.node)
}

stretch.block <- function(model.tree, node, operationID) {
  node.row <- model.tree$tree[node,]
  if (!identical(node.row$operationID, operationID)) {
    stop("block root node is not of the specified operation")
  }
  current.node <- node
  terminal.node <- find.block.terminal.node(model.tree, current.node, operationID)
  while (identical(model.tree$tree[current.node, "operationID"], operationID)) {
    left.daughter <- model.tree$tree[current.node, "leftDaughter"]
    if (identical(model.tree$tree[left.daughter, "operationID"], operationID)) {
      # We need a non-matching op as left daughter, so we'll swap the left
      # daughter with the right daughter of the terminal node
      terminal.node.right.daughter <- model.tree$tree[terminal.node, "rightDaughter"]
      model.tree$tree[current.node, "leftDaughter"] <- terminal.node.right.daughter
      model.tree$tree[terminal.node, "rightDaughter"] <- left.daughter

      # Now we'll have a new terminal node, so let's go find it
      terminal.node <- find.block.terminal.node(model.tree, current.node, operationID)
    }

    # We're done with this node, so let's move on.
    current.node <- model.tree$tree[current.node, "rightDaughter"]
  }
  return(model.tree)
}

stretch.tree <- function(model.tree) {
  root.node <- find.root.node(model.tree)
  node.queue <- c(root.node)
  while (length(node.queue) > 0) {
    #Pop off of the front of our wildly inefficient "queue"
    current.node <- node.queue[1]
    node.queue <- node.queue[-1]

    current.op <- model.tree$tree[current.node, "operationID"]

    if (!block.of.ops.is.stretched(model.tree, current.node, current.op)) {
      model.tree <- stretch.block(model.tree, current.node, current.op)
    }

    if (!is.na(current.op)) {
      node.queue <- c(node.queue, as.numeric(model.tree$tree[current.node, c("leftDaughter", "rightDaughter")]))
    }
  }
  return(model.tree)
}

order.and.return.string <- function(model.tree, node) {
  # We assume the tree has been stretched i.e. commutable operations proceed
  # in blocks down the right daughters
  node.op <- model.tree$tree[node, "operationID"]
  if (!is.na(node.op)) {
    left.daughter <- model.tree$tree[node, "leftDaughter"]
    right.daughter <- model.tree$tree[node, "rightDaughter"]
    right.daughter.op <- model.tree$tree[right.daughter, "operationID"]
    if (!identical(right.daughter.op, node.op)) {
      # Terminal string of a block - get strings and re-order if necessary
      left.temp <- order.and.return.string(model.tree, left.daughter)
      left.daughter.string <- left.temp$string
      model.tree <- left.temp$model.tree

      right.temp <- order.and.return.string(model.tree, right.daughter)
      right.daughter.string <- right.temp$string
      model.tree <- right.temp$model.tree
      if (left.daughter.string > right.daughter.string) {
        model.tree$tree[node, "leftDaughter"] <- right.daughter
        model.tree$tree[node, "rightDaughter"] <- left.daughter
      }
    } else {
      # We're in a block of ops - need to look at the ops around us.
      # If we're in a block then "node" must be the root node of the
      # block, since we never call this function on right daughters in
      # the block.

      # Get the block's left daughter strings
      current.block.node <- node
      block.nodes <- c()
      block.left.daughters <- c()
      block.left.daughter.strings <- c()
      while (identical(model.tree$tree[current.block.node, "operationID"], node.op)) {
        block.nodes <- c(block.nodes, current.block.node)
        block.left.daughter <- model.tree$tree[current.block.node, "leftDaughter"]
        block.left.daughters <- c(block.left.daughters, block.left.daughter)
        block.left.daughter.temp <- order.and.return.string(model.tree,
                                                            block.left.daughter)
        model.tree <- block.left.daughter.temp$model.tree
        block.left.daughter.strings <- c(block.left.daughter.strings,
                                         block.left.daughter.temp$string)
        current.block.node <- model.tree$tree[current.block.node, "rightDaughter"]
      }

      # current.block.node now contains the right daughter of the terminal node -
      # we need to include this in the ordering
      terminal.right.daughter <- current.block.node
      block.left.daughters <- c(block.left.daughters, terminal.right.daughter)
      block.left.daughter.temp <- order.and.return.string(model.tree,
                                                          terminal.right.daughter)
      model.tree <- block.left.daughter.temp$model.tree
      block.left.daughter.strings <- c(block.left.daughter.strings,
                                       block.left.daughter.temp$string)

      correct.daughter.order <- order(block.left.daughter.strings)
      ordered.daughters <- block.left.daughters[correct.daughter.order]
      for (block.node.index in seq_along(block.nodes)) {
        block.node <- block.nodes[block.node.index]
        new.left.daughter <- ordered.daughters[block.node.index]
        model.tree$tree[block.node, "leftDaughter"] <- new.left.daughter
      }

      # Now we update the right daughter of the terminal node
      block.node <- tail(block.nodes, 1)
      new.right.daughter <- tail(ordered.daughters, 1)
      model.tree$tree[block.node, "rightDaughter"] <- new.right.daughter
    }
  }
  return(list(string=node.to.string(model.tree, node, show.node.ids=FALSE),
              model.tree=model.tree))
}

order.tree <- function(model.tree) {
  root.node <- find.root.node(model.tree)
  model.temp <- order.and.return.string(model.tree, root.node)
  model.tree <- model.temp$model.tree
  return(model.tree)
}

reduce.to.canonical.tree <- function(model.tree) {
  # First we stretch the tree, then we order it.
  model.tree <- stretch.tree(model.tree)
  model.tree <- order.tree(model.tree)
  return(model.tree)
}
mattdneal/gaussianProcess documentation built on May 21, 2019, 12:58 p.m.