directedENA coverage - 26.83%

Files
Source

#' Generate directed ENA network
#'
#' @param data data.frame
#' @param ... additional parameters. `accumulated_data` to bypass accumulation, also passed to plot
#' @param optimize_node_positions logical
#' @param plotted_points data.frame
#' @param node_positions data.frame
#' @param print_plot logical
#' @param only_sending_centroid logical
#' @param units TBD
#' @param codes TBD
#' @param conversations TBD
#' @param rotation_on TBD
#' @param optimize_on TBD
#'
#' @return object of type ena.set.directed
#' @export
#'
#' @examples
#'
#' data("RS.data", package = "rENA")
#' rs_code_cols <- colnames(RS.data)[15:20]
#' rs_set_one <- directed_ena(RS.data,
#'    units = c("Condition", "GroupName", "UserName"),
#'    conversations = c("Condition", "GroupName", "ActivityNumber"),
#'    codes = rs_code_cols,
#'    binary = TRUE,
#'    rotation_on = "response",
#'    optimize_on = "response",
#'    windowSize = 4, print_plot = FALSE
#' )
directed_ena = function(
  data, units, conversations, codes,
  ...,
  optimize_node_positions = TRUE,
  only_sending_centroid = FALSE,
  plotted_points = NULL,
  node_positions = NULL,
  rotation_on = "response",
  optimize_on = c("ground", "response"),
  print_plot = FALSE
) {
  args <- list(...);

  set <- ena.set.directed(data, units, conversations, codes);

  if(!is.null(args$accumulated_data)) {
    data <- args$accumulated_data;
    if(is.data.frame(data$weights[[1]])) {
      adjacency_matrices <- data$weights;
    }
    else {
      adjacency_matrices = format_adjacency_matrix(data$weights);
    }
    set$adjacency_matrix = generate_adjacency_matrix(adjacency_matrices);
  }
  else {
    # This should create set$connection.counts and set$model$row.connection.counts
    set <- directed_accumulation(x = set, units = units, conversations = conversations, codes = codes, ...);
  }

  if(is.null(plotted_points)) {
    set <- directed_model(set, rotation_on = rotation_on, optimize_on = optimize_on, ...);
  }
  else {
    set$points <- plotted_points;
  }

  if(is.null(node_positions)) {
    if(optimize_node_positions == FALSE) {
      set$rotation$nodes = generate_node_coordinates(codes, total_node_connections_calculator(abs((as.square.data.frame(colMeans(as.matrix(set$connection.counts)), names = set$`_function.params`$codes)))))[];
    }
    else {
      # pass abs(adjacency_matrix) in order to account for subtraction matrix
      # set$nodes = generate_node_coordinates(set$codes, NULL);

      # adjacency_matrices <- list(as.square.data.frame(as.numeric(set$connection.counts[1,2:10]), names = set$codes));
      # if (only_sending_centroid) {
      #   set$centroid_coordinates = calculate_centroid(adjacency_matrices);
      # }
      # else {
      #   # NOTE - directed network centroid is midpoint between tail and head of centroid vector
      #   set$centroid_coordinates = calculate_centroid_vectors(adjacency_matrices);
      # }
      # set$rotation$nodes$node_radius <- max_node_diameter(set$rotation$nodes) / 2;


      # if(!is.null(connections_calculator)) {
      #   set$rotation$nodes[, size := sapply(code, connections_calculator)];
      # }
      # else {
        # set$rotation$nodes[, size := 1];
      # }

      # max_size = max(set$rotation$nodes$size);

      # dims <- set$rotation$nodes[, rENA::find_dimension_cols(set$rotation$nodes), with = F];
      # set$rotation$nodes$x0 <- dims[, 1] - node_radius
      # set$rotation$nodes$x1 <- dims[, 1] + node_radius
      # set$rotation$nodes$y0 <- dims[, 2] - node_radius
      # set$rotation$nodes$y1 <- dims[, 2] + node_radius
    }
  }
  else {
    set$nodes <- node_positions;
  }

  if(print_plot == TRUE) {
    set <- plot.ena.directed.set(x = set, ...);
    invisible(set);
  }
  else {
    return(set);
  }
}

#' Create directed ena set object
#'
#' @param data data.frame
#' @param ... Not implemented
#' @param units TBD
#' @param conversations TBD
#' @param codes TBD
#'
#' @return empty list of type ena.set.directed
#' @export
ena.set.directed <- function(data, units, conversations, codes, ...) {
  as.ena.directed.set(
    list(
      meta.data = NULL,
      model = list(
        model.type = "Directed",
        raw.input = data.table::as.data.table(data)
      ),
      rotation = list (
        codes = codes,
        nodes = NULL
      ),
      plots = list(),
      "_function.params" = list(
        units = units,
        conversations = conversations,
        codes = codes
      )
    )
  )
}



#' Plot a directed ena model
#'
#' @param x ena.set.directed to plot
#' @param ... parameters to pass along - Not yet implemented
#' @param only_sending_centroid logical
#' @param self_connections_in_center logical
#' @param node_size_dynamic logical
#' @param node_color Default: "#000000"
#' @param print_plot logical
#' @param multiplier logical
#' @param multiplier_nodes logical
#' @param multiplier_edges logical
#' @param units List
#' @param dimensions numeric vector of length two
#' @param with_mean logical
#' @param with_points logical
#' @param with_edges logical
#' @param with_nodes logical
#' @param with_node_labels logical
#' @param node_center_diameter character, either "self" (default), "response", or "ground"
#' @param node_outer_diameter character, either "self", "response" (default), or "ground"
#' @param node_shape character, see: https://plotly.com/r/reference/#scatter-marker-symbol
#' @param scale_points logcical, default: TRUE
#' @param colors matrix, see directedENA:::COLORS_HSV
#' @param title character, default is empty
#' @param plot_width numeric, default is NULL allowing for dynamic sizing
#' @param plot_height see `plot_width`
#' @param edge_scale  TBD
#'
#' @return ENA model with plot attached
#' @export
#'
#' @examples
#'
#' data("RS.data", package = "rENA")
#' rs_code_cols <- colnames(RS.data)[15:20]
#' rs_set_one <- directed_ena(RS.data,
#'    units = c("Condition", "GroupName", "UserName"),
#'    conversations = c("Condition", "GroupName", "ActivityNumber"),
#'    codes = rs_code_cols,
#'    binary = TRUE,
#'    rotation_on = "response",
#'    optimize_on = "response",
#'    windowSize = 4, print_plot = FALSE
#' )
#' plot(rs_set_one)
plot.ena.directed.set <- function(
  # Model options
    x, ...,
    units = NULL,
    dimensions = 1:2,
    print_plot = TRUE,

  # Multipliers for scaling elements
    multiplier = 5,
    multiplier_nodes = multiplier,
    multiplier_edges = multiplier,

  # What to include
    with_mean = FALSE,
    with_points = FALSE,
    with_edges = TRUE,
    with_nodes = TRUE,
    with_node_labels = TRUE,

  # Configuration
    only_sending_centroid = FALSE,  # Deprecate if we remove support for custom centroid_coordinates on `x`
    self_connections_in_center = TRUE,
    node_center_diameter = "self",
    node_outer_diameter = "response",
    node_size_dynamic = TRUE,
    node_shape = "circle",
    node_color = "#000000",

    edge_scale = c(0.5,1),
    scale_points = TRUE,
    colors = NULL,
    title = NULL,
    plot_width = NULL,
    plot_height = plot_width
) {
  args <- list(...);

  if(is.null(args$print_plot)) {
    args$print_plot <- print_plot;
  }

  if(!is.null(units)) {
    if(inherits(units, "ena.points")) {
      adjacency_matrix <- colMeans(as.matrix(units))
      if(is.null(colors)) {
        colors <- COLORS_HSV[, 1, drop = FALSE]
      }
    }
    else if (inherits(units, "ena.line.weights")) {
      adjacency_matrix <- as.square.data.frame(
        matrix(colMeans(as.matrix(units)),nrow = 1),
        names = x$`_function.params`$codes
      )
      if(is.null(colors)) {
        colors <- COLORS_HSV[, 1, drop = FALSE]
      }
      units = list("Units" = x$line.weights$ENA_UNIT %in% units$ENA_UNIT)
    }
    else {
      adjacency_matrix <- Reduce(`-`, lapply(units, function(u) {
        as.square.data.frame(
          colMeans(
            as.matrix(x$line.weights[u & ENA_DIRECTION == "response"])
          ),
          names = x$`_function.params`$codes
        )
      }))
      if(is.null(colors)) {
        colors <- COLORS_HSV[, seq(length(units)), drop = FALSE]
      }
    }
  }
  else {
    adjacency_matrix <- as.square.data.frame(colMeans(as.matrix(x$line.weights[ENA_DIRECTION == "response",])), names = x$`_function.params`$codes);
    if(is.null(colors)) {
      colors <- COLORS_HSV
    }
    units <- list("Units" = rep(TRUE, nrow(x$line.weights)))
  }
  if(is.character(colors)) {
    colors <- apply(col2rgb(colors), 2, rgb2hsv)
  }

  nodes <- data.table::copy(x$rotation$nodes);
  nodes$size = 1
  node_radius = max_node_diameter(nodes) / 2;
  nodes[, ':=' (relative_size = size / max(size)) ];
  nodes[, ':=' (node_radius = relative_size * node_radius) ];
  all_dim_cols <- find_dimension_cols(nodes);
  dim_cols <- which(all_dim_cols)[dimensions];
  other_cols <- which(!all_dim_cols);
  nodes <- nodes[, .SD, .SDcols = c(names(other_cols), names(dim_cols)), with = TRUE]
  setnames(nodes, old = names(dim_cols), new = c("x", "y"))

  midpoints = generate_midpoints(nodes, abs(adjacency_matrix))[];
  directed_edges = generate_directed_edges(nodes, adjacency_matrix)[];
  sending_receiving = assign_directed_edge_coordinates(directed_edges$sending_receiving, nodes, midpoints, multiplier = multiplier_edges);
  max_axis_marker <- max(abs(nodes$x), abs(nodes$y)) + 0.5;

  # Pulling it all together now, create shapes to plot
  # center_node_connections <- NULL;
  shapes <- NULL
  if(node_size_dynamic) {
    if(!is.null(self_connections_in_center)) {

      if(node_outer_diameter %in% c("ground", "response")) {
        wh_sr_column = c("receiver")
        if(node_outer_diameter == "ground") {
          wh_sr_column = c("sender")
        }

        # browser()
        nodes[, c("size", "outer_size", "inner_size", "shape", "color", "sign") := {
          outer <- sum(directed_edges$sending_receiving[
              directed_edges$sending_receiving[, {
                  total_connections = apply(.SD, 2, function(col) {
                    col == code
                  })
                  if (length(total_connections) == 1) {
                    return(total_connections > 0)
                  } else{
                    rowSums(total_connections) > 0
                  }
              }, .SDcols = c(wh_sr_column)]
            , (connections)]
          );
          inner <- directed_edges$self_connections[sender == as.vector(.BY$code)][, (connections)];
          # multi <- multiplier_nodes;
          this_color <- node_color;
          if(outer >= 0) {
            this_sign <- "positive"
            this_color <- hsv(colors[1,1],colors[2,1], colors[3,1])
          }
          else {
            this_sign <- "negative"
            this_color <- hsv(colors[1,2], colors[2,2], colors[3,2])
          }
          list(
            (abs(outer) + abs(inner)),
            (abs(outer)),
            (abs(inner)),
            node_shape,
            this_color,
            this_sign
          );
        }, by = c("code")]
      }
      else if (node_outer_diameter == "self") {
        nodes[directed_edges$self_connections, outer_size := connections, on = c("code" = "sender")]
      }

      # if (self_connections_in_center) {
      #   center_node_connections = assign_center_node_coordinates(
      #      center_node_connections =  directed_edges$self_connections
      #     ,nodes = nodes
      #   );
      # }
      # else {
      #   center_node_connections = assign_center_node_coordinates(
      #      center_node_connections =  directed_edges$external_connections
      #     ,nodes = nodes
      #   );
      # }
    }

    # shapes <- generate_node_shapes(center_node_connections, nodes, node_size_dynamic = node_size_dynamic, colors = colors);
  }
  else {
    nodes[ ,c("size", "outer_size", "shape", "color") := list((size * multiplier_nodes), 0, node_shape, node_color), by = c("code")]
  }

  # shapes <- lapply(shapes_df, to_shapes_list);
  nodes_dims <- as.matrix(x$rotation$nodes);
  shapes <- list()

  # Create Plot Object -----
    plot <- plot_ly(width = plot_width, height = plot_height) %>%
      layout(
         autosize = TRUE
        ,showlegend = TRUE
        ,title = title
        ,xaxis = list(type = "linear", title = "", autorange = TRUE, range = c(-max_axis_marker, max_axis_marker), scaleanchor = "y")
        ,yaxis = list(type = "linear", title = "", autorange = TRUE, range = c(-max_axis_marker, max_axis_marker))
      )

  # Allow editing of the plot ----
    # NOTE: This doesn't allow renaming/moving annotations, so may note be worth keeping
    if( !is.null(args$editable) ) {
      plot <- config( p = plot, editable = args$editable )
    }

  # Add node shapes ----
    # TODO: Revisit this, it may not be needed, with the nodes being added below
    if(!is.null(shapes)) {
      plot <- layout( p = plot ) #, shapes = shapes )
    }

  # Add points to plot ----
    # browser()
    plot <- plot_directed_points(
      plot = plot, set = x, units = units,
      with_mean = with_mean, with_points = with_points,
      scale_units = scale_points,
      colors = colors, ...
    );

  # Add edges ----
    if(with_edges == TRUE) {
      edge_shapes <- create_edge_shapes(
        plot = plot, sending_receiving, midpoints = midpoints,
        colors = colors,
        thickness_range = edge_scale,
        set = x,
        multiplier = multiplier_edges
      )
      shapes <- c(shapes, edge_shapes)
    }

  # Add nodes and their labels ------
    if(with_nodes == TRUE) {
      node_shapes <- create_node_shapes(plot = plot, nodes = nodes, colors = node_color, set = x, multiplier = multiplier_nodes)
      shapes <- c(shapes, node_shapes)
      # plot <- add_markers(
      #    p = plot
      #   ,x = ~x
      #   ,y = ~y
      #   ,data = nodes
      #   ,name = "Nodes (response)"
      #   ,legendgroup = "nodes-received"
      #   ,hovertemplate = paste('(%{x}, %{y})<extra></extra>')
      #   ,marker = list(
      #      size = ~abs(size) * multiplier_nodes
      #     ,color = node_color
      #     ,symbol = node_shape
      #     ,opacity = 0.9
      #     ,sizemode = rep("diameter", nrow(nodes))
      #     ,sizeref =  rep(1, nrow(nodes))
      #     ,sizemin =  rep(1, nrow(nodes))
      #    )
      #   ,text = ~code
      #   ,legendgroup = "nodes"
      # ) %>% add_markers(
      #   ,x = ~x
      #   ,y = ~y
      #   ,data = nodes
      #   ,name = "Nodes (self)"
      #   ,legendgroup = "nodes-self"
      #   ,hovertemplate = paste('(%{x}, %{y})<extra></extra>')
      #   ,marker = list(
      #      size = ~abs(inner_size) * multiplier_nodes
      #     ,color = ~color
      #     ,symbol = node_shape
      #     ,opacity = 0.9
      #     ,sizemode = rep("diameter", nrow(nodes))
      #     ,sizeref =  rep(1, nrow(nodes))
      #     ,sizemin =  rep(1, nrow(nodes))
      #    )
      #   ,text = ~code
      #   ,legendgroup = "nodes"
      # )
    }
    if(with_node_labels) {
      plot <- add_text(
         p = plot
        ,x = nodes_dims[, dimensions[1]]
        ,y = nodes_dims[, dimensions[2]]
        ,text = nodes$code
        ,showlegend = FALSE
        ,legendgroup = "nodes"
        ,textposition = "top right"
      )
    }

  # Custom centroids ----
    ## TODO: Deprecate
    if(!is.null(x$centroid_coordinates)) {
      if (only_sending_centroid) {
        plot <- plot_centroids(plot, x$centroid_coordinates);
      }
      else {
        plot <- plot_centroid_vectors(plot, x$centroid_coordinates);
      }
    }

  # Add shapes ----
    plot <- layout(p = plot, shapes = shapes)

  # Save plot to model ----
    x$plots[[length(x$plots) + 1]] <- plot;

  # Finish and return ----
    if(args$print_plot == TRUE) {
      x$`_plot_op` = TRUE
      return(plot);
      # invisible(x);
    }
    else {
      return(x);
    }

}

color_saturation = function(connections) {
  normalized_range = c(min(abs(connections)), max(abs(connections)));
  return(
    rescale(x=abs(connections), from=normalized_range, to=c(0.2,1))
  );
}

network_color = function(network_element, colors = COLORS_HSV[, 1:2]) {
  if (network_element$connections >= 0) {
    color = colors[,1];
  }
  else {
    color = colors[,2];
  }

  # return(hsv(color[1], network_element$saturation, color[3]));
  return(hsv(color[1], color[2] * network_element$saturation, color[3]));
}

generate_outer_node_shapes = function(nodes, center_node_connections) {
  color = COLORS_HSV[, 1];
  if (is_empty(center_node_connections)) {
    outer_nodes_df = nodes[, .(code, size, x0, x1, y0, y1)][, node_proportions := 1]
  }
  else {
    outer_nodes_df = merge(
       nodes[, .(code, size)]
      ,center_node_connections[, .(sender, connections, x0, x1, y0, y1)]
      ,by.x="code"
      ,by.y="sender"
    )[,node_proportions := (size-abs(connections))/size];
  }

  # outer_node_shapes = outer_nodes_df[, .(
  #    xref = 'x', x0, x1
  #   ,yref = 'y', y0, y1
  #   ,type = 'circle'
  #   ,fillcolor = hsv(color[1], node_proportions, color[3])
  #   ,opacity = 1
  #   ,text = size
  #   ,line = list(width = 1, color = "black")
  # )] %>% split(nodes$code) %>% unname();
  outer_node_shapes <- lapply(outer_nodes_df$code, function(cd) {
    code <- outer_nodes_df[code == cd,]
    list(
         xref = 'x', x0 = code$x0, x1 = code$x1
        ,yref = 'y', y0 = code$y0, y1 = code$y1
        ,type = 'circle'
        ,fillcolor = "black" # hsv(color[1], code$node_proportions, color[3])
        ,opacity = 1
        ,text = code$size
        ,line = list(width = 1, color = "black")
    )
  })
  return(outer_node_shapes);
}

# Default: there is no outer_node_shapes; nodes are plotted just as black dots
generate_node_shapes <- function(
  center_node_connections,
  nodes = NULL,
  node_size_dynamic = FALSE,
  center_diameter = "self",    #c("self", "response", "ground")
  total_diameter = "response", #c("self", "response", "ground")
  colors = COLORS_HSV[, 1:2]
) {
    outer_node_shapes = list();
    if (is_empty(center_node_connections)) {
      outer_node_shapes <- generate_outer_node_shapes(nodes, NULL)
      return (outer_node_shapes)
    }
    else {
      center_node_shapes = list();
      browser()
      for( i in 1:length(center_node_connections$connections) ) {
        row = center_node_connections[i, ];
        center_node_shapes[[i]] = data.table(
           xref = 'x', x0 = row$x0, x1 = row$x1
          ,yref = 'y', y0 = row$y0, y1 = row$y1
          ,type='circle'
          ,fillcolor = network_color(row, colors = colors)
          ,text = abs(row$connections)
          ,hoveron = 'fills'
          ,hoverinfo = 'text'
          #,line = list(color='black', width=3)
        )
      }
      return(append(outer_node_shapes, center_node_shapes));
    }
  # }

}

to_shapes_list = function(shape_row) {
  if (is_empty(shape_row)) { return() };
  shape_list = as.list(shape_row)
  shape_list[['line']] = list(width=1, color='black')
  return(shape_list);
}

triangle_path <- function(send, recv, width, multiplier = 1, adjuster = 600) {
  x0 <- recv$x;
  y0 <- recv$y;
  x1 <- send$x;
  y1 <- send$y;

  slope <- orthogonal_slope(x0, y0, x1, y1);
  adjustment <- adjuster / multiplier;

  sine   <- sin_theta_two_nodes(send, recv);
  cosine <- cos_theta_two_nodes(send, recv);

  half_width <- width / 2;
  width_adj <- half_width / adjustment;

  # x0 = send$x - distance_from_endpoint_center*sin_theta_two_nodes(send, receive);
  # y0 = send$y + distance_from_endpoint_center*cos_theta_two_nodes(send, receive);
  # x1 = send$x + distance_from_endpoint_center*sin_theta_two_nodes(send, receive);
  # y1 = send$y - distance_from_endpoint_center*cos_theta_two_nodes(send, receive);
  path <- paste0("M ", x0, " ", y0, " ",
                 "L ", (x1 + (width_adj * sine)), " ", (y1 - (width_adj * cosine)), " ",
                 "L ", (x1 - (width_adj * sine)), " ", (y1 + (width_adj * cosine)), " ",
                 "Z")

  return(path)
}

create_node_shapes <- function(
  plot,
  nodes,
  colors = COLORS_HSV[, 1:2],
  set = NULL,
  multiplier = 1,
  adjuster = 600
) {
  # Shape definition ----
    shape = list(
      type = "circle",
      xanchor = "x",
      yanchor = "y",
      xref = "x",
      yref = "y",
      fillcolor = "black",
      line = list (color = "black")
    )

  # Generate shapes ----
    shapes = list()
    # browser()
    adjustment <- adjuster / multiplier;
    # nodes$size <- c(0.1066974, 0.1280369, nodes$size[3], 0)
    for(i in 1:nrow(nodes)) {
      node = nodes[i, ];

      # node_color <- node$color;
      this_shape <- shape;
      half_width <- node$size / 2;
      width_adj <- half_width / adjustment;

      this_shape$x0 <- as.numeric(node$x) - width_adj
      this_shape$y0 <- as.numeric(node$y) - width_adj
      this_shape$x1 <- as.numeric(node$x) + width_adj
      this_shape$y1 <- as.numeric(node$y) + width_adj
      shapes[[length(shapes) + 1]] <- this_shape

      if(!is.null(node$inner_size) && node$inner_size > 0) {
        this_shape <- shape;
        node_color <- node$color;
        half_width <- node$inner_size / 2;
        width_adj <- half_width / adjustment;
        this_shape$fillcolor <- node_color;
        this_shape$line$color <- node_color;
        this_shape$x0 <- as.numeric(node$x) - width_adj
        this_shape$y0 <- as.numeric(node$y) - width_adj
        this_shape$x1 <- as.numeric(node$x) + width_adj
        this_shape$y1 <- as.numeric(node$y) + width_adj
        shapes[[length(shapes) + 1]] <- this_shape
      }
    }

  # Return updated plot ----
    # plot <- layout(p = plot, shapes = shapes)
    # return(plot);
    return(shapes);
}

create_edge_shapes = function(
  plot,
  edges, midpoints,
  colors = COLORS_HSV[, 1:2],
  thickness_range = c(ifelse(min(edges$connections) == 0, 0, 0.1), 1),
  set = NULL,
  multiplier = 1,
  adjustment = 600
) {
  # Setup ----
    # thickness = c(min(abs(edges$connections)), max(abs(edges$connections)));
    opacity = rescale(x=abs(edges$connections), from=c(0,1), to = thickness_range);

    # message("Note 5: I think we may want the saturation to be multiplied by color[2]")

  # Hidden point ----
  ### - specifically for use with the legend
    plot = add_trace(
         plot
        ,legendgroup = "network"
        ,type = "scatter"
        ,hoverinfo = NULL
        ,mode = "markers"
        ,name = "Network"
        ,opacity = 1
        ,x = 0
        ,y = 0
        ,marker = list ( color = "#000", size = 0.1 )
    )

  # Edge Shapes ----
    shape = list(
      type = "path",
      xanchor = "x",
      yanchor = "y",
      xref = "x",
      fillcolor = "black",
      yref = "y",
      line = list (color = "black")
    )

    shapes = list()
    for(i in 1:nrow(edges)) {
      row = edges[i, ];
      edge_color <- network_color(row, colors = colors);
      send <- set$rotation$nodes[code == row$sender]
      recv <- set$rotation$nodes[code == row$receiver];
      mid  <- midpoints[sending %in% c(send$code, recv$code) & receiving %in% c(send$code, recv$code)]
      if(nrow(mid) == 0) {
        midpoint <- list(x = recv$SVD1, y = recv$SVD2);
      } else {
        midpoint <- list(x = mid$midpoint_x, y = mid$midpoint_y);
      }
      # print(row$connections)
      this_shape <- shape;
      this_shape$path <- triangle_path(
        list(x=send$SVD1, y=send$SVD2),
        midpoint,
        width = row$connections,
        multiplier = multiplier,
        adjuster = adjustment
      )
      this_shape$line$color <- edge_color;
      this_shape$fillcolor <- edge_color;
      this_shape$opacity <- opacity[i]
      this_shape$name <- abs(row$connections)
      this_shape$legendgroup = "network"

      shapes[[length(shapes) + 1]] <- this_shape
      # plot = add_trace(
      #    plot
      #   ,fill = 'toself'
      #   ,fillcolor = edge_color
      #   ,legendgroup = "network"
      #   ,hoverinfo = 'text'
      #   ,hoveron = 'fills'
      #   ,marker = list(opacity = 0) # don't plot midpoints
      #   ,mode = "markers"
      #   ,opacity = opacity[i]
      #   ,text = abs(row$connections)
      #   ,name = NULL
      #   ,type = 'scatter'
      #   ,showlegend = FALSE
      #   ,x = as.numeric(row[, grep(x = colnames(row), pattern = "^x_path"), with = FALSE])
      #   ,y = as.numeric(row[, grep(x = colnames(row), pattern = "^y_path"), with = FALSE])
      # )
    }

    # plot <- layout(p = plot, shapes = shapes)

  # Return updated plot ----
    # return(plot);
    return(shapes);
}

plot_centroids = function(plot, centroid_info) {
  for (i in 1:length(centroid_info$x)) {
    color = centroid_info$color[[i]];
    plot = add_trace(
      plot
      , x= centroid_info$x[i]
      , y = centroid_info$y[i]
      , text='centroid'
      , marker = list(
        color = color
      )
    )
  }
  return(plot);
}

plot_directed_points <- function(
  plot, set, units,
  labels = names(units), dimensions = 1:2,
  scale_units = TRUE,
  opacity_range = c(0.25,0.75),
  with_mean = TRUE,
  with_points = TRUE,
  just_vectors = TRUE,
  colors = COLORS_HSV
) {

  for(i in seq(length(units))) {
    # browser()
    point_color_hsv <- colors[, i];
    point_color <- hsv(point_color_hsv[1], point_color_hsv[2], point_color_hsv[3])
    wh_units <- units[[i]]
    is_mean <- length(wh_units) > 1
    resp <- as.matrix(set$points[wh_units & ENA_DIRECTION == "response",])
    grnd <- as.matrix(set$points[wh_units & ENA_DIRECTION == "ground",])

    resp_pts <- resp[, dimensions, drop = FALSE]
    grnd_pts <- grnd[, dimensions, drop = FALSE]
    scaleFactor = 1.0

    if( scale_units == TRUE ) {
      np.min.x = abs(min(as.matrix(set$rotation$nodes)[, dimensions[1]]));
      np.min.y = abs(min(as.matrix(set$rotation$nodes)[, dimensions[2]]));
      rp.min.x = abs(min(as.matrix(set$points)[, dimensions[1]]));
      rp.min.y = abs(min(as.matrix(set$points)[, dimensions[2]]));
      maxMin = abs(max(np.min.x / rp.min.x, np.min.y / rp.min.y));

      np.max.x = abs(max(as.matrix(set$rotation$nodes)[, dimensions[1]]));
      np.max.y = abs(max(as.matrix(set$rotation$nodes)[, dimensions[2]]));
      rp.max.x = abs(max(as.matrix(set$points)[, dimensions[1]]));
      rp.max.y = abs(max(as.matrix(set$points)[, dimensions[2]]));
      maxMax = abs(max(np.max.x / rp.max.x, np.max.y / rp.max.y));
      scaleFactor = min(maxMin, maxMax);
    }
    else if (is.numeric(scale_units)) {
      scaleFactor = scale_units;
    }
    # browser()
    resp_pts <- resp_pts * scaleFactor
    grnd_pts <- grnd_pts * scaleFactor

    pts_df <- data.frame(x = c(grnd_pts[,1], resp_pts[,1]), y = c(grnd_pts[,2], resp_pts[,2]));
    distances <- sapply(seq(nrow(grnd_pts)), function(r) dist(rbind(grnd_pts[r,], resp_pts[r,])));
    arrow_colors <- matrix(rep(point_color_hsv, length(distances)), ncol = 3, byrow = TRUE);

    # arrow_colors[, 2] <- arrow_colors[, 2] * (distances);
    # opacity <- scales::rescale(x = distances, to = c(0, point_color_hsv[2]))
    opacity <- scales::rescale(x = distances, to = c(opacity_range[1], opacity_range[2]))

    # arrow_color_hex <- apply(arrow_colors, 1, function(x) hsv(x[1], x[2], x[3]));
    arrow_color_hex <- sapply(seq(nrow(arrow_colors)), function(i) {
      x <- arrow_colors[i, ];
      hsv(x[1], x[2], x[3], opacity[i])
    });
    if(with_points == TRUE) {
      if(just_vectors == FALSE) {
        plot <- add_markers(
           plot
          ,x = pts_df$x
          ,y = pts_df$y
          ,text = labels[i]
          ,legendgroup = "points"
          ,name = paste0(labels[i], " Points")
          ,marker = list(
             color = point_color
            ,size = ~size * multiplier
            ,sizemode = rep("diameter", nrow(pts_df))
            ,sizeref =  rep(1, nrow(pts_df))
            ,sizemin =  rep(1, nrow(pts_df))
          )
        )
      }
      plot <- add_annotations(
         plot
        ,ax = grnd_pts[,1]
        ,x = resp_pts[,1]
        ,axref = "x"
        ,ay = grnd_pts[,2]
        ,y = resp_pts[,2]
        ,ayref = "y"
        ,showarrow = TRUE
        ,text = ''
        ,xref = 'x'
        ,yref = 'y'
        ,arrowcolor = arrow_color_hex
        ,legendgroup = "points"
      )
    }

    if(with_mean == TRUE && is_mean == TRUE && nrow(grnd_pts) > 1) {
      mean_color <- hsv(point_color_hsv[1], point_color_hsv[2], point_color_hsv[3], alpha = 1.0);
      plot <- add_markers(
         plot
        ,x = c(mean(grnd_pts[, 1])) #, mean(resp_pts[, 1]))
        ,y = c(mean(grnd_pts[, 2])) #, mean(resp_pts[, 2]))
        # ,x = c(mean(resp_pts[, 1])) #, mean(resp_pts[, 1]))
        # ,y = c(mean(resp_pts[, 2])) #, mean(resp_pts[, 2]))
        ,text = paste0(labels[i], " Mean")
        ,marker = list( color = mean_color, symbol = "square" )
        ,legendgroup = "means"
        ,name = paste0(labels[i], " Mean")
      )
      # browser()
      plot <- add_annotations(
         plot
        ,ax = mean(grnd_pts[,1])
        ,ay = mean(grnd_pts[,2])
        ,x = mean(resp_pts[,1])
        ,y = mean(resp_pts[,2])
        ,axref = "x"
        ,ayref = "y"
        ,showarrow = TRUE
        ,text = ''
        # ,xref = 'x'
        # ,yref = 'y'
        ,arrowcolor = mean_color
        ,legendgroup = "means"
      )
      # plot <- add_annotations(
      #    plot
      #   ,ax = 0
      #   ,x = 1
      #   ,axref = "x"
      #   ,ay = 0
      #   ,y = 1
      #   ,ayref = "y"
      #   ,showarrow = TRUE
      #   ,text = ''
      #   # ,xref = 'x'
      #   # ,yref = 'y'
      #   ,arrowcolor = arrow_color_hex
      #   ,legendgroup = "means"
      # )
    }
  }
  return(plot)
}

plot_centroid_vectors = function(plot, vector_data) {
  for (i in 1:length(vector_data)) {
    data = vector_data[i] %>% unlist;
    if (data['sending_x'] == data['receiving_x'] && data['sending_y'] == data['receiving_y']) {
      plot = add_trace(
         plot
        ,x = as.numeric(data['sending_x'])
        ,y = as.numeric(data['sending_y'])
        ,text = 'centroid'
        ,marker = list( color = data['color'])
      )
    }
    else {
      plot = add_annotations(
         plot
        ,ax = as.numeric(data['sending_x'])
        ,x = as.numeric(data['receiving_x'])
        ,axref = "x"
        ,ay = as.numeric(data['sending_y'])
        ,y = as.numeric(data['receiving_y'])
        ,ayref = "y"
        ,showarrow = TRUE
        ,text = ''
        ,xref = 'x'
        ,yref = 'y'
        ,arrowcolor = data['color']
      )
    }
  }
  return(plot);
}

to_square <- function(x) {
  n <- sqrt(length(x))

  dim(x) <- c(n, n);
  x
}

#' Re-class vector as ena.co.occurrence
#'
#' @param x Vector to re-class
#'
#' @return re-classed vector
#' @export
as.ena.co.occurrence <- function(x) {
  if(is.factor(x)) {
    x = as.character(x)
  }
  class(x) = c("ena.co.occurrence", class(x))
  x
}

#' Title
#'
#' @param x directed ENA model
#'
#' @return matrix
#' @export
as.directed.matrix <- function(x) {
  return(x);
}

#' Title
#'
#' @param x directed ENA model
#' @param names character
#'
#' @return square matrix
#' @export
as.square <- function(x, names = NULL) {
  UseMethod("as.square", x);
}

#' Title
#'
#' @param x directed ENA model
#' @param ... unused
#'
#' @return square matrix
#' @export
as.square.ena.connections <- function(x, ...) {
  x_ <- as.matrix(x);
  square_matrices <- lapply(seq(nrow(x_)), function(y) {
    y <- as.ena.matrix(to_square(x_[y, ]), "ena.directed.connections")
  })
  square_matrices
}

#' Title
#'
#' @param x directed ENA model
#' @param names character
#'
#' @return square matrix
#' @export
as.square.matrix <- function(x, names = NULL) {
  x_ <- as.data.frame(as.matrix(x))
  as.square(x_)
}

#' Title
#'
#' @param x vector
#' @param names character vector
#'
#' @return square data.frame
#' @export
as.square.data.frame <- function(x, names = NULL) {
  x_ <- as.matrix(x)
  n <- sqrt(length(x_))

  dim(x_) <- c(n, n);
  df <- as.data.frame(x_)

  if(is.null(names)) {
    names <- paste("V", seq(n), sep = "")
  }
  rownames(df) <- colnames(df) <- names;
  df <- as.ena.matrix(df, "ena.directed.connections")
  df
}

#' Convert object to a directed set
#'
#' @param x Object to convert
#'
#' @return list as ena.directed.set
#' @export
as.ena.directed.set <- function(x) {
  to <- c("ena.directed.set")
  if(!inherits(x = x, what = "ena.set")) {
    to <- c(to, "ena.set")
  }
  class(x) <- c(to, class(x));
  return(x);
}

#' Print a directed ena model
#'
#' @param x model to print
#' @param ... not implemented
#' @param plot logical if TRUE, print the plots (default is FALSE)
#'
#' @export
print.ena.set.directed <- function(x, ..., plot = FALSE) {
  x.unclass <- unclass(x)

  if(
    !is.null(x.unclass$`_plot_op`) &&
    x.unclass$`_plot_op` == T
  ) {
    base::print(x.unclass$plots)
  }
  else {
    if(plot == FALSE) {
      x.unclass$plots <- NULL
    }
    base::print(x.unclass)
  }
}


is_empty <- function(x) {
  return(is.null(x) || nrow(x) == 0 || ncol(x) == 0);
}

to_datatable = function(appended_vectors, ncol, nrow) {
  # because we're representing the data as a vector of row vectors, we need to
  # reverse the matrix ncol/nrow and transpose the matrix to get the correct output
  return(data.table(t(matrix(appended_vectors, nrow = ncol, ncol = nrow))));
}

two_point_weighted_average = function(point_1, point_2, relative_weight) {
  # distance formula from:
  # https://www.khanacademy.org/partner-content/pixar/environment-modeling-2/mathematics-of-parabolas2-ver2/v/weighted-average-two-points
  return((1-relative_weight)*point_1 + (relative_weight*point_2));
}

orthogonal_slope = function(receiver_x, receiver_y, sender_x, sender_y) {
  sender_to_receiver_slope = (sender_y - receiver_y)/(sender_x - receiver_x);
  return(-(1/sender_to_receiver_slope));
}

cos_theta_two_nodes = function (sender, receiver) {
  hypotenus = sqrt((receiver$x - sender$x)^2 + (receiver$y - sender$y)^2);
  adjacent = receiver$x - sender$x;
  return(adjacent/hypotenus);
}

sin_theta_two_nodes = function (sender, receiver) {
  hypotenus = sqrt((receiver$x - sender$x)^2 + (receiver$y - sender$y)^2);
  opposite = receiver$y - sender$y;
  return(opposite/hypotenus);
}

is_vertical_line = function(slope) {
  return(abs(slope) == Inf);
}

is_horizontal_line = function(slope) {
  return(is.nan(slope));
}

format_adjacency_matrix = function(raw_data) {
  # node_names = raw_data$node_names;
  # nrows = length(node_names);
  # return(function(raw_matrix) {
  #   matrix = as.data.frame(to_datatable(raw_matrix, nrows, nrows));
  #   colnames(matrix) = rownames(matrix) = node_names;
  #   return(matrix);
  # })
  node_names <- raw_data$node_names;
  node_len <- length(node_names);
  lapply(raw_data$raw_data, function(x) {
    as.data.frame(matrix(x, ncol = node_len, nrow = node_len, dimnames = list(node_names, node_names), byrow = TRUE))
  })
}

raw_coordinates_to_datatable = function(names) {
  nrows = length(names);
  pts.circle <- t(sapply(1:nrows, function(r) c(cos(2*r*pi/nrows), sin(2*r*pi/nrows))))
  colnames(pts.circle) = c('x', 'y');

  # ncols = 3; # 3 columns, nodename, x, and y
  # nodes = to_datatable(raw_data$coordinates, ncols, nrows);
  nodes <- data.table("name" = names);
  nodes <- cbind(nodes, pts.circle);

  # nodes[, ':=' (x = as.double(x), y = as.double(y))];
  return(nodes);
}

# -------------------------------------
# DATA TRANSFORMATIONS
# -------------------------------------

generate_adjacency_matrix = function(matrices) {
  # guard clauses
  if (length(matrices) > 2) stop("cannot pass this function more than 2 adjacency matrices");
  if (typeof(matrices) != 'list') stop('raw data for adjacency matrix must be in a list')

  # business logic
  if ( length(matrices) == 2 ) {
    # message("Note 4: need to remove reduce, use Reduce?")
    adjacency_matrix = reduce(matrices, `-`);
  }
  else {
    adjacency_matrix = matrices[[1]];
  }

  return(adjacency_matrix);
}

total_node_connections_calculator = function(adjacency_matrix) {
  # return "partially applied" anonymous function in order to generate column values
  return(
    function(node_name) {
      # total node connections should just be sending connections
      total = rowSums(adjacency_matrix[node_name,]);
      # HACK: make sure node is not invisible, set minimum to 1;
      if (total <= 0) { total = 1 };
      return(total);
    }
  )
}

max_node_diameter = function(nodes){
  pts <- as.matrix(nodes)
  plot_area = (max(pts[,1]) - min(pts[,1])) * (max(pts[,2]) - min(pts[,2]));
  # assume that multiplying by the max scalar gives us the area of a square, which should encompass
  # the outer bounds of the largest plotted node.
  # take the square root to get the base of the square, which we will use as the max
  # diameter size, translated to our plot's cartesian plane.
  # all other nodes and edges will be pegged to this max node radius
  return(sqrt(MAX_NODE_DIAMETER_SCALAR * plot_area));
}

generate_node_coordinates = function(codes, connections_calculator = NULL) {
  nodes = raw_coordinates_to_datatable(codes);
  node_radius = max_node_diameter(nodes) / 2;

  if(!is.null(connections_calculator)) {
    nodes[, size := sapply(name, connections_calculator)];
  }
  else {
    nodes[, size := 1];
  }

  max_size = max(nodes$size);
  nodes[, ':=' (relative_size = size / max_size) ];
  nodes[, ':=' (node_radius = relative_size * node_radius) ];

  nodes[, ':=' (
     x0 = x - node_radius
    ,x1 = x + node_radius
    ,y0 = y - node_radius
    ,y1 = y + node_radius
  )];
  return(nodes);
};

# NOTE: technically, these labels are somewhat misleading since this matrix represents
# the total edge weights of both sender and receiver. However, in order to identify
# the weighted midpoint between two nodes, we arbitrarily assign one node to be a
# sender and another to be a receiver.
# This allows us to calculate the proportion of edge weight from the sender to receiver out
# of the total connections between the two nodes.
generate_midpoints = function(nodes, adjacency_matrix) {
  midpoints = data.table((t(combn(nodes$code, 2))), stringsAsFactors = FALSE);
  colnames(midpoints) = c('sending', 'receiving');

  # browser()
  midpoints[, c("sender_receiver_connections","total_connections", "proportion", "midpoint_x", "midpoint_y") := {
    dim_cols <- which(find_dimension_cols(nodes))
    send <- adjacency_matrix[sending, receiving];
    recv <- adjacency_matrix[receiving, sending];
    prop <- send / (send + recv)
    pts_s <- as.numeric(nodes[code == sending, dim_cols, with = FALSE])
    pts_r <- as.numeric(nodes[code == receiving, dim_cols, with = FALSE])

    mid_x <- two_point_weighted_average( point_1 = pts_s[1], point_2 = pts_r[1], relative_weight = prop);
    mid_y <- two_point_weighted_average( point_1 = pts_s[2], point_2 = pts_r[2], relative_weight = prop);

    list(send, send + recv, prop, mid_x, mid_y);
  }, by = c("sending", "receiving")];

  # this will also remove rows where proportion is NaN
  return(midpoints[proportion != 0]);
}

generate_directed_edges = function(nodes, adjacency_matrix) {
  # directed_edges = data.table(
  #   permutations(n=length(nodes$name), r=2, v=nodes$name, repeats.allowed = TRUE)
  # )
  cols_sr <- c("sender", "receiver");
  directed_edges <- data.table(
    expand.grid("receiver"=nodes$code[order(nodes$code)],
                "sender"=nodes$code[order(nodes$code)], stringsAsFactors = FALSE)
  )[, cols_sr, with = FALSE];
  directed_edges[, connections := adjacency_matrix[sender, receiver], by = cols_sr];
  directed_edges[, saturation := color_saturation(connections)]

  sending_receiving = copy(directed_edges[sender != receiver & abs(connections) > 0]);
  self_connections = copy(directed_edges[sender == receiver]);
  external_connections = copy(directed_edges[sender == receiver]);

  if (!is_empty(external_connections)) {
    # message("Note 1: `external_connections` wasn't saving the mutation, seems questionable. Also, the
    #           resulting connections were all the same.")
    # external_connections %>% mutate(connections = nodes[name == sender]$size - connections);
    external_connections[, connections := {
        nodes[code == sender]$size - connections
      },
      by = sender
    ]
  }
  edges_list <- list(
    sending_receiving = sending_receiving
   ,self_connections = self_connections
   ,external_connections = external_connections
  );

  edges = lapply(edges_list, function(edge_df) {
    edge_df[, relative_size := { abs(connections) / nodes[code == sender]$size }, by = sender];
    edge_df[, distance_from_endpoint_center := { relative_size * nodes[code == sender]$node_radius }, by = sender];

    return(copy(edge_df));
  });

  return(edges);
}

assign_center_node_coordinates = function(center_node_connections, nodes) {
  if(is_empty(center_node_connections)) return(center_node_connections);

  if("distance_from_endpoint_center" %in% colnames(center_node_connections)) {
    setnames(center_node_connections, old = "distance_from_endpoint_center", new = "radius")
  }

  # xy_cols <- colnames(as.matrix(nodes))[1:2];
  xy_cols <- names(which(find_dimension_cols(nodes)))

  center_node_connections <- center_node_connections[nodes,
                                                     c(colnames(center_node_connections), xy_cols),
                                                     on = c(sender = "code"),
                                                     with = FALSE][order(nodes$code)]
  # browser()
  dims <- center_node_connections[, rENA::find_dimension_cols(center_node_connections), with = F];
  center_node_connections$x0 <- dims[, 1] - nodes$node_radius
  center_node_connections$x1 <- dims[, 1] + nodes$node_radius
  center_node_connections$y0 <- dims[, 2] - nodes$node_radius
  center_node_connections$y1 <- dims[, 2] + nodes$node_radius

  center_node_connections
};

calculate_centroid_vectors = function(matrices) {
  # matrices = format_adjacency_matrix(raw_data);

  nodes = raw_coordinates_to_datatable(colnames(matrices[[1]]));
  mass_vectors = list();
  for (i in 1:length(matrices)) {
    # From above, this appears to assume each item in matrices is already a matrix, but it's still a function
    adjacency_matrix = matrices[[i]];
    total_mass = sum(adjacency_matrix);
    node_sending = rowSums(adjacency_matrix);
    all_sending = data.table(name = names(node_sending), mass = as.vector(node_sending));
    sending_mass = merge(all_sending, nodes[, .(name,x,y)], by="name")[, ':=' (
      mass_x = mass*x
      , mass_y = mass*y
    )];
    node_receiving = colSums(adjacency_matrix);
    all_receiving = data.table(name = names(node_receiving), mass = as.vector(node_receiving));
    receiving_mass = merge(all_receiving, nodes[, .(name,x,y)], by="name")[, ':=' (
      mass_x = mass*x
      , mass_y = mass*y
    )];
    color = COLORS_HSV[, i+1]
    sending_x = sum(sending_mass$mass_x) / total_mass;
    sending_y = sum(sending_mass$mass_y) / total_mass;
    receiving_x = sum(receiving_mass$mass_x) / total_mass;
    receiving_y = sum(receiving_mass$mass_y) / total_mass;

    mass_vectors[[i]] = list(
      sending_x = sending_x
      , sending_y = sending_y
      , receiving_x = receiving_x
      , receiving_y = receiving_y
      , color = hsv(color['h'], 0.93, color['v'])
    )
  }
  return(mass_vectors);
}

calculate_centroid = function(raw_data) {
  # matrices = lapply(data$raw_data, format_adjacency_matrix(data));
  matrices = format_adjacency_matrix(raw_data);

  nodes = raw_coordinates_to_datatable(raw_data);
  x_coords = c();
  y_coords = c();
  color = list();
  for (i in 1:length(matrices)) {
    centroid_color = COLORS_HSV[, i+1]
    adjacency_matrix = matrices[[i]];
    node_sending = rowSums(adjacency_matrix);
    total_mass = sum(node_sending);
    all_sending = data.table(name = names(node_sending), mass = as.vector(node_sending));
    node_mass = merge(all_sending, nodes[, .(name,x,y)], by="name")[, ':=' (
      mass_x = mass*x
      , mass_y = mass*y
    )];
    x_coords[i] = sum(node_mass$mass_x) / total_mass;
    y_coords[i] = sum(node_mass$mass_y) / total_mass;
    # #F01186
    color[[i]] = hsv(centroid_color['h'], 0.93, centroid_color['v'])
  }
  return(list(x = x_coords , y = y_coords, color = color));
}

assign_directed_edge_coordinates = function(sending_receiving, nodes, midpoints, multiplier = 1) {
  # message("Note 2: Making a copy to remove DT warning for now") # I can't find the automatic R copy (clm)
  sending_receiving <- copy(sending_receiving);
  # nodes <- copy(nodes);
  # midpoints <- copy(midpoints);
  sr_cols <- c("sender", "receiver")
  colvector <- c("endpoint_coordinates_x", "endpoint_coordinates_y")

  sending_receiving[, (colvector) := {
    midpoint_present_for_nodes <- midpoints[
      (midpoints$sending %in% .SD$sender   & midpoints$receiving %in% .SD$receiver) |
      (midpoints$sending %in% .SD$receiver & midpoints$receiving %in% .SD$sender)
    ];
    if (is_empty(midpoint_present_for_nodes)) {
      list(x = nodes[code == receiver]$x, y = nodes[code == receiver]$y);
    }
    else {
      list(x = midpoint_present_for_nodes$midpoint_x, y = midpoint_present_for_nodes$midpoint_y);
    }
  }, by = sr_cols, .SDcols = sr_cols]

  colvector_tri <- c("triangle_base_coordinates_x0", "triangle_base_coordinates_y0",
                     "triangle_base_coordinates_x1", "triangle_base_coordinates_y1")
  sending_receiving[, (colvector_tri) := {
    send = nodes[code == .SD$sender];
    receive = nodes[code == .SD$receiver];
    slope = orthogonal_slope(receive$x, receive$y, send$x, send$y);
    sine = sin_theta_two_nodes(send, receive);
    cosine = cos_theta_two_nodes(send, receive);
    distance_from_endpoint_center = distance_from_endpoint_center * multiplier;
    x0 = send$x - distance_from_endpoint_center*sin_theta_two_nodes(send, receive);
    y0 = send$y + distance_from_endpoint_center*cos_theta_two_nodes(send, receive);
    x1 = send$x + distance_from_endpoint_center*sin_theta_two_nodes(send, receive);
    y1 = send$y - distance_from_endpoint_center*cos_theta_two_nodes(send, receive);
    if (is_horizontal_line(slope)) {
      y0 = y1 = send$y;
      x0 = send$x + distance_from_endpoint_center;
      x1 = send$x - distance_from_endpoint_center;
    }
    else if (is_vertical_line(slope)) {
      x0 = x1 = send$x;
      y0 = send$y + distance_from_endpoint_center;
      y1 = send$y - distance_from_endpoint_center;
    }
    list(x0 = x0, y0 = y0, x1 = x1, y1 = y1)
  }, by = sr_cols, .SDcols = sr_cols]

  colvector_paths <- c("x_path_ep_1", "x_path_tr_1", "x_path_tr_2", "x_path_ep_2",
                       "y_path_ep_1", "y_path_tr_1", "y_path_tr_2", "y_path_ep_2")
  sending_receiving[, (colvector_paths) := {
    list(
       .SD$endpoint_coordinates_x
      ,.SD$triangle_base_coordinates_x1
      ,.SD$triangle_base_coordinates_x0
      ,.SD$endpoint_coordinates_x
      ,.SD$endpoint_coordinates_y
      ,.SD$triangle_base_coordinates_y1
      ,.SD$triangle_base_coordinates_y0
      ,.SD$endpoint_coordinates_y
    )
  }, by = sr_cols]
  return(sending_receiving[]);
}


#' Accumulation for Directed Networks
#'
#' @param x object to accumulate
#' @param units character vector of columns to use as units
#' @param conversations character vector of columns to use as conversations
#' @param codes character vector of columns to use as codes
#' @param binary logical TRUE (default) will converat row accumulations to binary
#' @param weight_using TBD
#' @param windowSize numeric integer representing window size
#' @param ... TBD
#'
#' @return ena model object
#' @export
directed_accumulation <- function(
  x,
  units, conversations, codes,
  binary = TRUE, weight_using = NULL,
  windowSize = 4, ...
) {
  #### Prepare input -----
    dat <- NULL
    set <- NULL

    if(inherits(x = x, what = "data.frame")) {
      set <- ena.set.directed(x, units, conversations, codes);
    }
    else if(inherits(x = x, what = "ena.set")) {
      set <- x;
    }
    else {
      stop("`x` must be an existing ena.set.directed or a data.frame or data.table in which to create a set from")
    }

    dat <- data.table::copy(set$model$raw.input);

  #### Calcuate row.connection.counts -----
    all.cols <- c(units, codes)
    new.cols <- paste("V", seq(length(codes) ^ 2), sep = "")
    row.connection.counts <- dat[,
      (new.cols) := accum_stanza_window(df = .SD[, .SD, .SDcols = codes], windowSize = windowSize, binary = binary),
      by = conversations,
      .SDcols = all.cols
    ]
    for(i in new.cols) {
      set(row.connection.counts, j = i, value = as.ena.co.occurrence(row.connection.counts[[i]]))
    }

  #### Set the column classes -----
    for(i in colnames(row.connection.counts)) {
      if(i %in% c(units, conversations))
        set(row.connection.counts, j = i, value = rENA::as.ena.metadata(row.connection.counts[[i]]))
      else if (i %in% c(codes))
        set(row.connection.counts, j = i, value = rENA:::as.ena.code(row.connection.counts[[i]]))
    }
    row.connection.counts <- as.ena.matrix(x = row.connection.counts, "row.connections")

  #### Additional columns added -----
    row.connection.counts$ENA_UNIT <- rENA::as.ena.metadata(rENA::merge_columns_c(df = row.connection.counts, cols = units, sep = "."));

  #### Calculate connection.counts (Sum rows by unit)
    connection.counts <- row.connection.counts[, lapply(.SD, sum), .SDcols = new.cols, by = c("ENA_UNIT", units)]
    for(i in which(colnames(connection.counts) %in% c("ENA_UNIT", units))) {
      set(connection.counts, j = i, value = rENA::as.ena.metadata(connection.counts[[i]]))
    }
    connection.counts <- as.ena.matrix(x = connection.counts, "ena.connections")

  ##### Check weighting of accumulation ----
    if( binary == FALSE && !is.null(weight_using) ) {
      browser()
    }

  #### Store values on the ena.set object -----
    set$connection.counts <- connection.counts;
    set$model$row.connection.counts <- row.connection.counts;

    set$meta.data <- connection.counts[,rENA::find_meta_cols(connection.counts), with = FALSE];
    set$`_function.params`$binary <- binary
    set$`_function.params`$windowSize <- windowSize

  #### Done -----
    return(set)
}

#' Normalize and optimize a directed accumulation
#'
#' @param x a ena.set.directed with a set of accumulated adjacency matrices
#' @param norm.by TBD
#' @param rotate.using TBD
#' @param rotation.params TBD
#' @param rotation.set TBD
#' @param rotation_on TBD
#' @param optimize_on TBD
#' @param node.position.method TBD
#' @param ... TBD
#'
#' @return ena.set.directed with calcualted points and node positiosn
#' @export
directed_model <- function(
  x, ...,
  norm.by = rENA::fun_sphere_norm,
  node.position.method = directed_node_optimization,
  rotate.using = rENA::ena.svd,
  rotation.params = NULL, rotation.set = NULL,
  rotation_on = "response", #c("ground", "response"),
  optimize_on = c("ground", "response")
) {

  ##### Generate data.frame with a ground and response row for each unit -----
    ##### > Prep ----
      meta_cols <- which(rENA::find_meta_cols(x$connection.counts))

    ##### > Generate ----
      df <- data.table::rbindlist(lapply(seq(nrow(x$connection.counts)), function(r) {
        row <- x$connection.counts[r, ];
        meta <- row[, meta_cols, with = FALSE]
        mat <- to_square(as.matrix(row))

        ground <- as.vector(t(mat))
        response <- as.vector(mat)

        df <- data.frame(matrix(c(ground, response), nrow = 2, byrow = TRUE, dimnames = list(NULL, colnames(row)[find_code_cols(row)])))
        df <- cbind(meta, ENA_DIRECTION = c("ground", "response"), df)
        df
      }))
      # setorder(x = df, cols = "ENA_DIRECTION")

    ##### > Set column classes ----
      class(df) <- class(x$connection.counts);
      for(i in colnames(df)) {
        if(i %in% c(names(meta_cols), "ENA_DIRECTION"))
          set(df, j = i, value = rENA::as.ena.metadata(df[[i]]))
        else
          set(df, j = i, value = as.ena.co.occurrence(df[[i]]))
      }

  ##### Model preperation -----
    all_meta_cols <- df[, find_meta_cols(df), with = FALSE]
    code_columns <- colnames(df)[rENA::find_code_cols(df)];

  ##### Normalize the ground/response adjacency vectors ----
    ### > Normalize adjacency vectors ----
      line.weights <- norm.by(as.matrix(df));

    ### > Classify objects ----
      colnames(line.weights) <- code_columns

      line.weights.dt <- as.data.table(line.weights)
      for (i in seq(ncol(line.weights.dt))) {
        set(line.weights.dt, j = i, value = as.ena.co.occurrence(line.weights.dt[[i]]))
      }

      x$line.weights <- cbind(all_meta_cols, line.weights.dt)
      for(i in which(!find_code_cols(x$line.weights))) {
        set(x$line.weights, j = i, value = rENA::as.ena.metadata(x$line.weights[[i]]))
      }
      class(x$line.weights) <- c("ena.line.weights", class(x$line.weights))


  ##### Center the Normalized data ----
    ### > Center the line.weights ----
      points.for.projection <- center_data_c(line.weights)

    ### > Classify objects ----
      colnames(points.for.projection) <- code_columns;
      x$model$points.for.projection = as.data.table(points.for.projection)

      for (i in seq(ncol(x$model$points.for.projection))) {
        set(x$model$points.for.projection, j = i, value = as.ena.co.occurrence(x$model$points.for.projection[[i]]))
      }
      x$model$points.for.projection <- as.ena.matrix(cbind( all_meta_cols, x$model$points.for.projection), "ena.points")

      for(i in which(!find_code_cols(x$model$points.for.projection))) {
        set(x$model$points.for.projection, j = i, value = rENA::as.ena.metadata(x$model$points.for.projection[[i]]))
      }

  ##### SVD -----

    ### > New rotation ----
      if (!is.null(rotate.using) && is.null(rotation.set)) {
        ### >> Row preparation ----
          rotation_rows <- x$model$points.for.projection$ENA_DIRECTION %in% rotation_on
          # resp_rows <- !grnd_rows
          # rotation <- do.call(rotate.using, list(x, rotation.params))
          # rotation_ground <- prcomp(points.for.projection[grnd_rows,], retx=FALSE, scale=FALSE, center=FALSE, tol=0)
          # rotation_response <- prcomp(points.for.projection[resp_rows,], retx=FALSE, scale=FALSE, center=FALSE, tol=0)

          pts = as.matrix(x$model$points.for.projection)[rotation_rows,]

        ### >> Perform SVD on normed, centered data ----
          pcaResults = prcomp(pts, retx=FALSE, scale=FALSE, center=FALSE, tol=0)

        ### >> Classify Results ----
          colnames(pcaResults$rotation) = c(
            paste('SVD',as.character(1:ncol(pcaResults$rotation)), sep='')
          );

          # rotationSet = ENARotationSet$new(rotation = pcaResults$pca, codes = enaset$codes, node.positions = NULL, eigenvalues = pcaResults$latent)
          rotation = ENARotationSet$new(
            rotation = pcaResults$rotation,
            codes = x$codes,
            node.positions = NULL,
            eigenvalues = pcaResults$sdev^2
          )

          x$rotation.matrix <- as.data.table(rotation$rotation, keep.rownames = "codes")
          for (i in seq(ncol(x$rotation.matrix))) {
            if(i == 1) {
              set(x$rotation.matrix, j = i, value = as.ena.metadata(x$rotation.matrix[[i]]))
            } else {
              set(x$rotation.matrix, j = i, value = rENA:::as.ena.dimension(x$rotation.matrix[[i]]))
            }
          }
          class(x$rotation.matrix) <- c("ena.rotation.matrix", class(x$rotation.matrix))

          x$rotation$rotation.matrix <- x$rotation.matrix
          x$rotation$eigenvalues <- rotation$eigenvalues;
      }
    ### > Custom rotation provided ----
      else if (!is.null(rotation.set)) {
        if (is(rotation.set, "ena.rotation.set")) {
          x$rotation.matrix <- rotation.set$rotation.matrix
          x$rotation$rotation.matrix <- rotation.set$rotation.matrix
          x$rotation$nodes <- rotation.set$nodes;
          x$rotation$eigenvalues <- rotation.set$eigenvalues
        }
        else {
          stop("Supplied rotation.set is not an instance of ENARotationSet")
        }
      }
      else {
        stop("Unable to find or create a rotation set")
      }

  ##### Generate the rotated points ----
    if (!is.null(x$rotation.matrix)) {
      points <- points.for.projection %*% as.matrix(x$rotation.matrix)
      points.dt <- as.data.table(points)
      for (i in seq(ncol(points.dt))) {
        set(points.dt, j = i, value = rENA:::as.ena.dimension(points.dt[[i]]))
      }
      if(grepl(x = x$model$model.type, pattern = "Trajectory")) {
        x$points <- cbind(x$trajectories, points.dt)
      }
      else {
        x$points <- cbind(all_meta_cols, points.dt)
      }
      x$points <- as.ena.matrix(x$points, "ena.points")
      for(i in which(!find_dimension_cols(x$points))) {
        set(x$points, j = i, value = rENA::as.ena.metadata(x$points[[i]]))
      }
    }
    else {
      stop(paste0("There is no rotation matrix, if you supplied a custom ",
        "rotation.set, be sure it contains a rotation.matrix"))
    }

  ##### Calculate node positions -----
    #  - The supplied methoed is responsible is expected to return a list
    #    with two keys, "node.positions" and "centroids"
    if (exists("rotation") && !is.null(rotation) && is.null(rotation.set)) {
      optimization_rows <- x$model$points.for.projection$ENA_DIRECTION %in% optimize_on
      positions <- node.position.method(x, filter_rows = optimization_rows);

      if (all(names(positions) %in% c("node.positions", "centroids"))) {
        x$rotation$nodes <- as.data.table(positions$node.positions)
        colnames(x$rotation$nodes) <- colnames(points)
        rownames(x$rotation$nodes) <- x$rotation$codes

        for (i in seq(ncol(x$rotation$nodes))) {
          set(x$rotation$nodes, j = i, value = rENA:::as.ena.dimension(x$rotation$nodes[[i]]))
        }
        x$rotation$nodes <- data.table(
          code = structure(x$rotation$codes, class = c("code", class(x$rotation$codes))),
          x$rotation$nodes
        )
        class(x$rotation$nodes) = c("ena.nodes", class(x$rotation$nodes))

        x$model$centroids <- as.data.table(positions$centroids)
        for (i in seq(ncol(x$model$centroids))) {
          set(x$model$centroids, j = i,
            value = rENA:::as.ena.dimension(x$model$centroids[[i]])
          )
        }
        colnames(x$model$centroids) <- colnames(as.matrix(x$rotation.matrix))
        x$model$centroids = cbind(
          data.table(unit = x$model$unit.labels),
          x$model$centroids
        )
        # set(x$model$centroids, j = 1L,
        #   value = as.ena.metadata(x$model$centroids[[1L]])
        # )
        x$model$centroids <- as.ena.matrix(x$model$centroids)
      }
      else {
        stop(paste0("The node position method didn't return back the ",
          "expected objects:\n",
          "\tExpected: c('node.positions','centroids')\n",
          "\tReceived: ", names(positions), sep = ""))
      }
    }
    else if (!is.null(rotation.set)) {
      x$rotation$nodes <- rotation.set$nodes
    }

    if (is.null(x$rotation$nodes)) {
      stop("Unable to determine the node positions either by calculating
                  them using `node.position.method` or using a supplied
                  `rotation.set`")
    }

  # Class setting
    class(x$rotation) <- c("ena.rotation.set", class(x$rotation))

  # Variance ----
    var_rot_data <- var(points)
    diagonal_variance <- as.vector(diag(var_rot_data))
    x$model$variance <- diagonal_variance / sum(diagonal_variance)
    names(x$model$variance) <- colnames(x$rotation$rotation.matrix)[-1]

  # Plot object ----
    x$plots <- list() #default = ena.plot(enadata, ...))
    # class(enadata$model$plot) <- c("ena.plot", class(enadata$model$plot))

  # Additional parameters stored ----
    x$`_function.params`$norm.by <- norm.by

  return(x)
}

#' Significance Test for DirectedENA Points
#'
#' @param x directedENA model
#' @param wh logical vector of length nrow(x$points), TRUE in group1, FALSE in group2
#' @param dim TBD
#'
#' @return list of t.test results
#' @export
significance.test <- function(x, wh, dim = 1) {
  grp1_grnd <- as.matrix(x$points[ wh & ENA_DIRECTION == c("ground"),])[, dim]
  grp1_resp <- as.matrix(x$points[ wh & ENA_DIRECTION == c("response"),])[, dim]
  grp2_grnd <- as.matrix(x$points[ !wh & ENA_DIRECTION == c("ground"),])[, dim]
  grp2_resp <- as.matrix(x$points[ !wh & ENA_DIRECTION == c("response"),])[, dim]

  list(
    ground   = stats::t.test(grp1_grnd, grp2_grnd),
    response =  stats::t.test(grp1_resp, grp2_resp),
    ground_response_group1 =  stats::t.test(grp1_grnd, grp1_resp),
    ground_response_group2 =  stats::t.test(grp2_grnd, grp2_resp),
    ground_response_diff =  stats::t.test(grp1_grnd - grp1_resp, grp2_grnd - grp2_resp)
  )
}

#' Calculate node positions for a directed ena model
#'
#' @param set ena model
#' @param filter_rows logical indicating rows from model points to use
#'
#' @return list with matrix of node.positions and centroids
#'
#' @export
directed_node_optimization <- function(set, filter_rows = rep(TRUE, nrow(set$model$points.for.projection))) {
  # stop("Node optimization isn't implemented yet.")

  points = as.matrix(set$points)[filter_rows, ];
  weights = as.matrix(set$line.weights)[filter_rows, ];
  # positions = directed_node_positions(weights, points, ncol(points));
  # positions = directed_node_positions_with_ground_response_added(weights, points, ncol(points));
  if (length(unique(set$points[filter_rows,]$ENA_DIRECTION)) == 2) {
    positions = directed_node_positions_with_ground_response_added(weights, points, ncol(points));
  }
  else {
    positions = directed_node_positions(weights, points, ncol(points));
  }

  node.positions = positions$nodes;
  rownames(node.positions) = set$enadata$codes;

  return(list("node.positions" = node.positions, "centroids" = positions$centroids))
}

// [[Rcpp::depends(RcppArmadillo)]]

#include <RcppArmadillo.h>

using namespace Rcpp;

//' Upper Triangle from Vector
//'
//' @title vector to upper triangle
//' @description TBD
//' @param v TBD
//' @export
// [[Rcpp::export]]
arma::rowvec vector_to_ut(arma::mat v) {
  int vL = v.size();
  int vS = ( (vL * (vL + 1)) / 2) - vL;

  arma::rowvec vR2( vS, arma::fill::zeros );
  int s = 0;
  for( int i = 2; i <= vL; i++ ) {
    for (int j = 0; j < i-1; j++ ) {
      vR2[s] = v[j] * v[i-1];
      s++;
    }
  }
  return vR2;
}

//' Acumulate Stanza Window
//'
//' @title accum_stanza_window
//' @name accum_stanza_window
//' @description TBD
//' @param df A dataframe
//' @param windowSize Integer for number of rows in the stanza window
//' @param binary Logical, treat codes as binary or leave as weighted
//' @export
// [[Rcpp::export]]
DataFrame accum_stanza_window(
  DataFrame df,
  float windowSize = 1,
  bool binary = true
) {
  int dfRows = df.nrows();
  int dfCols = df.size();
  int numCoOccurences = dfCols * dfCols;

  arma::mat df_CoOccurred(dfRows, numCoOccurences, arma::fill::zeros);
  arma::mat df_AsMatrix(dfRows, dfCols, arma::fill::zeros);

  for (int i=0; i<dfCols;i++) {
    df_AsMatrix.col(i) = Rcpp::as<arma::vec>(df[i]);
  }

  for(int row = 0; row < dfRows; row++) {
    int earliestRow = 0, lastRow = row;

    if (windowSize == std::numeric_limits<double>::infinity()) {
      earliestRow = 0;
    }
    else if ( windowSize == 0 ) {
      earliestRow = row;
    }
    else if ( row - (windowSize-1) >= 0 ) {
      earliestRow = row - (windowSize - 1);
    }

    arma::mat currRows2 = df_AsMatrix( arma::span(earliestRow, lastRow), arma::span::all );
    arma::mat currRowsSummed = arma::sum(currRows2);
    arma::mat r = df_AsMatrix( lastRow, arma::span::all );
    arma::mat currRowAdj (dfCols, dfCols, arma::fill::zeros);
    currRowAdj.diag() = r;

    arma::mat W (dfCols, dfCols, arma::fill::zeros);


    W = (currRowsSummed.t() * r) - currRowAdj;
    // Rcpp::Rcout << "r" << std::endl;
    // Rcpp::Rcout << r << std::endl;
    // Rcpp::Rcout << "currRows2" << std::endl;
    // Rcpp::Rcout << currRows2 << std::endl;
    // Rcpp::Rcout << "currRowsSummed" << std::endl;
    // Rcpp::Rcout << currRowsSummed << std::endl;
    // Rcpp::Rcout << "currRowAdj" << std::endl;
    // Rcpp::Rcout << currRowAdj << std::endl;
    // Rcpp::Rcout << "W" << std::endl;
    // Rcpp::Rcout << W << std::endl;
    df_CoOccurred.row(row) = W.as_col().t();
  }
  if(binary == true) {
    df_CoOccurred.elem( find(df_CoOccurred > 0) ).ones();
  }

  return wrap(df_CoOccurred);
}

//' Multiobjective, Component by Component, with Ellipsoidal Scaling, for directed ENA
//'
//' @title Multiobjective, Component by Component, with Ellipsoidal Scaling, for directed ENA
//' @description TBD
//' @param line_weights TBD
//' @param points TBD
//' @param numDims TBD
//' @export
// [[Rcpp::export]]
Rcpp::List directed_node_positions(arma::mat line_weights, arma::mat points, int numDims) { //, bool by_column = true) { // = R_NilValue ) {
  int numNodes = ceil(std::sqrt(static_cast<double>(line_weights.n_cols)));

  arma::mat node_weights = arma::mat(line_weights.n_rows, numNodes, arma::fill::zeros); // zc: added an extra column

  int row_count = line_weights.n_rows;
  for (int k = 0; k < row_count; k++) {
    arma::mat currAdj = line_weights.row(k);

    int z = 0;
    for(int x = 0; x < numNodes; x++) {
      for(int y = 0; y < numNodes; y++) {
        node_weights(k,x) = node_weights(k,x) + currAdj(z);
        z = z + 1;
      }
    }
  }

  for (int k = 0; k < row_count; k++) {
    double length = 0;
    for(int i = 0; i < numNodes; i++) {
      length = length + std::abs(node_weights(k,i));
    }
    if(length < 0.0001) {
      length = 0.0001;
    }
    for(int i = 0; i < numNodes; i++) {
      node_weights(k,i) = node_weights(k,i) / length;
    }
  }

  arma::mat ssX = arma::mat(numDims, numNodes, arma::fill::zeros);
  arma::mat ssA = node_weights.t() * node_weights;
  arma::mat ssb;
  for(int i = 0; i < numDims; i++) {
    ssb = node_weights.t() * points.col(i);
    ssX.row(i) = arma::solve(ssA, ssb, arma::solve_opts::equilibrate  ).t();
  }

  arma::mat centroids = (ssX * node_weights.t()).t();

  return Rcpp::List::create(
    _("nodes") = ssX.t(),
    //_("correlations") = compute_difference_correlations(centroids, t),
    _("centroids") = centroids,
    _("weights") = node_weights, // zc: remember that the last column is all 1
    _("points") = points
  );
}

//' Node position optimization with ground and response weights/points added
//'
//' @title Node position optimization with ground and response weights/points added
//' @description TBD
//' @param line_weights TBD
//' @param points TBD
//' @param numDims TBD
//' @export
// [[Rcpp::export]]
Rcpp::List directed_node_positions_with_ground_response_added(arma::mat line_weights, arma::mat points, int numDims) { //, bool by_column = true) { // = R_NilValue ) {
  int numNodes = ceil(std::sqrt(static_cast<double>(line_weights.n_cols)));

  arma::mat node_weights = arma::mat(line_weights.n_rows, numNodes, arma::fill::zeros);

  int row_count = line_weights.n_rows;
  for (int k = 0; k < row_count; k++) {
    arma::mat currAdj = line_weights.row(k);

    int z = 0;
    for(int x = 0; x < numNodes; x++) {
      for(int y = 0; y < numNodes; y++) {
        node_weights(k,x) = node_weights(k,x) + currAdj(z);
        z = z + 1;
      }
    }
  }

  for (int k = 0; k < row_count; k++) {
    double length = 0;
    for(int i = 0; i < numNodes; i++) {
      length = length + std::abs(node_weights(k,i));
    }
    if(length < 0.0001) {
      length = 0.0001;
    }
    for(int i = 0; i < numNodes; i++) {
      node_weights(k,i) = node_weights(k,i) / length;
    }
  }
  // the following block is to add ground and response node weights/points
  arma::mat node_weights_added = arma::mat(line_weights.n_rows/2, numNodes, arma::fill::zeros);
  arma::mat points_added = arma::mat(line_weights.n_rows/2, numDims, arma::fill::zeros);

  for(int k=0;k<row_count;k+=2)
  {
    for(int i=0;i<numNodes;i++)
      node_weights_added(k/2,i)=node_weights(k,i)+node_weights(k+1,i);
    for(int i=0;i<numDims;i++)
      points_added(k/2,i)=points(k,i)+points(k+1,i);
  }
  arma::mat ssX = arma::mat(numDims, numNodes, arma::fill::zeros);
  arma::mat ssA = node_weights_added.t() * node_weights_added;
  arma::mat ssb;
  for(int i = 0; i < numDims; i++) {
    ssb = node_weights_added.t() * points_added.col(i);
    ssX.row(i) = arma::solve(ssA, ssb, arma::solve_opts::equilibrate  ).t();
  }

  arma::mat centroids = (ssX * node_weights.t()).t();

  return Rcpp::List::create(
    _("nodes") = ssX.t(),
    //_("correlations") = compute_difference_correlations(centroids, t),
    _("centroids") = centroids,
    _("weights") = node_weights,
    _("points") = points
  );
}

// [[Rcpp::export]]
Rcpp::NumericMatrix center_data_c(arma::mat values) {
  arma::mat centered = values.each_row() - mean(values);
  // arma::mat centered2(values.n_rows, values.n_cols);
  // arma::colvec m = mean(values);
  // for(int i = 0; i < values.n_rows; i++) {
  //   centered2.row(i) = values.row(i) - m;
  // }
  return Rcpp::wrap(centered);
}

/*** R
# dat <- data.table::data.table(
#   Name=c("J","Z"),
#   Day=c(1,1,1,1,1,1,2,2,2,2,2,2),
#   c1=c(1,1,1,1,1,0,0,1,1,0,0,1),
#   c2=c(1,1,1,0,0,1,0,1,0,1,0,0),
#   c3=c(0,0,1,0,1,0,1,0,0,0,1,0),
#   c4=c(1,1,1,0,0,1,0,1,0,1,0,0)
# );
# dat_acc <- dat[, {
#     accum_stanza_window(.SD, windowSize = 2, binary = TRUE)
#   },
#   by = c("Day"),
#   .SDcols = c("c1", "c2", "c3", "c4")
# ]
#
# lines2 <- matrix(c(0,1,0,1,1,0,1,0,0,0,0,1), ncol = 3, dimnames = list(NULL, LETTERS[1:3]))
# accum_stanza_window(lines2, windowSize = 4, binary = FALSE)

# dat <- read.csv(system.file("extdata", "devils_advocate_data.csv", package = "directedENA"))
# code_cols <- colnames(dat)[8:13];
#
# set_resp_rot <- directed_ena(dat,
#   units = c("Devils.Advocate", "Group", "Speaker"),
#   conversations = c("Group", "Round", "Time"),
#   codes = code_cols,
#   binary = TRUE,
#   windowSize = 4,
#   rotation_on = "response"
# )
#
# points = as.matrix(set_resp_rot$points) #[filter_rows, ];
# weights = as.matrix(set_resp_rot$line.weights) #[filter_rows, ];
# positions = directed_node_positions(weights, points, ncol(points))
# positions2 = directed_node_positions_2(weights, points, ncol(points))
#
# rows <- set_resp_rot$points$ENA_DIRECTION == "ground"
# correlations(pts = points[rows, ], cts = positions$centroids[rows, ], direction = "ground")
# correlations(pts = points[rows, ], cts = positions2$centroids[rows, ], direction = "ground")
#
# rows <- set_resp_rot$points$ENA_DIRECTION == "response"
# correlations(pts = points[rows, ], cts = positions$centroids[rows, ], direction = "response")
# correlations(pts = points[rows, ], cts = positions2$centroids[rows, ], direction = "response")
#
# set_grnd_rot <- directed_ena(dat,
#   units = c("Devils.Advocate", "Group", "Speaker"),
#   conversations = c("Group", "Round", "Time"),
#   codes = code_cols,
#   binary = TRUE,
#   windowSize = 4,
#   rotation_on = "ground"
# )
# # filter_rows <- set$line.weights$ENA_DIRECTION == "response"
# points2 = as.matrix(set_grnd_rot$points) #[filter_rows, ];
# weights2 = as.matrix(set_grnd_rot$line.weights) #[filter_rows, ];
# positions21 = directed_node_positions(weights2, points2, ncol(points2))
# positions22 = directed_node_positions_2(weights2, points2, ncol(points2))
#
# rows2 <- set_grnd_rot$points$ENA_DIRECTION == "ground"
# correlations(pts = points2[rows2, ], cts = positions2$centroids[rows2, ], direction = "ground")
#
# rows2 <- set_grnd_rot$points$ENA_DIRECTION == "response"
# correlations(pts = points[rows2, ], cts = positions2$centroids[rows2, ], direction = "response")
#
# microbenchmark::microbenchmark(
#     positions1 = directed_node_positions(weights, points, ncol(points)),
#     positions2 = directed_node_positions_2(weights2, points2, ncol(points2))
# )

print("Done")
*/

1		#' Accumulation for Directed Networks
2		#'
3		#' @param x object to accumulate
4		#' @param units character vector of columns to use as units
5		#' @param conversations character vector of columns to use as conversations
6		#' @param codes character vector of columns to use as codes
7		#' @param binary logical TRUE (default) will converat row accumulations to binary
8		#' @param weight_using TBD
9		#' @param windowSize numeric integer representing window size
10		#' @param ... TBD
11		#'
12		#' @return ena model object
13		#' @export
14		directed_accumulation <- function(
15		x,
16		units, conversations, codes,
17		binary = TRUE, weight_using = NULL,
18		windowSize = 4, ...
19		) {
20		#### Prepare input -----
21	27x	dat <- NULL
22	27x	set <- NULL
23
24	27x	if(inherits(x = x, what = "data.frame")) {
25	27x	set <- ena.set.directed(x, units, conversations, codes);
26		}
27	!	else if(inherits(x = x, what = "ena.set")) {
28	!	set <- x;
29		}
30		else {
31	!	stop("`x` must be an existing ena.set.directed or a data.frame or data.table in which to create a set from")
32		}
33
34	27x	dat <- data.table::copy(set$model$raw.input);
35
36		#### Calcuate row.connection.counts -----
37	27x	all.cols <- c(units, codes)
38	27x	new.cols <- paste("V", seq(length(codes) ^ 2), sep = "")
39	27x	row.connection.counts <- dat[,
40	27x	(new.cols) := accum_stanza_window(df = .SD[, .SD, .SDcols = codes], windowSize = windowSize, binary = binary),
41	27x	by = conversations,
42	27x	.SDcols = all.cols
43		]
44	27x	for(i in new.cols) {
45	275x	set(row.connection.counts, j = i, value = as.ena.co.occurrence(row.connection.counts[[i]]))
46		}
47
48		#### Set the column classes -----
49	27x	for(i in colnames(row.connection.counts)) {
50	482x	if(i %in% c(units, conversations))
51	60x	set(row.connection.counts, j = i, value = rENA::as.ena.metadata(row.connection.counts[[i]]))
52	422x	else if (i %in% c(codes))
53	85x	set(row.connection.counts, j = i, value = rENA:::as.ena.code(row.connection.counts[[i]]))
54		}
55	27x	row.connection.counts <- as.ena.matrix(x = row.connection.counts, "row.connections")
56
57		#### Additional columns added -----
58	27x	row.connection.counts$ENA_UNIT <- rENA::as.ena.metadata(rENA::merge_columns_c(df = row.connection.counts, cols = units, sep = "."));
59
60		#### Calculate connection.counts (Sum rows by unit)
61	27x	connection.counts <- row.connection.counts[, lapply(.SD, sum), .SDcols = new.cols, by = c("ENA_UNIT", units)]
62	27x	for(i in which(colnames(connection.counts) %in% c("ENA_UNIT", units))) {
63	60x	set(connection.counts, j = i, value = rENA::as.ena.metadata(connection.counts[[i]]))
64		}
65	27x	connection.counts <- as.ena.matrix(x = connection.counts, "ena.connections")
66
67		##### Check weighting of accumulation ----
68	27x	if( binary == FALSE && !is.null(weight_using) ) {
69	!	browser()
70		}
71
72		#### Store values on the ena.set object -----
73	27x	set$connection.counts <- connection.counts;
74	27x	set$model$row.connection.counts <- row.connection.counts;
75
76	27x	set$meta.data <- connection.counts[,rENA::find_meta_cols(connection.counts), with = FALSE];
77	27x	set$`_function.params`$binary <- binary
78	27x	set$`_function.params`$windowSize <- windowSize
79
80		#### Done -----
81	27x	return(set)
82		}

1		.onLoad <- function(libname, pkgname) {
2	!	utils::globalVariables(c(".", "ENA_DIRECTION", "code", "col2rgb", "combn", "connections", "dist",
3	!	"distance_from_endpoint_center", "hsv", "is", "mass", "name", "node_proportions",
4	!	"outer_size", "prcomp", "proportion", "receiver", "receiving", "reduce", "relative_size",
5	!	"rgb2hsv", "saturation", "sender", "sending", "size", "var", "x", "x0", "x1", "y", "y0", "y1"))
6		}
7

1		#' Generate directed ENA network
2		#'
3		#' @param data data.frame
4		#' @param ... additional parameters. `accumulated_data` to bypass accumulation, also passed to plot
5		#' @param optimize_node_positions logical
6		#' @param plotted_points data.frame
7		#' @param node_positions data.frame
8		#' @param print_plot logical
9		#' @param only_sending_centroid logical
10		#' @param units TBD
11		#' @param codes TBD
12		#' @param conversations TBD
13		#' @param rotation_on TBD
14		#' @param optimize_on TBD
15		#'
16		#' @return object of type ena.set.directed
17		#' @export
18		#'
19		#' @examples
20		#'
21		#' data("RS.data", package = "rENA")
22		#' rs_code_cols <- colnames(RS.data)[15:20]
23		#' rs_set_one <- directed_ena(RS.data,
24		#' units = c("Condition", "GroupName", "UserName"),
25		#' conversations = c("Condition", "GroupName", "ActivityNumber"),
26		#' codes = rs_code_cols,
27		#' binary = TRUE,
28		#' rotation_on = "response",
29		#' optimize_on = "response",
30		#' windowSize = 4, print_plot = FALSE
31		#' )
32		directed_ena = function(
33		data, units, conversations, codes,
34		...,
35		optimize_node_positions = TRUE,
36		only_sending_centroid = FALSE,
37		plotted_points = NULL,
38		node_positions = NULL,
39		rotation_on = "response",
40		optimize_on = c("ground", "response"),
41		print_plot = FALSE
42		) {
43	!	args <- list(...);
44
45	!	set <- ena.set.directed(data, units, conversations, codes);
46
47	!	if(!is.null(args$accumulated_data)) {
48	!	data <- args$accumulated_data;
49	!	if(is.data.frame(data$weights[[1]])) {
50	!	adjacency_matrices <- data$weights;
51		}
52		else {
53	!	adjacency_matrices = format_adjacency_matrix(data$weights);
54		}
55	!	set$adjacency_matrix = generate_adjacency_matrix(adjacency_matrices);
56		}
57		else {
58		# This should create set$connection.counts and set$model$row.connection.counts
59	!	set <- directed_accumulation(x = set, units = units, conversations = conversations, codes = codes, ...);
60		}
61
62	!	if(is.null(plotted_points)) {
63	!	set <- directed_model(set, rotation_on = rotation_on, optimize_on = optimize_on, ...);
64		}
65		else {
66	!	set$points <- plotted_points;
67		}
68
69	!	if(is.null(node_positions)) {
70	!	if(optimize_node_positions == FALSE) {
71	!	set$rotation$nodes = generate_node_coordinates(codes, total_node_connections_calculator(abs((as.square.data.frame(colMeans(as.matrix(set$connection.counts)), names = set$`_function.params`$codes)))))[];
72		}
73		else {
74		# pass abs(adjacency_matrix) in order to account for subtraction matrix
75		# set$nodes = generate_node_coordinates(set$codes, NULL);
76
77		# adjacency_matrices <- list(as.square.data.frame(as.numeric(set$connection.counts[1,2:10]), names = set$codes));
78		# if (only_sending_centroid) {
79		# set$centroid_coordinates = calculate_centroid(adjacency_matrices);
80		# }
81		# else {
82		# # NOTE - directed network centroid is midpoint between tail and head of centroid vector
83		# set$centroid_coordinates = calculate_centroid_vectors(adjacency_matrices);
84		# }
85		# set$rotation$nodes$node_radius <- max_node_diameter(set$rotation$nodes) / 2;
86
87
88		# if(!is.null(connections_calculator)) {
89		# set$rotation$nodes[, size := sapply(code, connections_calculator)];
90		# }
91		# else {
92		# set$rotation$nodes[, size := 1];
93		# }
94
95		# max_size = max(set$rotation$nodes$size);
96
97		# dims <- set$rotation$nodes[, rENA::find_dimension_cols(set$rotation$nodes), with = F];
98		# set$rotation$nodes$x0 <- dims[, 1] - node_radius
99		# set$rotation$nodes$x1 <- dims[, 1] + node_radius
100		# set$rotation$nodes$y0 <- dims[, 2] - node_radius
101		# set$rotation$nodes$y1 <- dims[, 2] + node_radius
102		}
103		}
104		else {
105	!	set$nodes <- node_positions;
106		}
107
108	!	if(print_plot == TRUE) {
109	!	set <- plot.ena.directed.set(x = set, ...);
110	!	invisible(set);
111		}
112		else {
113	!	return(set);
114		}
115		}
116
117		#' Create directed ena set object
118		#'
119		#' @param data data.frame
120		#' @param ... Not implemented
121		#' @param units TBD
122		#' @param conversations TBD
123		#' @param codes TBD
124		#'
125		#' @return empty list of type ena.set.directed
126		#' @export
127		ena.set.directed <- function(data, units, conversations, codes, ...) {
128	27x	as.ena.directed.set(
129	27x	list(
130	27x	meta.data = NULL,
131	27x	model = list(
132	27x	model.type = "Directed",
133	27x	raw.input = data.table::as.data.table(data)
134		),
135	27x	rotation = list (
136	27x	codes = codes,
137	27x	nodes = NULL
138		),
139	27x	plots = list(),
140	27x	"_function.params" = list(
141	27x	units = units,
142	27x	conversations = conversations,
143	27x	codes = codes
144		)
145		)
146		)
147		}
148
149

1		#' Plot a directed ena model
2		#'
3		#' @param x ena.set.directed to plot
4		#' @param ... parameters to pass along - Not yet implemented
5		#' @param only_sending_centroid logical
6		#' @param self_connections_in_center logical
7		#' @param node_size_dynamic logical
8		#' @param node_color Default: "#000000"
9		#' @param print_plot logical
10		#' @param multiplier logical
11		#' @param multiplier_nodes logical
12		#' @param multiplier_edges logical
13		#' @param units List
14		#' @param dimensions numeric vector of length two
15		#' @param with_mean logical
16		#' @param with_points logical
17		#' @param with_edges logical
18		#' @param with_nodes logical
19		#' @param with_node_labels logical
20		#' @param node_center_diameter character, either "self" (default), "response", or "ground"
21		#' @param node_outer_diameter character, either "self", "response" (default), or "ground"
22		#' @param node_shape character, see: https://plotly.com/r/reference/#scatter-marker-symbol
23		#' @param scale_points logcical, default: TRUE
24		#' @param colors matrix, see directedENA:::COLORS_HSV
25		#' @param title character, default is empty
26		#' @param plot_width numeric, default is NULL allowing for dynamic sizing
27		#' @param plot_height see `plot_width`
28		#' @param edge_scale TBD
29		#'
30		#' @return ENA model with plot attached
31		#' @export
32		#'
33		#' @examples
34		#'
35		#' data("RS.data", package = "rENA")
36		#' rs_code_cols <- colnames(RS.data)[15:20]
37		#' rs_set_one <- directed_ena(RS.data,
38		#' units = c("Condition", "GroupName", "UserName"),
39		#' conversations = c("Condition", "GroupName", "ActivityNumber"),
40		#' codes = rs_code_cols,
41		#' binary = TRUE,
42		#' rotation_on = "response",
43		#' optimize_on = "response",
44		#' windowSize = 4, print_plot = FALSE
45		#' )
46		#' plot(rs_set_one)
47		plot.ena.directed.set <- function(
48		# Model options
49		x, ...,
50		units = NULL,
51		dimensions = 1:2,
52		print_plot = TRUE,
53
54		# Multipliers for scaling elements
55		multiplier = 5,
56		multiplier_nodes = multiplier,
57		multiplier_edges = multiplier,
58
59		# What to include
60		with_mean = FALSE,
61		with_points = FALSE,
62		with_edges = TRUE,
63		with_nodes = TRUE,
64		with_node_labels = TRUE,
65
66		# Configuration
67		only_sending_centroid = FALSE, # Deprecate if we remove support for custom centroid_coordinates on `x`
68		self_connections_in_center = TRUE,
69		node_center_diameter = "self",
70		node_outer_diameter = "response",
71		node_size_dynamic = TRUE,
72		node_shape = "circle",
73		node_color = "#000000",
74
75		edge_scale = c(0.5,1),
76		scale_points = TRUE,
77		colors = NULL,
78		title = NULL,
79		plot_width = NULL,
80		plot_height = plot_width
81		) {
82	!	args <- list(...);
83
84	!	if(is.null(args$print_plot)) {
85	!	args$print_plot <- print_plot;
86		}
87
88	!	if(!is.null(units)) {
89	!	if(inherits(units, "ena.points")) {
90	!	adjacency_matrix <- colMeans(as.matrix(units))
91	!	if(is.null(colors)) {
92	!	colors <- COLORS_HSV[, 1, drop = FALSE]
93		}
94		}
95	!	else if (inherits(units, "ena.line.weights")) {
96	!	adjacency_matrix <- as.square.data.frame(
97	!	matrix(colMeans(as.matrix(units)),nrow = 1),
98	!	names = x$`_function.params`$codes
99		)
100	!	if(is.null(colors)) {
101	!	colors <- COLORS_HSV[, 1, drop = FALSE]
102		}
103	!	units = list("Units" = x$line.weights$ENA_UNIT %in% units$ENA_UNIT)
104		}
105		else {
106	!	adjacency_matrix <- Reduce(`-`, lapply(units, function(u) {
107	!	as.square.data.frame(
108	!	colMeans(
109	!	as.matrix(x$line.weights[u & ENA_DIRECTION == "response"])
110		),
111	!	names = x$`_function.params`$codes
112		)
113		}))
114	!	if(is.null(colors)) {
115	!	colors <- COLORS_HSV[, seq(length(units)), drop = FALSE]
116		}
117		}
118		}
119		else {
120	!	adjacency_matrix <- as.square.data.frame(colMeans(as.matrix(x$line.weights[ENA_DIRECTION == "response",])), names = x$`_function.params`$codes);
121	!	if(is.null(colors)) {
122	!	colors <- COLORS_HSV
123		}
124	!	units <- list("Units" = rep(TRUE, nrow(x$line.weights)))
125		}
126	!	if(is.character(colors)) {
127	!	colors <- apply(col2rgb(colors), 2, rgb2hsv)
128		}
129
130	!	nodes <- data.table::copy(x$rotation$nodes);
131	!	nodes$size = 1
132	!	node_radius = max_node_diameter(nodes) / 2;
133	!	nodes[, ':=' (relative_size = size / max(size)) ];
134	!	nodes[, ':=' (node_radius = relative_size * node_radius) ];
135	!	all_dim_cols <- find_dimension_cols(nodes);
136	!	dim_cols <- which(all_dim_cols)[dimensions];
137	!	other_cols <- which(!all_dim_cols);
138	!	nodes <- nodes[, .SD, .SDcols = c(names(other_cols), names(dim_cols)), with = TRUE]
139	!	setnames(nodes, old = names(dim_cols), new = c("x", "y"))
140
141	!	midpoints = generate_midpoints(nodes, abs(adjacency_matrix))[];
142	!	directed_edges = generate_directed_edges(nodes, adjacency_matrix)[];
143	!	sending_receiving = assign_directed_edge_coordinates(directed_edges$sending_receiving, nodes, midpoints, multiplier = multiplier_edges);
144	!	max_axis_marker <- max(abs(nodes$x), abs(nodes$y)) + 0.5;
145
146		# Pulling it all together now, create shapes to plot
147		# center_node_connections <- NULL;
148	!	shapes <- NULL
149	!	if(node_size_dynamic) {
150	!	if(!is.null(self_connections_in_center)) {
151
152	!	if(node_outer_diameter %in% c("ground", "response")) {
153	!	wh_sr_column = c("receiver")
154	!	if(node_outer_diameter == "ground") {
155	!	wh_sr_column = c("sender")
156		}
157
158		# browser()
159	!	nodes[, c("size", "outer_size", "inner_size", "shape", "color", "sign") := {
160	!	outer <- sum(directed_edges$sending_receiving[
161	!	directed_edges$sending_receiving[, {
162	!	total_connections = apply(.SD, 2, function(col) {
163	!	col == code
164		})
165	!	if (length(total_connections) == 1) {
166	!	return(total_connections > 0)
167		} else{
168	!	rowSums(total_connections) > 0
169		}
170	!	}, .SDcols = c(wh_sr_column)]
171	!	, (connections)]
172		);
173	!	inner <- directed_edges$self_connections[sender == as.vector(.BY$code)][, (connections)];
174		# multi <- multiplier_nodes;
175	!	this_color <- node_color;
176	!	if(outer >= 0) {
177	!	this_sign <- "positive"
178	!	this_color <- hsv(colors[1,1],colors[2,1], colors[3,1])
179		}
180		else {
181	!	this_sign <- "negative"
182	!	this_color <- hsv(colors[1,2], colors[2,2], colors[3,2])
183		}
184	!	list(
185	!	(abs(outer) + abs(inner)),
186	!	(abs(outer)),
187	!	(abs(inner)),
188	!	node_shape,
189	!	this_color,
190	!	this_sign
191		);
192	!	}, by = c("code")]
193		}
194	!	else if (node_outer_diameter == "self") {
195	!	nodes[directed_edges$self_connections, outer_size := connections, on = c("code" = "sender")]
196		}
197
198		# if (self_connections_in_center) {
199		# center_node_connections = assign_center_node_coordinates(
200		# center_node_connections = directed_edges$self_connections
201		# ,nodes = nodes
202		# );
203		# }
204		# else {
205		# center_node_connections = assign_center_node_coordinates(
206		# center_node_connections = directed_edges$external_connections
207		# ,nodes = nodes
208		# );
209		# }
210		}
211
212		# shapes <- generate_node_shapes(center_node_connections, nodes, node_size_dynamic = node_size_dynamic, colors = colors);
213		}
214		else {
215	!	nodes[ ,c("size", "outer_size", "shape", "color") := list((size * multiplier_nodes), 0, node_shape, node_color), by = c("code")]
216		}
217
218		# shapes <- lapply(shapes_df, to_shapes_list);
219	!	nodes_dims <- as.matrix(x$rotation$nodes);
220	!	shapes <- list()
221
222		# Create Plot Object -----
223	!	plot <- plot_ly(width = plot_width, height = plot_height) %>%
224	!	layout(
225	!	autosize = TRUE
226	!	,showlegend = TRUE
227	!	,title = title
228	!	,xaxis = list(type = "linear", title = "", autorange = TRUE, range = c(-max_axis_marker, max_axis_marker), scaleanchor = "y")
229	!	,yaxis = list(type = "linear", title = "", autorange = TRUE, range = c(-max_axis_marker, max_axis_marker))
230		)
231
232		# Allow editing of the plot ----
233		# NOTE: This doesn't allow renaming/moving annotations, so may note be worth keeping
234	!	if( !is.null(args$editable) ) {
235	!	plot <- config( p = plot, editable = args$editable )
236		}
237
238		# Add node shapes ----
239		# TODO: Revisit this, it may not be needed, with the nodes being added below
240	!	if(!is.null(shapes)) {
241	!	plot <- layout( p = plot ) #, shapes = shapes )
242		}
243
244		# Add points to plot ----
245		# browser()
246	!	plot <- plot_directed_points(
247	!	plot = plot, set = x, units = units,
248	!	with_mean = with_mean, with_points = with_points,
249	!	scale_units = scale_points,
250	!	colors = colors, ...
251		);
252
253		# Add edges ----
254	!	if(with_edges == TRUE) {
255	!	edge_shapes <- create_edge_shapes(
256	!	plot = plot, sending_receiving, midpoints = midpoints,
257	!	colors = colors,
258	!	thickness_range = edge_scale,
259	!	set = x,
260	!	multiplier = multiplier_edges
261		)
262	!	shapes <- c(shapes, edge_shapes)
263		}
264
265		# Add nodes and their labels ------
266	!	if(with_nodes == TRUE) {
267	!	node_shapes <- create_node_shapes(plot = plot, nodes = nodes, colors = node_color, set = x, multiplier = multiplier_nodes)
268	!	shapes <- c(shapes, node_shapes)
269		# plot <- add_markers(
270		# p = plot
271		# ,x = ~x
272		# ,y = ~y
273		# ,data = nodes
274		# ,name = "Nodes (response)"
275		# ,legendgroup = "nodes-received"
276		# ,hovertemplate = paste('(%{x}, %{y})<extra></extra>')
277		# ,marker = list(
278		# size = ~abs(size) * multiplier_nodes
279		# ,color = node_color
280		# ,symbol = node_shape
281		# ,opacity = 0.9
282		# ,sizemode = rep("diameter", nrow(nodes))
283		# ,sizeref = rep(1, nrow(nodes))
284		# ,sizemin = rep(1, nrow(nodes))
285		# )
286		# ,text = ~code
287		# ,legendgroup = "nodes"
288		# ) %>% add_markers(
289		# ,x = ~x
290		# ,y = ~y
291		# ,data = nodes
292		# ,name = "Nodes (self)"
293		# ,legendgroup = "nodes-self"
294		# ,hovertemplate = paste('(%{x}, %{y})<extra></extra>')
295		# ,marker = list(
296		# size = ~abs(inner_size) * multiplier_nodes
297		# ,color = ~color
298		# ,symbol = node_shape
299		# ,opacity = 0.9
300		# ,sizemode = rep("diameter", nrow(nodes))
301		# ,sizeref = rep(1, nrow(nodes))
302		# ,sizemin = rep(1, nrow(nodes))
303		# )
304		# ,text = ~code
305		# ,legendgroup = "nodes"
306		# )
307		}
308	!	if(with_node_labels) {
309	!	plot <- add_text(
310	!	p = plot
311	!	,x = nodes_dims[, dimensions[1]]
312	!	,y = nodes_dims[, dimensions[2]]
313	!	,text = nodes$code
314	!	,showlegend = FALSE
315	!	,legendgroup = "nodes"
316	!	,textposition = "top right"
317		)
318		}
319
320		# Custom centroids ----
321		## TODO: Deprecate
322	!	if(!is.null(x$centroid_coordinates)) {
323	!	if (only_sending_centroid) {
324	!	plot <- plot_centroids(plot, x$centroid_coordinates);
325		}
326		else {
327	!	plot <- plot_centroid_vectors(plot, x$centroid_coordinates);
328		}
329		}
330
331		# Add shapes ----
332	!	plot <- layout(p = plot, shapes = shapes)
333
334		# Save plot to model ----
335	!	x$plots[[length(x$plots) + 1]] <- plot;
336
337		# Finish and return ----
338	!	if(args$print_plot == TRUE) {
339	!	x$`_plot_op` = TRUE
340	!	return(plot);
341		# invisible(x);
342		}
343		else {
344	!	return(x);
345		}
346
347		}

1		color_saturation = function(connections) {
2	!	normalized_range = c(min(abs(connections)), max(abs(connections)));
3	!	return(
4	!	rescale(x=abs(connections), from=normalized_range, to=c(0.2,1))
5		);
6		}
7
8		network_color = function(network_element, colors = COLORS_HSV[, 1:2]) {
9	!	if (network_element$connections >= 0) {
10	!	color = colors[,1];
11		}
12		else {
13	!	color = colors[,2];
14		}
15
16		# return(hsv(color[1], network_element$saturation, color[3]));
17	!	return(hsv(color[1], color[2] * network_element$saturation, color[3]));
18		}
19
20		generate_outer_node_shapes = function(nodes, center_node_connections) {
21	!	color = COLORS_HSV[, 1];
22	!	if (is_empty(center_node_connections)) {
23	!	outer_nodes_df = nodes[, .(code, size, x0, x1, y0, y1)][, node_proportions := 1]
24		}
25		else {
26	!	outer_nodes_df = merge(
27	!	nodes[, .(code, size)]
28	!	,center_node_connections[, .(sender, connections, x0, x1, y0, y1)]
29	!	,by.x="code"
30	!	,by.y="sender"
31	!	)[,node_proportions := (size-abs(connections))/size];
32		}
33
34		# outer_node_shapes = outer_nodes_df[, .(
35		# xref = 'x', x0, x1
36		# ,yref = 'y', y0, y1
37		# ,type = 'circle'
38		# ,fillcolor = hsv(color[1], node_proportions, color[3])
39		# ,opacity = 1
40		# ,text = size
41		# ,line = list(width = 1, color = "black")
42		# )] %>% split(nodes$code) %>% unname();
43	!	outer_node_shapes <- lapply(outer_nodes_df$code, function(cd) {
44	!	code <- outer_nodes_df[code == cd,]
45	!	list(
46	!	xref = 'x', x0 = code$x0, x1 = code$x1
47	!	,yref = 'y', y0 = code$y0, y1 = code$y1
48	!	,type = 'circle'
49	!	,fillcolor = "black" # hsv(color[1], code$node_proportions, color[3])
50	!	,opacity = 1
51	!	,text = code$size
52	!	,line = list(width = 1, color = "black")
53		)
54		})
55	!	return(outer_node_shapes);
56		}
57
58		# Default: there is no outer_node_shapes; nodes are plotted just as black dots
59		generate_node_shapes <- function(
60		center_node_connections,
61		nodes = NULL,
62		node_size_dynamic = FALSE,
63		center_diameter = "self", #c("self", "response", "ground")
64		total_diameter = "response", #c("self", "response", "ground")
65		colors = COLORS_HSV[, 1:2]
66		) {
67	!	outer_node_shapes = list();
68	!	if (is_empty(center_node_connections)) {
69	!	outer_node_shapes <- generate_outer_node_shapes(nodes, NULL)
70	!	return (outer_node_shapes)
71		}
72		else {
73	!	center_node_shapes = list();
74	!	browser()
75	!	for( i in 1:length(center_node_connections$connections) ) {
76	!	row = center_node_connections[i, ];
77	!	center_node_shapes[[i]] = data.table(
78	!	xref = 'x', x0 = row$x0, x1 = row$x1
79	!	,yref = 'y', y0 = row$y0, y1 = row$y1
80	!	,type='circle'
81	!	,fillcolor = network_color(row, colors = colors)
82	!	,text = abs(row$connections)
83	!	,hoveron = 'fills'
84	!	,hoverinfo = 'text'
85		#,line = list(color='black', width=3)
86		)
87		}
88	!	return(append(outer_node_shapes, center_node_shapes));
89		}
90		# }
91
92		}
93
94		to_shapes_list = function(shape_row) {
95	!	if (is_empty(shape_row)) { return() };
96	!	shape_list = as.list(shape_row)
97	!	shape_list[['line']] = list(width=1, color='black')
98	!	return(shape_list);
99		}
100
101		triangle_path <- function(send, recv, width, multiplier = 1, adjuster = 600) {
102	!	x0 <- recv$x;
103	!	y0 <- recv$y;
104	!	x1 <- send$x;
105	!	y1 <- send$y;
106
107	!	slope <- orthogonal_slope(x0, y0, x1, y1);
108	!	adjustment <- adjuster / multiplier;
109
110	!	sine <- sin_theta_two_nodes(send, recv);
111	!	cosine <- cos_theta_two_nodes(send, recv);
112
113	!	half_width <- width / 2;
114	!	width_adj <- half_width / adjustment;
115
116		# x0 = send$x - distance_from_endpoint_center*sin_theta_two_nodes(send, receive);
117		# y0 = send$y + distance_from_endpoint_center*cos_theta_two_nodes(send, receive);
118		# x1 = send$x + distance_from_endpoint_center*sin_theta_two_nodes(send, receive);
119		# y1 = send$y - distance_from_endpoint_center*cos_theta_two_nodes(send, receive);
120	!	path <- paste0("M ", x0, " ", y0, " ",
121	!	"L ", (x1 + (width_adj * sine)), " ", (y1 - (width_adj * cosine)), " ",
122	!	"L ", (x1 - (width_adj * sine)), " ", (y1 + (width_adj * cosine)), " ",
123	!	"Z")
124
125	!	return(path)
126		}
127
128		create_node_shapes <- function(
129		plot,
130		nodes,
131		colors = COLORS_HSV[, 1:2],
132		set = NULL,
133		multiplier = 1,
134		adjuster = 600
135		) {
136		# Shape definition ----
137	!	shape = list(
138	!	type = "circle",
139	!	xanchor = "x",
140	!	yanchor = "y",
141	!	xref = "x",
142	!	yref = "y",
143	!	fillcolor = "black",
144	!	line = list (color = "black")
145		)
146
147		# Generate shapes ----
148	!	shapes = list()
149		# browser()
150	!	adjustment <- adjuster / multiplier;
151		# nodes$size <- c(0.1066974, 0.1280369, nodes$size[3], 0)
152	!	for(i in 1:nrow(nodes)) {
153	!	node = nodes[i, ];
154
155		# node_color <- node$color;
156	!	this_shape <- shape;
157	!	half_width <- node$size / 2;
158	!	width_adj <- half_width / adjustment;
159
160	!	this_shape$x0 <- as.numeric(node$x) - width_adj
161	!	this_shape$y0 <- as.numeric(node$y) - width_adj
162	!	this_shape$x1 <- as.numeric(node$x) + width_adj
163	!	this_shape$y1 <- as.numeric(node$y) + width_adj
164	!	shapes[[length(shapes) + 1]] <- this_shape
165
166	!	if(!is.null(node$inner_size) && node$inner_size > 0) {
167	!	this_shape <- shape;
168	!	node_color <- node$color;
169	!	half_width <- node$inner_size / 2;
170	!	width_adj <- half_width / adjustment;
171	!	this_shape$fillcolor <- node_color;
172	!	this_shape$line$color <- node_color;
173	!	this_shape$x0 <- as.numeric(node$x) - width_adj
174	!	this_shape$y0 <- as.numeric(node$y) - width_adj
175	!	this_shape$x1 <- as.numeric(node$x) + width_adj
176	!	this_shape$y1 <- as.numeric(node$y) + width_adj
177	!	shapes[[length(shapes) + 1]] <- this_shape
178		}
179		}
180
181		# Return updated plot ----
182		# plot <- layout(p = plot, shapes = shapes)
183		# return(plot);
184	!	return(shapes);
185		}
186
187		create_edge_shapes = function(
188		plot,
189		edges, midpoints,
190		colors = COLORS_HSV[, 1:2],
191		thickness_range = c(ifelse(min(edges$connections) == 0, 0, 0.1), 1),
192		set = NULL,
193		multiplier = 1,
194		adjustment = 600
195		) {
196		# Setup ----
197		# thickness = c(min(abs(edges$connections)), max(abs(edges$connections)));
198	!	opacity = rescale(x=abs(edges$connections), from=c(0,1), to = thickness_range);
199
200		# message("Note 5: I think we may want the saturation to be multiplied by color[2]")
201
202		# Hidden point ----
203		### - specifically for use with the legend
204	!	plot = add_trace(
205	!	plot
206	!	,legendgroup = "network"
207	!	,type = "scatter"
208	!	,hoverinfo = NULL
209	!	,mode = "markers"
210	!	,name = "Network"
211	!	,opacity = 1
212	!	,x = 0
213	!	,y = 0
214	!	,marker = list ( color = "#000", size = 0.1 )
215		)
216
217		# Edge Shapes ----
218	!	shape = list(
219	!	type = "path",
220	!	xanchor = "x",
221	!	yanchor = "y",
222	!	xref = "x",
223	!	fillcolor = "black",
224	!	yref = "y",
225	!	line = list (color = "black")
226		)
227
228	!	shapes = list()
229	!	for(i in 1:nrow(edges)) {
230	!	row = edges[i, ];
231	!	edge_color <- network_color(row, colors = colors);
232	!	send <- set$rotation$nodes[code == row$sender]
233	!	recv <- set$rotation$nodes[code == row$receiver];
234	!	mid <- midpoints[sending %in% c(send$code, recv$code) & receiving %in% c(send$code, recv$code)]
235	!	if(nrow(mid) == 0) {
236	!	midpoint <- list(x = recv$SVD1, y = recv$SVD2);
237		} else {
238	!	midpoint <- list(x = mid$midpoint_x, y = mid$midpoint_y);
239		}
240		# print(row$connections)
241	!	this_shape <- shape;
242	!	this_shape$path <- triangle_path(
243	!	list(x=send$SVD1, y=send$SVD2),
244	!	midpoint,
245	!	width = row$connections,
246	!	multiplier = multiplier,
247	!	adjuster = adjustment
248		)
249	!	this_shape$line$color <- edge_color;
250	!	this_shape$fillcolor <- edge_color;
251	!	this_shape$opacity <- opacity[i]
252	!	this_shape$name <- abs(row$connections)
253	!	this_shape$legendgroup = "network"
254
255	!	shapes[[length(shapes) + 1]] <- this_shape
256		# plot = add_trace(
257		# plot
258		# ,fill = 'toself'
259		# ,fillcolor = edge_color
260		# ,legendgroup = "network"
261		# ,hoverinfo = 'text'
262		# ,hoveron = 'fills'
263		# ,marker = list(opacity = 0) # don't plot midpoints
264		# ,mode = "markers"
265		# ,opacity = opacity[i]
266		# ,text = abs(row$connections)
267		# ,name = NULL
268		# ,type = 'scatter'
269		# ,showlegend = FALSE
270		# ,x = as.numeric(row[, grep(x = colnames(row), pattern = "^x_path"), with = FALSE])
271		# ,y = as.numeric(row[, grep(x = colnames(row), pattern = "^y_path"), with = FALSE])
272		# )
273		}
274
275		# plot <- layout(p = plot, shapes = shapes)
276
277		# Return updated plot ----
278		# return(plot);
279	!	return(shapes);
280		}
281
282		plot_centroids = function(plot, centroid_info) {
283	!	for (i in 1:length(centroid_info$x)) {
284	!	color = centroid_info$color[[i]];
285	!	plot = add_trace(
286	!	plot
287	!	, x= centroid_info$x[i]
288	!	, y = centroid_info$y[i]
289	!	, text='centroid'
290	!	, marker = list(
291	!	color = color
292		)
293		)
294		}
295	!	return(plot);
296		}
297
298		plot_directed_points <- function(
299		plot, set, units,
300		labels = names(units), dimensions = 1:2,
301		scale_units = TRUE,
302		opacity_range = c(0.25,0.75),
303		with_mean = TRUE,
304		with_points = TRUE,
305		just_vectors = TRUE,
306		colors = COLORS_HSV
307		) {
308
309	!	for(i in seq(length(units))) {
310		# browser()
311	!	point_color_hsv <- colors[, i];
312	!	point_color <- hsv(point_color_hsv[1], point_color_hsv[2], point_color_hsv[3])
313	!	wh_units <- units[[i]]
314	!	is_mean <- length(wh_units) > 1
315	!	resp <- as.matrix(set$points[wh_units & ENA_DIRECTION == "response",])
316	!	grnd <- as.matrix(set$points[wh_units & ENA_DIRECTION == "ground",])
317
318	!	resp_pts <- resp[, dimensions, drop = FALSE]
319	!	grnd_pts <- grnd[, dimensions, drop = FALSE]
320	!	scaleFactor = 1.0
321
322	!	if( scale_units == TRUE ) {
323	!	np.min.x = abs(min(as.matrix(set$rotation$nodes)[, dimensions[1]]));
324	!	np.min.y = abs(min(as.matrix(set$rotation$nodes)[, dimensions[2]]));
325	!	rp.min.x = abs(min(as.matrix(set$points)[, dimensions[1]]));
326	!	rp.min.y = abs(min(as.matrix(set$points)[, dimensions[2]]));
327	!	maxMin = abs(max(np.min.x / rp.min.x, np.min.y / rp.min.y));
328
329	!	np.max.x = abs(max(as.matrix(set$rotation$nodes)[, dimensions[1]]));
330	!	np.max.y = abs(max(as.matrix(set$rotation$nodes)[, dimensions[2]]));
331	!	rp.max.x = abs(max(as.matrix(set$points)[, dimensions[1]]));
332	!	rp.max.y = abs(max(as.matrix(set$points)[, dimensions[2]]));
333	!	maxMax = abs(max(np.max.x / rp.max.x, np.max.y / rp.max.y));
334	!	scaleFactor = min(maxMin, maxMax);
335		}
336	!	else if (is.numeric(scale_units)) {
337	!	scaleFactor = scale_units;
338		}
339		# browser()
340	!	resp_pts <- resp_pts * scaleFactor
341	!	grnd_pts <- grnd_pts * scaleFactor
342
343	!	pts_df <- data.frame(x = c(grnd_pts[,1], resp_pts[,1]), y = c(grnd_pts[,2], resp_pts[,2]));
344	!	distances <- sapply(seq(nrow(grnd_pts)), function(r) dist(rbind(grnd_pts[r,], resp_pts[r,])));
345	!	arrow_colors <- matrix(rep(point_color_hsv, length(distances)), ncol = 3, byrow = TRUE);
346
347		# arrow_colors[, 2] <- arrow_colors[, 2] * (distances);
348		# opacity <- scales::rescale(x = distances, to = c(0, point_color_hsv[2]))
349	!	opacity <- scales::rescale(x = distances, to = c(opacity_range[1], opacity_range[2]))
350
351		# arrow_color_hex <- apply(arrow_colors, 1, function(x) hsv(x[1], x[2], x[3]));
352	!	arrow_color_hex <- sapply(seq(nrow(arrow_colors)), function(i) {
353	!	x <- arrow_colors[i, ];
354	!	hsv(x[1], x[2], x[3], opacity[i])
355		});
356	!	if(with_points == TRUE) {
357	!	if(just_vectors == FALSE) {
358	!	plot <- add_markers(
359	!	plot
360	!	,x = pts_df$x
361	!	,y = pts_df$y
362	!	,text = labels[i]
363	!	,legendgroup = "points"
364	!	,name = paste0(labels[i], " Points")
365	!	,marker = list(
366	!	color = point_color
367	!	,size = ~size * multiplier
368	!	,sizemode = rep("diameter", nrow(pts_df))
369	!	,sizeref = rep(1, nrow(pts_df))
370	!	,sizemin = rep(1, nrow(pts_df))
371		)
372		)
373		}
374	!	plot <- add_annotations(
375	!	plot
376	!	,ax = grnd_pts[,1]
377	!	,x = resp_pts[,1]
378	!	,axref = "x"
379	!	,ay = grnd_pts[,2]
380	!	,y = resp_pts[,2]
381	!	,ayref = "y"
382	!	,showarrow = TRUE
383	!	,text = ''
384	!	,xref = 'x'
385	!	,yref = 'y'
386	!	,arrowcolor = arrow_color_hex
387	!	,legendgroup = "points"
388		)
389		}
390
391	!	if(with_mean == TRUE && is_mean == TRUE && nrow(grnd_pts) > 1) {
392	!	mean_color <- hsv(point_color_hsv[1], point_color_hsv[2], point_color_hsv[3], alpha = 1.0);
393	!	plot <- add_markers(
394	!	plot
395	!	,x = c(mean(grnd_pts[, 1])) #, mean(resp_pts[, 1]))
396	!	,y = c(mean(grnd_pts[, 2])) #, mean(resp_pts[, 2]))
397		# ,x = c(mean(resp_pts[, 1])) #, mean(resp_pts[, 1]))
398		# ,y = c(mean(resp_pts[, 2])) #, mean(resp_pts[, 2]))
399	!	,text = paste0(labels[i], " Mean")
400	!	,marker = list( color = mean_color, symbol = "square" )
401	!	,legendgroup = "means"
402	!	,name = paste0(labels[i], " Mean")
403		)
404		# browser()
405	!	plot <- add_annotations(
406	!	plot
407	!	,ax = mean(grnd_pts[,1])
408	!	,ay = mean(grnd_pts[,2])
409	!	,x = mean(resp_pts[,1])
410	!	,y = mean(resp_pts[,2])
411	!	,axref = "x"
412	!	,ayref = "y"
413	!	,showarrow = TRUE
414	!	,text = ''
415		# ,xref = 'x'
416		# ,yref = 'y'
417	!	,arrowcolor = mean_color
418	!	,legendgroup = "means"
419		)
420		# plot <- add_annotations(
421		# plot
422		# ,ax = 0
423		# ,x = 1
424		# ,axref = "x"
425		# ,ay = 0
426		# ,y = 1
427		# ,ayref = "y"
428		# ,showarrow = TRUE
429		# ,text = ''
430		# # ,xref = 'x'
431		# # ,yref = 'y'
432		# ,arrowcolor = arrow_color_hex
433		# ,legendgroup = "means"
434		# )
435		}
436		}
437	!	return(plot)
438		}
439
440		plot_centroid_vectors = function(plot, vector_data) {
441	!	for (i in 1:length(vector_data)) {
442	!	data = vector_data[i] %>% unlist;
443	!	if (data['sending_x'] == data['receiving_x'] && data['sending_y'] == data['receiving_y']) {
444	!	plot = add_trace(
445	!	plot
446	!	,x = as.numeric(data['sending_x'])
447	!	,y = as.numeric(data['sending_y'])
448	!	,text = 'centroid'
449	!	,marker = list( color = data['color'])
450		)
451		}
452		else {
453	!	plot = add_annotations(
454	!	plot
455	!	,ax = as.numeric(data['sending_x'])
456	!	,x = as.numeric(data['receiving_x'])
457	!	,axref = "x"
458	!	,ay = as.numeric(data['sending_y'])
459	!	,y = as.numeric(data['receiving_y'])
460	!	,ayref = "y"
461	!	,showarrow = TRUE
462	!	,text = ''
463	!	,xref = 'x'
464	!	,yref = 'y'
465	!	,arrowcolor = data['color']
466		)
467		}
468		}
469	!	return(plot);
470		}

1		to_square <- function(x) {
2	626x	n <- sqrt(length(x))
3
4	626x	dim(x) <- c(n, n);
5	626x	x
6		}
7
8		#' Re-class vector as ena.co.occurrence
9		#'
10		#' @param x Vector to re-class
11		#'
12		#' @return re-classed vector
13		#' @export
14		as.ena.co.occurrence <- function(x) {
15	425x	if(is.factor(x)) {
16	!	x = as.character(x)
17		}
18	425x	class(x) = c("ena.co.occurrence", class(x))
19	425x	x
20		}
21
22		#' Title
23		#'
24		#' @param x directed ENA model
25		#'
26		#' @return matrix
27		#' @export
28		as.directed.matrix <- function(x) {
29	!	return(x);
30		}
31
32		#' Title
33		#'
34		#' @param x directed ENA model
35		#' @param names character
36		#'
37		#' @return square matrix
38		#' @export
39		as.square <- function(x, names = NULL) {
40	!	UseMethod("as.square", x);
41		}
42
43		#' Title
44		#'
45		#' @param x directed ENA model
46		#' @param ... unused
47		#'
48		#' @return square matrix
49		#' @export
50		as.square.ena.connections <- function(x, ...) {
51	!	x_ <- as.matrix(x);
52	!	square_matrices <- lapply(seq(nrow(x_)), function(y) {
53	!	y <- as.ena.matrix(to_square(x_[y, ]), "ena.directed.connections")
54		})
55	!	square_matrices
56		}
57
58		#' Title
59		#'
60		#' @param x directed ENA model
61		#' @param names character
62		#'
63		#' @return square matrix
64		#' @export
65		as.square.matrix <- function(x, names = NULL) {
66	!	x_ <- as.data.frame(as.matrix(x))
67	!	as.square(x_)
68		}
69
70		#' Title
71		#'
72		#' @param x vector
73		#' @param names character vector
74		#'
75		#' @return square data.frame
76		#' @export
77		as.square.data.frame <- function(x, names = NULL) {
78	!	x_ <- as.matrix(x)
79	!	n <- sqrt(length(x_))
80
81	!	dim(x_) <- c(n, n);
82	!	df <- as.data.frame(x_)
83
84	!	if(is.null(names)) {
85	!	names <- paste("V", seq(n), sep = "")
86		}
87	!	rownames(df) <- colnames(df) <- names;
88	!	df <- as.ena.matrix(df, "ena.directed.connections")
89	!	df
90		}
91
92		#' Convert object to a directed set
93		#'
94		#' @param x Object to convert
95		#'
96		#' @return list as ena.directed.set
97		#' @export
98		as.ena.directed.set <- function(x) {
99	27x	to <- c("ena.directed.set")
100	27x	if(!inherits(x = x, what = "ena.set")) {
101	27x	to <- c(to, "ena.set")
102		}
103	27x	class(x) <- c(to, class(x));
104	27x	return(x);
105		}
106
107		#' Print a directed ena model
108		#'
109		#' @param x model to print
110		#' @param ... not implemented
111		#' @param plot logical if TRUE, print the plots (default is FALSE)
112		#'
113		#' @export
114		print.ena.set.directed <- function(x, ..., plot = FALSE) {
115	!	x.unclass <- unclass(x)
116
117		if(
118	!	!is.null(x.unclass$`_plot_op`) &&
119	!	x.unclass$`_plot_op` == T
120		) {
121	!	base::print(x.unclass$plots)
122		}
123		else {
124	!	if(plot == FALSE) {
125	!	x.unclass$plots <- NULL
126		}
127	!	base::print(x.unclass)
128		}
129		}
130
131
132		is_empty <- function(x) {
133	!	return(is.null(x) \|\| nrow(x) == 0 \|\| ncol(x) == 0);
134		}
135
136		to_datatable = function(appended_vectors, ncol, nrow) {
137		# because we're representing the data as a vector of row vectors, we need to
138		# reverse the matrix ncol/nrow and transpose the matrix to get the correct output
139	!	return(data.table(t(matrix(appended_vectors, nrow = ncol, ncol = nrow))));
140		}
141
142		two_point_weighted_average = function(point_1, point_2, relative_weight) {
143		# distance formula from:
144		# https://www.khanacademy.org/partner-content/pixar/environment-modeling-2/mathematics-of-parabolas2-ver2/v/weighted-average-two-points
145	!	return((1-relative_weight)point_1 + (relative_weightpoint_2));
146		}
147
148		orthogonal_slope = function(receiver_x, receiver_y, sender_x, sender_y) {
149	!	sender_to_receiver_slope = (sender_y - receiver_y)/(sender_x - receiver_x);
150	!	return(-(1/sender_to_receiver_slope));
151		}
152
153		cos_theta_two_nodes = function (sender, receiver) {
154	!	hypotenus = sqrt((receiver$x - sender$x)^2 + (receiver$y - sender$y)^2);
155	!	adjacent = receiver$x - sender$x;
156	!	return(adjacent/hypotenus);
157		}
158
159		sin_theta_two_nodes = function (sender, receiver) {
160	!	hypotenus = sqrt((receiver$x - sender$x)^2 + (receiver$y - sender$y)^2);
161	!	opposite = receiver$y - sender$y;
162	!	return(opposite/hypotenus);
163		}
164
165		is_vertical_line = function(slope) {
166	!	return(abs(slope) == Inf);
167		}
168
169		is_horizontal_line = function(slope) {
170	!	return(is.nan(slope));
171		}
172
173		format_adjacency_matrix = function(raw_data) {
174		# node_names = raw_data$node_names;
175		# nrows = length(node_names);
176		# return(function(raw_matrix) {
177		# matrix = as.data.frame(to_datatable(raw_matrix, nrows, nrows));
178		# colnames(matrix) = rownames(matrix) = node_names;
179		# return(matrix);
180		# })
181	!	node_names <- raw_data$node_names;
182	!	node_len <- length(node_names);
183	!	lapply(raw_data$raw_data, function(x) {
184	!	as.data.frame(matrix(x, ncol = node_len, nrow = node_len, dimnames = list(node_names, node_names), byrow = TRUE))
185		})
186		}
187
188		raw_coordinates_to_datatable = function(names) {
189	!	nrows = length(names);
190	!	pts.circle <- t(sapply(1:nrows, function(r) c(cos(2rpi/nrows), sin(2rpi/nrows))))
191	!	colnames(pts.circle) = c('x', 'y');
192
193		# ncols = 3; # 3 columns, nodename, x, and y
194		# nodes = to_datatable(raw_data$coordinates, ncols, nrows);
195	!	nodes <- data.table("name" = names);
196	!	nodes <- cbind(nodes, pts.circle);
197
198		# nodes[, ':=' (x = as.double(x), y = as.double(y))];
199	!	return(nodes);
200		}
201
202		# -------------------------------------
203		# DATA TRANSFORMATIONS
204		# -------------------------------------
205
206		generate_adjacency_matrix = function(matrices) {
207		# guard clauses
208	!	if (length(matrices) > 2) stop("cannot pass this function more than 2 adjacency matrices");
209	!	if (typeof(matrices) != 'list') stop('raw data for adjacency matrix must be in a list')
210
211		# business logic
212	!	if ( length(matrices) == 2 ) {
213		# message("Note 4: need to remove reduce, use Reduce?")
214	!	adjacency_matrix = reduce(matrices, `-`);
215		}
216		else {
217	!	adjacency_matrix = matrices[[1]];
218		}
219
220	!	return(adjacency_matrix);
221		}
222
223		total_node_connections_calculator = function(adjacency_matrix) {
224		# return "partially applied" anonymous function in order to generate column values
225	!	return(
226	!	function(node_name) {
227		# total node connections should just be sending connections
228	!	total = rowSums(adjacency_matrix[node_name,]);
229		# HACK: make sure node is not invisible, set minimum to 1;
230	!	if (total <= 0) { total = 1 };
231	!	return(total);
232		}
233		)
234		}
235
236		max_node_diameter = function(nodes){
237	!	pts <- as.matrix(nodes)
238	!	plot_area = (max(pts[,1]) - min(pts[,1])) * (max(pts[,2]) - min(pts[,2]));
239		# assume that multiplying by the max scalar gives us the area of a square, which should encompass
240		# the outer bounds of the largest plotted node.
241		# take the square root to get the base of the square, which we will use as the max
242		# diameter size, translated to our plot's cartesian plane.
243		# all other nodes and edges will be pegged to this max node radius
244	!	return(sqrt(MAX_NODE_DIAMETER_SCALAR * plot_area));
245		}
246
247		generate_node_coordinates = function(codes, connections_calculator = NULL) {
248	!	nodes = raw_coordinates_to_datatable(codes);
249	!	node_radius = max_node_diameter(nodes) / 2;
250
251	!	if(!is.null(connections_calculator)) {
252	!	nodes[, size := sapply(name, connections_calculator)];
253		}
254		else {
255	!	nodes[, size := 1];
256		}
257
258	!	max_size = max(nodes$size);
259	!	nodes[, ':=' (relative_size = size / max_size) ];
260	!	nodes[, ':=' (node_radius = relative_size * node_radius) ];
261
262	!	nodes[, ':=' (
263	!	x0 = x - node_radius
264	!	,x1 = x + node_radius
265	!	,y0 = y - node_radius
266	!	,y1 = y + node_radius
267		)];
268	!	return(nodes);
269		};
270
271		# NOTE: technically, these labels are somewhat misleading since this matrix represents
272		# the total edge weights of both sender and receiver. However, in order to identify
273		# the weighted midpoint between two nodes, we arbitrarily assign one node to be a
274		# sender and another to be a receiver.
275		# This allows us to calculate the proportion of edge weight from the sender to receiver out
276		# of the total connections between the two nodes.
277		generate_midpoints = function(nodes, adjacency_matrix) {
278	!	midpoints = data.table((t(combn(nodes$code, 2))), stringsAsFactors = FALSE);
279	!	colnames(midpoints) = c('sending', 'receiving');
280
281		# browser()
282	!	midpoints[, c("sender_receiver_connections","total_connections", "proportion", "midpoint_x", "midpoint_y") := {
283	!	dim_cols <- which(find_dimension_cols(nodes))
284	!	send <- adjacency_matrix[sending, receiving];
285	!	recv <- adjacency_matrix[receiving, sending];
286	!	prop <- send / (send + recv)
287	!	pts_s <- as.numeric(nodes[code == sending, dim_cols, with = FALSE])
288	!	pts_r <- as.numeric(nodes[code == receiving, dim_cols, with = FALSE])
289
290	!	mid_x <- two_point_weighted_average( point_1 = pts_s[1], point_2 = pts_r[1], relative_weight = prop);
291	!	mid_y <- two_point_weighted_average( point_1 = pts_s[2], point_2 = pts_r[2], relative_weight = prop);
292
293	!	list(send, send + recv, prop, mid_x, mid_y);
294	!	}, by = c("sending", "receiving")];
295
296		# this will also remove rows where proportion is NaN
297	!	return(midpoints[proportion != 0]);
298		}
299
300		generate_directed_edges = function(nodes, adjacency_matrix) {
301		# directed_edges = data.table(
302		# permutations(n=length(nodes$name), r=2, v=nodes$name, repeats.allowed = TRUE)
303		# )
304	!	cols_sr <- c("sender", "receiver");
305	!	directed_edges <- data.table(
306	!	expand.grid("receiver"=nodes$code[order(nodes$code)],
307	!	"sender"=nodes$code[order(nodes$code)], stringsAsFactors = FALSE)
308	!	)[, cols_sr, with = FALSE];
309	!	directed_edges[, connections := adjacency_matrix[sender, receiver], by = cols_sr];
310	!	directed_edges[, saturation := color_saturation(connections)]
311
312	!	sending_receiving = copy(directed_edges[sender != receiver & abs(connections) > 0]);
313	!	self_connections = copy(directed_edges[sender == receiver]);
314	!	external_connections = copy(directed_edges[sender == receiver]);
315
316	!	if (!is_empty(external_connections)) {
317		# message("Note 1: `external_connections` wasn't saving the mutation, seems questionable. Also, the
318		# resulting connections were all the same.")
319		# external_connections %>% mutate(connections = nodes[name == sender]$size - connections);
320	!	external_connections[, connections := {
321	!	nodes[code == sender]$size - connections
322		},
323	!	by = sender
324		]
325		}
326	!	edges_list <- list(
327	!	sending_receiving = sending_receiving
328	!	,self_connections = self_connections
329	!	,external_connections = external_connections
330		);
331
332	!	edges = lapply(edges_list, function(edge_df) {
333	!	edge_df[, relative_size := { abs(connections) / nodes[code == sender]$size }, by = sender];
334	!	edge_df[, distance_from_endpoint_center := { relative_size * nodes[code == sender]$node_radius }, by = sender];
335
336	!	return(copy(edge_df));
337		});
338
339	!	return(edges);
340		}
341
342		assign_center_node_coordinates = function(center_node_connections, nodes) {
343	!	if(is_empty(center_node_connections)) return(center_node_connections);
344
345	!	if("distance_from_endpoint_center" %in% colnames(center_node_connections)) {
346	!	setnames(center_node_connections, old = "distance_from_endpoint_center", new = "radius")
347		}
348
349		# xy_cols <- colnames(as.matrix(nodes))[1:2];
350	!	xy_cols <- names(which(find_dimension_cols(nodes)))
351
352	!	center_node_connections <- center_node_connections[nodes,
353	!	c(colnames(center_node_connections), xy_cols),
354	!	on = c(sender = "code"),
355	!	with = FALSE][order(nodes$code)]
356		# browser()
357	!	dims <- center_node_connections[, rENA::find_dimension_cols(center_node_connections), with = F];
358	!	center_node_connections$x0 <- dims[, 1] - nodes$node_radius
359	!	center_node_connections$x1 <- dims[, 1] + nodes$node_radius
360	!	center_node_connections$y0 <- dims[, 2] - nodes$node_radius
361	!	center_node_connections$y1 <- dims[, 2] + nodes$node_radius
362
363	!	center_node_connections
364		};
365
366		calculate_centroid_vectors = function(matrices) {
367		# matrices = format_adjacency_matrix(raw_data);
368
369	!	nodes = raw_coordinates_to_datatable(colnames(matrices[[1]]));
370	!	mass_vectors = list();
371	!	for (i in 1:length(matrices)) {
372		# From above, this appears to assume each item in matrices is already a matrix, but it's still a function
373	!	adjacency_matrix = matrices[[i]];
374	!	total_mass = sum(adjacency_matrix);
375	!	node_sending = rowSums(adjacency_matrix);
376	!	all_sending = data.table(name = names(node_sending), mass = as.vector(node_sending));
377	!	sending_mass = merge(all_sending, nodes[, .(name,x,y)], by="name")[, ':=' (
378	!	mass_x = mass*x
379	!	, mass_y = mass*y
380		)];
381	!	node_receiving = colSums(adjacency_matrix);
382	!	all_receiving = data.table(name = names(node_receiving), mass = as.vector(node_receiving));
383	!	receiving_mass = merge(all_receiving, nodes[, .(name,x,y)], by="name")[, ':=' (
384	!	mass_x = mass*x
385	!	, mass_y = mass*y
386		)];
387	!	color = COLORS_HSV[, i+1]
388	!	sending_x = sum(sending_mass$mass_x) / total_mass;
389	!	sending_y = sum(sending_mass$mass_y) / total_mass;
390	!	receiving_x = sum(receiving_mass$mass_x) / total_mass;
391	!	receiving_y = sum(receiving_mass$mass_y) / total_mass;
392
393	!	mass_vectors[[i]] = list(
394	!	sending_x = sending_x
395	!	, sending_y = sending_y
396	!	, receiving_x = receiving_x
397	!	, receiving_y = receiving_y
398	!	, color = hsv(color['h'], 0.93, color['v'])
399		)
400		}
401	!	return(mass_vectors);
402		}
403
404		calculate_centroid = function(raw_data) {
405		# matrices = lapply(data$raw_data, format_adjacency_matrix(data));
406	!	matrices = format_adjacency_matrix(raw_data);
407
408	!	nodes = raw_coordinates_to_datatable(raw_data);
409	!	x_coords = c();
410	!	y_coords = c();
411	!	color = list();
412	!	for (i in 1:length(matrices)) {
413	!	centroid_color = COLORS_HSV[, i+1]
414	!	adjacency_matrix = matrices[[i]];
415	!	node_sending = rowSums(adjacency_matrix);
416	!	total_mass = sum(node_sending);
417	!	all_sending = data.table(name = names(node_sending), mass = as.vector(node_sending));
418	!	node_mass = merge(all_sending, nodes[, .(name,x,y)], by="name")[, ':=' (
419	!	mass_x = mass*x
420	!	, mass_y = mass*y
421		)];
422	!	x_coords[i] = sum(node_mass$mass_x) / total_mass;
423	!	y_coords[i] = sum(node_mass$mass_y) / total_mass;
424		# #F01186
425	!	color[[i]] = hsv(centroid_color['h'], 0.93, centroid_color['v'])
426		}
427	!	return(list(x = x_coords , y = y_coords, color = color));
428		}
429
430		assign_directed_edge_coordinates = function(sending_receiving, nodes, midpoints, multiplier = 1) {
431		# message("Note 2: Making a copy to remove DT warning for now") # I can't find the automatic R copy (clm)
432	!	sending_receiving <- copy(sending_receiving);
433		# nodes <- copy(nodes);
434		# midpoints <- copy(midpoints);
435	!	sr_cols <- c("sender", "receiver")
436	!	colvector <- c("endpoint_coordinates_x", "endpoint_coordinates_y")
437
438	!	sending_receiving[, (colvector) := {
439	!	midpoint_present_for_nodes <- midpoints[
440	!	(midpoints$sending %in% .SD$sender & midpoints$receiving %in% .SD$receiver) \|
441	!	(midpoints$sending %in% .SD$receiver & midpoints$receiving %in% .SD$sender)
442		];
443	!	if (is_empty(midpoint_present_for_nodes)) {
444	!	list(x = nodes[code == receiver]$x, y = nodes[code == receiver]$y);
445		}
446		else {
447	!	list(x = midpoint_present_for_nodes$midpoint_x, y = midpoint_present_for_nodes$midpoint_y);
448		}
449	!	}, by = sr_cols, .SDcols = sr_cols]
450
451	!	colvector_tri <- c("triangle_base_coordinates_x0", "triangle_base_coordinates_y0",
452	!	"triangle_base_coordinates_x1", "triangle_base_coordinates_y1")
453	!	sending_receiving[, (colvector_tri) := {
454	!	send = nodes[code == .SD$sender];
455	!	receive = nodes[code == .SD$receiver];
456	!	slope = orthogonal_slope(receive$x, receive$y, send$x, send$y);
457	!	sine = sin_theta_two_nodes(send, receive);
458	!	cosine = cos_theta_two_nodes(send, receive);
459	!	distance_from_endpoint_center = distance_from_endpoint_center * multiplier;
460	!	x0 = send$x - distance_from_endpoint_center*sin_theta_two_nodes(send, receive);
461	!	y0 = send$y + distance_from_endpoint_center*cos_theta_two_nodes(send, receive);
462	!	x1 = send$x + distance_from_endpoint_center*sin_theta_two_nodes(send, receive);
463	!	y1 = send$y - distance_from_endpoint_center*cos_theta_two_nodes(send, receive);
464	!	if (is_horizontal_line(slope)) {
465	!	y0 = y1 = send$y;
466	!	x0 = send$x + distance_from_endpoint_center;
467	!	x1 = send$x - distance_from_endpoint_center;
468		}
469	!	else if (is_vertical_line(slope)) {
470	!	x0 = x1 = send$x;
471	!	y0 = send$y + distance_from_endpoint_center;
472	!	y1 = send$y - distance_from_endpoint_center;
473		}
474	!	list(x0 = x0, y0 = y0, x1 = x1, y1 = y1)
475	!	}, by = sr_cols, .SDcols = sr_cols]
476
477	!	colvector_paths <- c("x_path_ep_1", "x_path_tr_1", "x_path_tr_2", "x_path_ep_2",
478	!	"y_path_ep_1", "y_path_tr_1", "y_path_tr_2", "y_path_ep_2")
479	!	sending_receiving[, (colvector_paths) := {
480	!	list(
481	!	.SD$endpoint_coordinates_x
482	!	,.SD$triangle_base_coordinates_x1
483	!	,.SD$triangle_base_coordinates_x0
484	!	,.SD$endpoint_coordinates_x
485	!	,.SD$endpoint_coordinates_y
486	!	,.SD$triangle_base_coordinates_y1
487	!	,.SD$triangle_base_coordinates_y0
488	!	,.SD$endpoint_coordinates_y
489		)
490	!	}, by = sr_cols]
491	!	return(sending_receiving[]);
492		}
493

1		#' Normalize and optimize a directed accumulation
2		#'
3		#' @param x a ena.set.directed with a set of accumulated adjacency matrices
4		#' @param norm.by TBD
5		#' @param rotate.using TBD
6		#' @param rotation.params TBD
7		#' @param rotation.set TBD
8		#' @param rotation_on TBD
9		#' @param optimize_on TBD
10		#' @param node.position.method TBD
11		#' @param ... TBD
12		#'
13		#' @return ena.set.directed with calcualted points and node positiosn
14		#' @export
15		directed_model <- function(
16		x, ...,
17		norm.by = rENA::fun_sphere_norm,
18		node.position.method = directed_node_optimization,
19		rotate.using = rENA::ena.svd,
20		rotation.params = NULL, rotation.set = NULL,
21		rotation_on = "response", #c("ground", "response"),
22		optimize_on = c("ground", "response")
23		) {
24
25		##### Generate data.frame with a ground and response row for each unit -----
26		##### > Prep ----
27	2x	meta_cols <- which(rENA::find_meta_cols(x$connection.counts))
28
29		##### > Generate ----
30	2x	df <- data.table::rbindlist(lapply(seq(nrow(x$connection.counts)), function(r) {
31	626x	row <- x$connection.counts[r, ];
32	626x	meta <- row[, meta_cols, with = FALSE]
33	626x	mat <- to_square(as.matrix(row))
34
35	626x	ground <- as.vector(t(mat))
36	626x	response <- as.vector(mat)
37
38	626x	df <- data.frame(matrix(c(ground, response), nrow = 2, byrow = TRUE, dimnames = list(NULL, colnames(row)[find_code_cols(row)])))
39	626x	df <- cbind(meta, ENA_DIRECTION = c("ground", "response"), df)
40	626x	df
41		}))
42		# setorder(x = df, cols = "ENA_DIRECTION")
43
44		##### > Set column classes ----
45	2x	class(df) <- class(x$connection.counts);
46	2x	for(i in colnames(df)) {
47	62x	if(i %in% c(names(meta_cols), "ENA_DIRECTION"))
48	12x	set(df, j = i, value = rENA::as.ena.metadata(df[[i]]))
49		else
50	50x	set(df, j = i, value = as.ena.co.occurrence(df[[i]]))
51		}
52
53		##### Model preperation -----
54	2x	all_meta_cols <- df[, find_meta_cols(df), with = FALSE]
55	2x	code_columns <- colnames(df)[rENA::find_code_cols(df)];
56
57		##### Normalize the ground/response adjacency vectors ----
58		### > Normalize adjacency vectors ----
59	2x	line.weights <- norm.by(as.matrix(df));
60
61		### > Classify objects ----
62	2x	colnames(line.weights) <- code_columns
63
64	2x	line.weights.dt <- as.data.table(line.weights)
65	2x	for (i in seq(ncol(line.weights.dt))) {
66	50x	set(line.weights.dt, j = i, value = as.ena.co.occurrence(line.weights.dt[[i]]))
67		}
68
69	2x	x$line.weights <- cbind(all_meta_cols, line.weights.dt)
70	2x	for(i in which(!find_code_cols(x$line.weights))) {
71	12x	set(x$line.weights, j = i, value = rENA::as.ena.metadata(x$line.weights[[i]]))
72		}
73	2x	class(x$line.weights) <- c("ena.line.weights", class(x$line.weights))
74
75
76		##### Center the Normalized data ----
77		### > Center the line.weights ----
78	2x	points.for.projection <- center_data_c(line.weights)
79
80		### > Classify objects ----
81	2x	colnames(points.for.projection) <- code_columns;
82	2x	x$model$points.for.projection = as.data.table(points.for.projection)
83
84	2x	for (i in seq(ncol(x$model$points.for.projection))) {
85	50x	set(x$model$points.for.projection, j = i, value = as.ena.co.occurrence(x$model$points.for.projection[[i]]))
86		}
87	2x	x$model$points.for.projection <- as.ena.matrix(cbind( all_meta_cols, x$model$points.for.projection), "ena.points")
88
89	2x	for(i in which(!find_code_cols(x$model$points.for.projection))) {
90	12x	set(x$model$points.for.projection, j = i, value = rENA::as.ena.metadata(x$model$points.for.projection[[i]]))
91		}
92
93		##### SVD -----
94
95		### > New rotation ----
96	2x	if (!is.null(rotate.using) && is.null(rotation.set)) {
97		### >> Row preparation ----
98	2x	rotation_rows <- x$model$points.for.projection$ENA_DIRECTION %in% rotation_on
99		# resp_rows <- !grnd_rows
100		# rotation <- do.call(rotate.using, list(x, rotation.params))
101		# rotation_ground <- prcomp(points.for.projection[grnd_rows,], retx=FALSE, scale=FALSE, center=FALSE, tol=0)
102		# rotation_response <- prcomp(points.for.projection[resp_rows,], retx=FALSE, scale=FALSE, center=FALSE, tol=0)
103
104	2x	pts = as.matrix(x$model$points.for.projection)[rotation_rows,]
105
106		### >> Perform SVD on normed, centered data ----
107	2x	pcaResults = prcomp(pts, retx=FALSE, scale=FALSE, center=FALSE, tol=0)
108
109		### >> Classify Results ----
110	2x	colnames(pcaResults$rotation) = c(
111	2x	paste('SVD',as.character(1:ncol(pcaResults$rotation)), sep='')
112		);
113
114		# rotationSet = ENARotationSet$new(rotation = pcaResults$pca, codes = enaset$codes, node.positions = NULL, eigenvalues = pcaResults$latent)
115	2x	rotation = ENARotationSet$new(
116	2x	rotation = pcaResults$rotation,
117	2x	codes = x$codes,
118	2x	node.positions = NULL,
119	2x	eigenvalues = pcaResults$sdev^2
120		)
121
122	2x	x$rotation.matrix <- as.data.table(rotation$rotation, keep.rownames = "codes")
123	2x	for (i in seq(ncol(x$rotation.matrix))) {
124	52x	if(i == 1) {
125	2x	set(x$rotation.matrix, j = i, value = as.ena.metadata(x$rotation.matrix[[i]]))
126		} else {
127	50x	set(x$rotation.matrix, j = i, value = rENA:::as.ena.dimension(x$rotation.matrix[[i]]))
128		}
129		}
130	2x	class(x$rotation.matrix) <- c("ena.rotation.matrix", class(x$rotation.matrix))
131
132	2x	x$rotation$rotation.matrix <- x$rotation.matrix
133	2x	x$rotation$eigenvalues <- rotation$eigenvalues;
134		}
135		### > Custom rotation provided ----
136	!	else if (!is.null(rotation.set)) {
137	!	if (is(rotation.set, "ena.rotation.set")) {
138	!	x$rotation.matrix <- rotation.set$rotation.matrix
139	!	x$rotation$rotation.matrix <- rotation.set$rotation.matrix
140	!	x$rotation$nodes <- rotation.set$nodes;
141	!	x$rotation$eigenvalues <- rotation.set$eigenvalues
142		}
143		else {
144	!	stop("Supplied rotation.set is not an instance of ENARotationSet")
145		}
146		}
147		else {
148	!	stop("Unable to find or create a rotation set")
149		}
150
151		##### Generate the rotated points ----
152	2x	if (!is.null(x$rotation.matrix)) {
153	2x	points <- points.for.projection %*% as.matrix(x$rotation.matrix)
154	2x	points.dt <- as.data.table(points)
155	2x	for (i in seq(ncol(points.dt))) {
156	50x	set(points.dt, j = i, value = rENA:::as.ena.dimension(points.dt[[i]]))
157		}
158	2x	if(grepl(x = x$model$model.type, pattern = "Trajectory")) {
159	!	x$points <- cbind(x$trajectories, points.dt)
160		}
161		else {
162	2x	x$points <- cbind(all_meta_cols, points.dt)
163		}
164	2x	x$points <- as.ena.matrix(x$points, "ena.points")
165	2x	for(i in which(!find_dimension_cols(x$points))) {
166	12x	set(x$points, j = i, value = rENA::as.ena.metadata(x$points[[i]]))
167		}
168		}
169		else {
170	!	stop(paste0("There is no rotation matrix, if you supplied a custom ",
171	!	"rotation.set, be sure it contains a rotation.matrix"))
172		}
173
174		##### Calculate node positions -----
175		# - The supplied methoed is responsible is expected to return a list
176		# with two keys, "node.positions" and "centroids"
177	2x	if (exists("rotation") && !is.null(rotation) && is.null(rotation.set)) {
178	2x	optimization_rows <- x$model$points.for.projection$ENA_DIRECTION %in% optimize_on
179	2x	positions <- node.position.method(x, filter_rows = optimization_rows);
180
181	2x	if (all(names(positions) %in% c("node.positions", "centroids"))) {
182	2x	x$rotation$nodes <- as.data.table(positions$node.positions)
183	2x	colnames(x$rotation$nodes) <- colnames(points)
184	2x	rownames(x$rotation$nodes) <- x$rotation$codes
185
186	2x	for (i in seq(ncol(x$rotation$nodes))) {
187	50x	set(x$rotation$nodes, j = i, value = rENA:::as.ena.dimension(x$rotation$nodes[[i]]))
188		}
189	2x	x$rotation$nodes <- data.table(
190	2x	code = structure(x$rotation$codes, class = c("code", class(x$rotation$codes))),
191	2x	x$rotation$nodes
192		)
193	2x	class(x$rotation$nodes) = c("ena.nodes", class(x$rotation$nodes))
194
195	2x	x$model$centroids <- as.data.table(positions$centroids)
196	2x	for (i in seq(ncol(x$model$centroids))) {
197	50x	set(x$model$centroids, j = i,
198	50x	value = rENA:::as.ena.dimension(x$model$centroids[[i]])
199		)
200		}
201	2x	colnames(x$model$centroids) <- colnames(as.matrix(x$rotation.matrix))
202	2x	x$model$centroids = cbind(
203	2x	data.table(unit = x$model$unit.labels),
204	2x	x$model$centroids
205		)
206		# set(x$model$centroids, j = 1L,
207		# value = as.ena.metadata(x$model$centroids[[1L]])
208		# )
209	2x	x$model$centroids <- as.ena.matrix(x$model$centroids)
210		}
211		else {
212	!	stop(paste0("The node position method didn't return back the ",
213	!	"expected objects:\n",
214	!	"\tExpected: c('node.positions','centroids')\n",
215	!	"\tReceived: ", names(positions), sep = ""))
216		}
217		}
218	!	else if (!is.null(rotation.set)) {
219	!	x$rotation$nodes <- rotation.set$nodes
220		}
221
222	2x	if (is.null(x$rotation$nodes)) {
223	!	stop("Unable to determine the node positions either by calculating
224	!	them using `node.position.method` or using a supplied
225	!	`rotation.set`")
226		}
227
228		# Class setting
229	2x	class(x$rotation) <- c("ena.rotation.set", class(x$rotation))
230
231		# Variance ----
232	2x	var_rot_data <- var(points)
233	2x	diagonal_variance <- as.vector(diag(var_rot_data))
234	2x	x$model$variance <- diagonal_variance / sum(diagonal_variance)
235	2x	names(x$model$variance) <- colnames(x$rotation$rotation.matrix)[-1]
236
237		# Plot object ----
238	2x	x$plots <- list() #default = ena.plot(enadata, ...))
239		# class(enadata$model$plot) <- c("ena.plot", class(enadata$model$plot))
240
241		# Additional parameters stored ----
242	2x	x$`_function.params`$norm.by <- norm.by
243
244	2x	return(x)
245		}

1		#' Significance Test for DirectedENA Points
2		#'
3		#' @param x directedENA model
4		#' @param wh logical vector of length nrow(x$points), TRUE in group1, FALSE in group2
5		#' @param dim TBD
6		#'
7		#' @return list of t.test results
8		#' @export
9		significance.test <- function(x, wh, dim = 1) {
10	!	grp1_grnd <- as.matrix(x$points[ wh & ENA_DIRECTION == c("ground"),])[, dim]
11	!	grp1_resp <- as.matrix(x$points[ wh & ENA_DIRECTION == c("response"),])[, dim]
12	!	grp2_grnd <- as.matrix(x$points[ !wh & ENA_DIRECTION == c("ground"),])[, dim]
13	!	grp2_resp <- as.matrix(x$points[ !wh & ENA_DIRECTION == c("response"),])[, dim]
14
15	!	list(
16	!	ground = stats::t.test(grp1_grnd, grp2_grnd),
17	!	response = stats::t.test(grp1_resp, grp2_resp),
18	!	ground_response_group1 = stats::t.test(grp1_grnd, grp1_resp),
19	!	ground_response_group2 = stats::t.test(grp2_grnd, grp2_resp),
20	!	ground_response_diff = stats::t.test(grp1_grnd - grp1_resp, grp2_grnd - grp2_resp)
21		)
22		}

1		#' Calculate node positions for a directed ena model
2		#'
3		#' @param set ena model
4		#' @param filter_rows logical indicating rows from model points to use
5		#'
6		#' @return list with matrix of node.positions and centroids
7		#'
8		#' @export
9		directed_node_optimization <- function(set, filter_rows = rep(TRUE, nrow(set$model$points.for.projection))) {
10		# stop("Node optimization isn't implemented yet.")
11
12	2x	points = as.matrix(set$points)[filter_rows, ];
13	2x	weights = as.matrix(set$line.weights)[filter_rows, ];
14		# positions = directed_node_positions(weights, points, ncol(points));
15		# positions = directed_node_positions_with_ground_response_added(weights, points, ncol(points));
16	2x	if (length(unique(set$points[filter_rows,]$ENA_DIRECTION)) == 2) {
17	1x	positions = directed_node_positions_with_ground_response_added(weights, points, ncol(points));
18		}
19		else {
20	1x	positions = directed_node_positions(weights, points, ncol(points));
21		}
22
23	2x	node.positions = positions$nodes;
24	2x	rownames(node.positions) = set$enadata$codes;
25
26	2x	return(list("node.positions" = node.positions, "centroids" = positions$centroids))
27		}

1		// [[Rcpp::depends(RcppArmadillo)]]
2
3		#include <RcppArmadillo.h>
4
5		using namespace Rcpp;
6
7		//' Upper Triangle from Vector
8		//'
9		//' @title vector to upper triangle
10		//' @description TBD
11		//' @param v TBD
12		//' @export
13		// [[Rcpp::export]]
14	!	arma::rowvec vector_to_ut(arma::mat v) {
15	!	int vL = v.size();
16	!	int vS = ( (vL * (vL + 1)) / 2) - vL;
17
18	!	arma::rowvec vR2( vS, arma::fill::zeros );
19	!	int s = 0;
20	!	for( int i = 2; i <= vL; i++ ) {
21	!	for (int j = 0; j < i-1; j++ ) {
22	!	vR2[s] = v[j] * v[i-1];
23	!	s++;
24		}
25		}
26	!	return vR2;
27		}
28
29		//' Acumulate Stanza Window
30		//'
31		//' @title accum_stanza_window
32		//' @name accum_stanza_window
33		//' @description TBD
34		//' @param df A dataframe
35		//' @param windowSize Integer for number of rows in the stanza window
36		//' @param binary Logical, treat codes as binary or leave as weighted
37		//' @export
38		// [[Rcpp::export]]
39	153x	DataFrame accum_stanza_window(
40		DataFrame df,
41		float windowSize = 1,
42		bool binary = true
43		) {
44	153x	int dfRows = df.nrows();
45	153x	int dfCols = df.size();
46	153x	int numCoOccurences = dfCols * dfCols;
47
48	306x	arma::mat df_CoOccurred(dfRows, numCoOccurences, arma::fill::zeros);
49	306x	arma::mat df_AsMatrix(dfRows, dfCols, arma::fill::zeros);
50
51	864x	for (int i=0; i<dfCols;i++) {
52	1422x	df_AsMatrix.col(i) = Rcpp::as<arma::vec>(df[i]);
53		}
54
55	24355x	for(int row = 0; row < dfRows; row++) {
56	24202x	int earliestRow = 0, lastRow = row;
57
58	24202x	if (windowSize == std::numeric_limits<double>::infinity()) {
59	!	earliestRow = 0;
60		}
61	24202x	else if ( windowSize == 0 ) {
62	!	earliestRow = row;
63		}
64	24202x	else if ( row - (windowSize-1) >= 0 ) {
65	23747x	earliestRow = row - (windowSize - 1);
66		}
67
68	72606x	arma::mat currRows2 = df_AsMatrix( arma::span(earliestRow, lastRow), arma::span::all );
69	48404x	arma::mat currRowsSummed = arma::sum(currRows2);
70	48404x	arma::mat r = df_AsMatrix( lastRow, arma::span::all );
71	48404x	arma::mat currRowAdj (dfCols, dfCols, arma::fill::zeros);
72	24202x	currRowAdj.diag() = r;
73
74	24202x	arma::mat W (dfCols, dfCols, arma::fill::zeros);
75
76
77	48404x	W = (currRowsSummed.t() * r) - currRowAdj;
78		// Rcpp::Rcout << "r" << std::endl;
79		// Rcpp::Rcout << r << std::endl;
80		// Rcpp::Rcout << "currRows2" << std::endl;
81		// Rcpp::Rcout << currRows2 << std::endl;
82		// Rcpp::Rcout << "currRowsSummed" << std::endl;
83		// Rcpp::Rcout << currRowsSummed << std::endl;
84		// Rcpp::Rcout << "currRowAdj" << std::endl;
85		// Rcpp::Rcout << currRowAdj << std::endl;
86		// Rcpp::Rcout << "W" << std::endl;
87		// Rcpp::Rcout << W << std::endl;
88	48404x	df_CoOccurred.row(row) = W.as_col().t();
89		}
90	153x	if(binary == true) {
91	258x	df_CoOccurred.elem( find(df_CoOccurred > 0) ).ones();
92		}
93
94	306x	return wrap(df_CoOccurred);
95		}
96
97		//' Multiobjective, Component by Component, with Ellipsoidal Scaling, for directed ENA
98		//'
99		//' @title Multiobjective, Component by Component, with Ellipsoidal Scaling, for directed ENA
100		//' @description TBD
101		//' @param line_weights TBD
102		//' @param points TBD
103		//' @param numDims TBD
104		//' @export
105		// [[Rcpp::export]]
106	1x	Rcpp::List directed_node_positions(arma::mat line_weights, arma::mat points, int numDims) { //, bool by_column = true) { // = R_NilValue ) {
107	1x	int numNodes = ceil(std::sqrt(static_cast<double>(line_weights.n_cols)));
108
109	2x	arma::mat node_weights = arma::mat(line_weights.n_rows, numNodes, arma::fill::zeros); // zc: added an extra column
110
111	1x	int row_count = line_weights.n_rows;
112	314x	for (int k = 0; k < row_count; k++) {
113	939x	arma::mat currAdj = line_weights.row(k);
114
115	313x	int z = 0;
116	1878x	for(int x = 0; x < numNodes; x++) {
117	9390x	for(int y = 0; y < numNodes; y++) {
118	23475x	node_weights(k,x) = node_weights(k,x) + currAdj(z);
119	7825x	z = z + 1;
120		}
121		}
122		}
123
124	314x	for (int k = 0; k < row_count; k++) {
125	313x	double length = 0;
126	1878x	for(int i = 0; i < numNodes; i++) {
127	3130x	length = length + std::abs(node_weights(k,i));
128		}
129	313x	if(length < 0.0001) {
130	74x	length = 0.0001;
131		}
132	1878x	for(int i = 0; i < numNodes; i++) {
133	4695x	node_weights(k,i) = node_weights(k,i) / length;
134		}
135		}
136
137	2x	arma::mat ssX = arma::mat(numDims, numNodes, arma::fill::zeros);
138	3x	arma::mat ssA = node_weights.t() * node_weights;
139	2x	arma::mat ssb;
140	26x	for(int i = 0; i < numDims; i++) {
141	75x	ssb = node_weights.t() * points.col(i);
142	75x	ssX.row(i) = arma::solve(ssA, ssb, arma::solve_opts::equilibrate ).t();
143		}
144
145	2x	arma::mat centroids = (ssX * node_weights.t()).t();
146
147		return Rcpp::List::create(
148	2x	_("nodes") = ssX.t(),
149		//_("correlations") = compute_difference_correlations(centroids, t),
150	2x	_("centroids") = centroids,
151	2x	_("weights") = node_weights, // zc: remember that the last column is all 1
152	2x	_("points") = points
153		);
154		}
155
156		//' Node position optimization with ground and response weights/points added
157		//'
158		//' @title Node position optimization with ground and response weights/points added
159		//' @description TBD
160		//' @param line_weights TBD
161		//' @param points TBD
162		//' @param numDims TBD
163		//' @export
164		// [[Rcpp::export]]
165	1x	Rcpp::List directed_node_positions_with_ground_response_added(arma::mat line_weights, arma::mat points, int numDims) { //, bool by_column = true) { // = R_NilValue ) {
166	1x	int numNodes = ceil(std::sqrt(static_cast<double>(line_weights.n_cols)));
167
168	2x	arma::mat node_weights = arma::mat(line_weights.n_rows, numNodes, arma::fill::zeros);
169
170	1x	int row_count = line_weights.n_rows;
171	627x	for (int k = 0; k < row_count; k++) {
172	1878x	arma::mat currAdj = line_weights.row(k);
173
174	626x	int z = 0;
175	3756x	for(int x = 0; x < numNodes; x++) {
176	18780x	for(int y = 0; y < numNodes; y++) {
177	46950x	node_weights(k,x) = node_weights(k,x) + currAdj(z);
178	15650x	z = z + 1;
179		}
180		}
181		}
182
183	627x	for (int k = 0; k < row_count; k++) {
184	626x	double length = 0;
185	3756x	for(int i = 0; i < numNodes; i++) {
186	6260x	length = length + std::abs(node_weights(k,i));
187		}
188	626x	if(length < 0.0001) {
189	148x	length = 0.0001;
190		}
191	3756x	for(int i = 0; i < numNodes; i++) {
192	9390x	node_weights(k,i) = node_weights(k,i) / length;
193		}
194		}
195		// the following block is to add ground and response node weights/points
196	2x	arma::mat node_weights_added = arma::mat(line_weights.n_rows/2, numNodes, arma::fill::zeros);
197	2x	arma::mat points_added = arma::mat(line_weights.n_rows/2, numDims, arma::fill::zeros);
198
199	314x	for(int k=0;k<row_count;k+=2)
200		{
201	1878x	for(int i=0;i<numNodes;i++)
202	6260x	node_weights_added(k/2,i)=node_weights(k,i)+node_weights(k+1,i);
203	8138x	for(int i=0;i<numDims;i++)
204	31300x	points_added(k/2,i)=points(k,i)+points(k+1,i);
205		}
206	2x	arma::mat ssX = arma::mat(numDims, numNodes, arma::fill::zeros);
207	3x	arma::mat ssA = node_weights_added.t() * node_weights_added;
208	2x	arma::mat ssb;
209	26x	for(int i = 0; i < numDims; i++) {
210	75x	ssb = node_weights_added.t() * points_added.col(i);
211	75x	ssX.row(i) = arma::solve(ssA, ssb, arma::solve_opts::equilibrate ).t();
212		}
213
214	2x	arma::mat centroids = (ssX * node_weights.t()).t();
215
216		return Rcpp::List::create(
217	2x	_("nodes") = ssX.t(),
218		//_("correlations") = compute_difference_correlations(centroids, t),
219	2x	_("centroids") = centroids,
220	2x	_("weights") = node_weights,
221	2x	_("points") = points
222		);
223		}
224
225		// [[Rcpp::export]]
226	2x	Rcpp::NumericMatrix center_data_c(arma::mat values) {
227	4x	arma::mat centered = values.each_row() - mean(values);
228		// arma::mat centered2(values.n_rows, values.n_cols);
229		// arma::colvec m = mean(values);
230		// for(int i = 0; i < values.n_rows; i++) {
231		// centered2.row(i) = values.row(i) - m;
232		// }
233	4x	return Rcpp::wrap(centered);
234		}
235
236		/*** R
237		# dat <- data.table::data.table(
238		# Name=c("J","Z"),
239		# Day=c(1,1,1,1,1,1,2,2,2,2,2,2),
240		# c1=c(1,1,1,1,1,0,0,1,1,0,0,1),
241		# c2=c(1,1,1,0,0,1,0,1,0,1,0,0),
242		# c3=c(0,0,1,0,1,0,1,0,0,0,1,0),
243		# c4=c(1,1,1,0,0,1,0,1,0,1,0,0)
244		# );
245		# dat_acc <- dat[, {
246		# accum_stanza_window(.SD, windowSize = 2, binary = TRUE)
247		# },
248		# by = c("Day"),
249		# .SDcols = c("c1", "c2", "c3", "c4")
250		# ]
251		#
252		# lines2 <- matrix(c(0,1,0,1,1,0,1,0,0,0,0,1), ncol = 3, dimnames = list(NULL, LETTERS[1:3]))
253		# accum_stanza_window(lines2, windowSize = 4, binary = FALSE)
254
255		# dat <- read.csv(system.file("extdata", "devils_advocate_data.csv", package = "directedENA"))
256		# code_cols <- colnames(dat)[8:13];
257		#
258		# set_resp_rot <- directed_ena(dat,
259		# units = c("Devils.Advocate", "Group", "Speaker"),
260		# conversations = c("Group", "Round", "Time"),
261		# codes = code_cols,
262		# binary = TRUE,
263		# windowSize = 4,
264		# rotation_on = "response"
265		# )
266		#
267		# points = as.matrix(set_resp_rot$points) #[filter_rows, ];
268		# weights = as.matrix(set_resp_rot$line.weights) #[filter_rows, ];
269		# positions = directed_node_positions(weights, points, ncol(points))
270		# positions2 = directed_node_positions_2(weights, points, ncol(points))
271		#
272		# rows <- set_resp_rot$points$ENA_DIRECTION == "ground"
273		# correlations(pts = points[rows, ], cts = positions$centroids[rows, ], direction = "ground")
274		# correlations(pts = points[rows, ], cts = positions2$centroids[rows, ], direction = "ground")
275		#
276		# rows <- set_resp_rot$points$ENA_DIRECTION == "response"
277		# correlations(pts = points[rows, ], cts = positions$centroids[rows, ], direction = "response")
278		# correlations(pts = points[rows, ], cts = positions2$centroids[rows, ], direction = "response")
279		#
280		# set_grnd_rot <- directed_ena(dat,
281		# units = c("Devils.Advocate", "Group", "Speaker"),
282		# conversations = c("Group", "Round", "Time"),
283		# codes = code_cols,
284		# binary = TRUE,
285		# windowSize = 4,
286		# rotation_on = "ground"
287		# )
288		# # filter_rows <- set$line.weights$ENA_DIRECTION == "response"
289		# points2 = as.matrix(set_grnd_rot$points) #[filter_rows, ];
290		# weights2 = as.matrix(set_grnd_rot$line.weights) #[filter_rows, ];
291		# positions21 = directed_node_positions(weights2, points2, ncol(points2))
292		# positions22 = directed_node_positions_2(weights2, points2, ncol(points2))
293		#
294		# rows2 <- set_grnd_rot$points$ENA_DIRECTION == "ground"
295		# correlations(pts = points2[rows2, ], cts = positions2$centroids[rows2, ], direction = "ground")
296		#
297		# rows2 <- set_grnd_rot$points$ENA_DIRECTION == "response"
298		# correlations(pts = points[rows2, ], cts = positions2$centroids[rows2, ], direction = "response")
299		#
300		# microbenchmark::microbenchmark(
301		# positions1 = directed_node_positions(weights, points, ncol(points)),
302		# positions2 = directed_node_positions_2(weights2, points2, ncol(points2))
303		# )
304
305		print("Done")
306		*/