ona coverage - 1.76%

Files
Source

alpha_to_hex <- function(alpha) {
  toupper((sprintf("%02s", as.hexmode(round(alpha * 255)))))
}
lighten_color <- function(col, amt = 1, alpha = 1) {
  rgb_white <- grDevices::col2rgb("#FFFFFF");

  col_as_rgb <- if(is.character(col))
      grDevices::col2rgb(col)
    else
      col_as_rgb <- col

  # rgb_white + (rgb_just_s - rgb_white) * (this_color_hsv$alpha * edge_arrow_saturation_multiplier)
  rgb_white + (col_as_rgb - rgb_white) * (alpha * amt)
}
darken_color <- function(col, amt = 1, alpha = 1) {
  rgb_white <- grDevices::col2rgb("#FFFFFF");

  col_as_rgb <- if(is.character(col))
      grDevices::col2rgb(col)
    else
      col_as_rgb <- col

  (rgb_white + (col_as_rgb - rgb_white) * (alpha)) * (1 - -(1 - amt))
}
CreatePointByDistance <- function(P1, P2, dist = 0.1) {
  P_dist <- LearnGeom::DistancePoints(P1, P2)
  if(P_dist == 0) {
    return(P1)
  }
  else {
    d_ratio <- (P_dist - dist) / P_dist
    xt <- (1 - d_ratio) * P1[1] + (d_ratio * P2[1])
    yt <- (1 - d_ratio) * P1[2] + (d_ratio * P2[2])

    P3 <- c(xt, yt)
    names(P3) <- c("X", "Y")

    return(P3)
  }
}
MoveLineAlong <- function(x1, y1, m) {
  x <- 0
  y <- m * x - m * x1 + y1

  Line <- c(m, y)
  names(Line) = c("slope", "intercept")
  class(Line) <- append(class(Line), "Line")
  Line
}
CalculateSegmentHeight <- function(radius, chord_length) {
  (radius) - ( 0.5 * sqrt(4*(radius ^ 2) - (chord_length ^ 2)) )
}
CalculateSegmentArea <- function(radius, height) {
  ((radius ^ 2) * acos((radius - height) / radius)) -
  ((radius - height) * sqrt((2 * radius * height) - (height ^ 2)))
}
edge_fun <- function(x, fun, y = NULL) {
  if(!is.matrix(x)) {
    if (!is.vector(x)) {
      stop("x is required as matrix or vector. If it's a vector, y is required to be the same.")
    }
    else if (is.null(y)) {
      stop("If x is a vector, y is required to be a vector.")
    }
  }
  else {
    y <- x[2,]
    x <- x[1,]
  }

  do.call(fun, list(x, y))
}
edge_center <- function(x, y = NULL) {
  mid <- edge_fun(x = x, fun = LearnGeom::MidPoint, y = y);
  names(mid) <- c("x", "y")

  mid
}
node_distance <- function(x, y = NULL) {
  dist <- edge_fun(x = x, fun = LearnGeom::DistancePoints, y = y);
  names(dist) <- c("distance")
  dist
}
edge_center_weighted <- function(x, weight, toward = 1) {
  edge_line <- edge_fun(x = x, fun = LearnGeom::CreateLinePoints);
  edge_mid <- edge_center(x)
  edge_len <- node_distance(x)

  pt_wh <- ifelse(weight > 0, 1, 2)
  pt_ref <- x[toward, ]

  v_dist <- edge_len * weight;
  if(edge_line[1] == Inf) {
    int <- matrix(c(
        rep(pt_ref[1], 2),
        pt_ref[2] + v_dist,
        pt_ref[2] - v_dist
      ),
      ncol = 2)
  }
  else if (edge_line[1] == 0) {
    # stop("No edge")
    int <- matrix(c(
      pt_ref[1] + v_dist,
      pt_ref[1] - v_dist,
      rep(pt_ref[2], 2)
    ),
    ncol = 2)
  }
  else {
    if(all(abs(diff(x)) < (.Machine$double.eps ^ 0.5))) {
      int <- x;
    }
    else {
      int <- LearnGeom::IntersectLineCircle(edge_line, pt_ref, v_dist)
    }
  }
  int_dists <- apply(int, 1, node_distance, x[-toward,])

  x_sorted <- sort(x[,1])
  y_sorted <- sort(x[,2])
  on_segment <- apply(int, 1, function(p) {
    p[1] >= x_sorted[1] && p[1] <= x_sorted[2] &&
    p[2] >= y_sorted[1] && p[2] <= y_sorted[2]
  })
  if(all(on_segment)) {
    int_closest <- which(int_dists == min(int_dists))
  }
  else {
    int_closest <- which(on_segment)
  }
  w_c <- int[int_closest, ]

  w_c
}

create_edge_matrix <- function(
  set, dims = 1:2,
  weights = colMeans(as.matrix(set$line.weights)),
  nodes = set$rotation$nodes,
  direction = 1,
  node_position_multiplier = 1
) {
  raw_weights <- weights;
  weights <- abs(weights);
  weights_matrix <- to_square(weights);
  raw_weights_matrix <- to_square(raw_weights);
  rownames(weights_matrix) <- colnames(weights_matrix) <- nodes$code

  nodes[, (seq(2, ncol(nodes))) := lapply(.SD, `*`, node_position_multiplier) , .SDcols = (seq(2, ncol(nodes)))];
  nodes_xy <- nodes[, c(1, dims+1), with = F]
  colnames(nodes_xy) <- c("code", "x", "y")
  node_combinations <- combn(length(nodes_xy[[1]]), 2);

  edge_names_matrix <- t(matrix(nodes$code[node_combinations], nrow = 2))
  edge_table <- data.table::data.table(node_1 = edge_names_matrix[, 1], node_2 = edge_names_matrix[, 2])
  edge_table[, label := paste(node_1, node_2, sep = "_")]
  edge_table <- edge_table[,
    c("slope","midpoint_x", "midpoint_y", "distance", "weight", "weight_1", "weight_2", "weighted_midpoint_x", "weighted_midpoint_y", "raw_weight_1", "raw_weight_2") := {
      node_indices <- node_combinations[,.I];
      node_1 <- .SD[["node_1"]];
      node_2 <- .SD[["node_2"]];
      points <- nodes_xy[code == node_1 | code == node_2, ];
      edge_line <- edge_fun(x = as.matrix(points), fun = LearnGeom::CreateLinePoints);
      cent <- edge_center(as.matrix(points));
      distance <- node_distance(as.matrix(points));

      # browser();
      if(direction == 1) {
        weight_2 <- weights_matrix[node_indices[1], node_indices[2]];
        weight_1 <- weights_matrix[node_indices[2], node_indices[1]];
        raw_weight_2 <- raw_weights_matrix[node_indices[1], node_indices[2]];
        raw_weight_1 <- raw_weights_matrix[node_indices[2], node_indices[1]];
      }
      else {
        weight_1 <- weights_matrix[node_indices[1], node_indices[2]];
        weight_2 <- weights_matrix[node_indices[2], node_indices[1]];
        raw_weight_1 <- raw_weights_matrix[node_indices[1], node_indices[2]];
        raw_weight_2 <- raw_weights_matrix[node_indices[2], node_indices[1]];
      }
      weight_toward = ifelse(weight_2 > weight_1, 2, 1);

      # weight_diff = ifelse(weight_1 - weight_2 > 0,  (weight_1 - weight_2) / (weight_1 + weight_2), 0)
      weighted_center <- cent
      weight_diff = max(weight_1, weight_2) / (weight_1 + weight_2)
      if(!is.nan(weight_diff)) {
        if( weight_diff == 1 ) {
          if(direction == 1) {
            if(weight_toward == 1) {
              weighted_center <- as.matrix(points[2,]);
            }
            else {
              weighted_center <- as.matrix(points[1,]);
            }
          }
          else {
            if(weight_toward == 2) {
              weighted_center <- as.matrix(points[1,]);
            }
            else {
              weighted_center <- as.matrix(points[2,]);
            }
          }
        }
        else if ( abs(weight_1 - weight_2) > 0.000000001 ) {
          weighted_center <- edge_center_weighted(as.matrix(points), weight_diff, toward = weight_toward)
          if(weight_diff > 0) {
          }
        }
        else {
          # print("Why are we here?")
          # we're here b/c equal weights, but are sure about self connections?
        }
      }
      else {
        # Note: means there isn't a connection?
        # print("Why are we here too?")
      }

      list(
        slope = ifelse(edge_line[1] == "Inf", Inf, edge_line[1]),
        midpoint_x = cent[1],
        midpoint_y = cent[2],
        distance = distance,
        weight = weight_1 + weight_2,
        weight_1 = weight_1,
        weight_2 = weight_2,
        weighted_midpoint_x = weighted_center[1],
        weighted_midpoint_y = weighted_center[2],
        raw_weight_1 = raw_weight_1,
        raw_weight_2 = raw_weight_2
      )
    },
    by = seq(nrow(edge_table))
  ];

  attr(edge_table, "nodes_xy") <- nodes_xy
  attr(edge_table, "node_combinations") <- node_combinations

  edge_table
}

edge_paths <- function(x,
  edge_color = NULL, edge_end = "outward",
  edge_size_multiplier = 1,
  desaturate_edges_by = "weight", lighten_edges_by = desaturate_edges_by,
  fake_alpha = FALSE,
  scale_edges_to = c(0,1), scale_edges_from = c(0, max(x[, c("weight_1", "weight_2")])),
  path_as = "svg", edge_arrow = TRUE, edge_arrow_offset = 0.5, edge_arrow_width = 0.05,
  edge_arrow_saturation_multiplier = 0.75, edge_halo = 0.5, edge_halo_multiplier = 1.04,
  edge_arrows_to_show = "both", arrow_halo = edge_halo,
  sender_direction = 1, edge_arrow_direction = 2
) {
  edge_arrows_to_show <- match.arg(edge_arrows_to_show, c("both", "min", "max"))
  edges <- data.table::copy(x);
  edges_ <- data.table::melt(edges, measure.vars = c("weight_1", "weight_2"));
  edges_[variable == "weight_2", c("node_1", "node_2") := { list(node_2, node_1) }];

  data.table::setorder(edges, weight);
  data.table::setorder(edges_, value);

  nodes_xy <- attr(x, "nodes_xy");
  node_combinations <- attr(x, "node_combinations");

  if(is.null(edge_color)) {
    edge_color <- rep("#000000", nrow(edges_));
  }
  else if(length(edge_color) == 1) {
    edge_color <- rep(edge_color, nrow(edges_));
  }
  else if(length(edge_color) == 2) {
    edge_color <- edges_[, { ifelse(.SD[[paste0("raw_", variable)]] < 0, edge_color[2], edge_color[1]) }, by = 1:nrow(edges_)][[2]]
  }

  rgb_white <- grDevices::col2rgb("#FFFFFF", alpha = F)
  paths <- list();
  # for(e_row in seq.int(nrow(edges_))) {
  for(e_row in seq.int(nrow(edges_))) {
    edge <- edges_[e_row,]
    mid_point <- c(edge$weighted_midpoint_x, edge$weighted_midpoint_y);
    node_col <- "node_1";
    weight_col <- "value"; #ifelse(e == 1, 2, 1));

    node <- nodes_xy[code == edge[[node_col]]];
    line_xy <- matrix(c(edge$weighted_midpoint_x, node$x, edge$weighted_midpoint_y, node$y), ncol = 2);
    edge_line <- edge_fun(line_xy, fun = LearnGeom::CreateLinePoints);

    rnd <- 40;
    mp <- 1;
    weight <- edge[[weight_col]] # * edge_size_multiplier;
    if(weight > 0) { # && (mid_point != as.matrix(node))) {
      radius <- edge[[weight_col]] / 2;

      path <- list(type = "path", weight = weight, line = list( width = edge_halo, color = edge_color[e_row] )) #"#FF0000" ) )
      if(!is.null(desaturate_edges_by)) {
        this_color_hsv <- as.list(t(grDevices::rgb2hsv(grDevices::col2rgb(edge_color[e_row]))));
        names(this_color_hsv) <- c("h", "s", "v");

        if (is.function(desaturate_edges_by)) {
          this_color_hsv$s <- do.call(desaturate_edges_by, list(weight))
        }
        else if(desaturate_edges_by == "weight") {
          this_color_hsv$s <- scales::rescale(x = weight, to = scale_edges_to, from = scale_edges_from);
        }

        if(!is.null(lighten_edges_by)) {
          if (is.function(lighten_edges_by)) {
            this_color_hsv$alpha <- do.call(lighten_edges_by, list(weight))
          }
          else if(lighten_edges_by == "weight") {
            # this_color_hsv$alpha <- this_color_hsv$s;
            this_color_hsv$alpha <- weight;
          }
        }

        col_just_s <- grDevices::hsv(this_color_hsv$h, this_color_hsv$s, this_color_hsv$v)
        rgb_just_s <- col2rgb(col_just_s)

        path$opacity <- 1
        if(fake_alpha == TRUE) {
          # browser()
          # rgb_s_and_a_equiv <- rgb_white + (rgb_just_s - rgb_white) * (this_color_hsv$alpha)
          rgb_s_and_a_equiv <- lighten_color(rgb_just_s, alpha = this_color_hsv$alpha);
          rgb_s_and_a_edge <- darken_color(rgb_just_s, alpha = this_color_hsv$alpha , amt = edge_halo_multiplier);

          path$fillcolor <- grDevices::rgb(rgb_s_and_a_equiv[1,1], rgb_s_and_a_equiv[2,1], rgb_s_and_a_equiv[3,1], maxColorValue = 255)
          # browser(expr = { path$fillcolor == "#FFCFCF" })
          path$line$color <- grDevices::rgb(rgb_s_and_a_edge[1,1], rgb_s_and_a_edge[2,1], rgb_s_and_a_edge[3,1], maxColorValue = 255)
        }
        else {
          # path$opacity <- this_color_hsv$alpha
          path$fillcolor <- do.call(grDevices::hsv, this_color_hsv)
        }
      }
      else {
        path$fillcolor <- edge_color[e_row];
        rgb_just_s <- col2rgb(path$fillcolor)
        this_color_hsv <- as.list(t(grDevices::rgb2hsv(grDevices::col2rgb(edge_color[e_row]))));
        names(this_color_hsv) <- c("h", "s", "v");
        this_color_hsv$alpha <- weight;
      }

      if(edge_end == "outward") {
        rotate_on <- c(x = node$x, y = node$y);
        pinch_on <- c(edge$weighted_midpoint_x, edge$weighted_midpoint_y);
        if (edge_line[1] == 0) {
          perp_line <- LearnGeom::CreateLinePoints(c(node$x, node$y), c(node$x+0.001, 0));
        }
        else if(is.finite(edge_line[1])){
          perp_line <- LearnGeom::Rotate(edge_line, as.matrix(node), 90);
        }
        else {
          perp_line <- LearnGeom::CreateLinePoints(c(node$x, node$y), c(0.001, node$y));
        }
        perp_base <- LearnGeom::IntersectLineCircle(perp_line, as.matrix(node), radius * edge_size_multiplier);
        # Find points of circle, with center at node, that intersect along the edge line.
        # The one that is furthest from the midpoint is used as the extension for the curve
        # edge_ext <- LearnGeom::IntersectLineCircle(edge_line, as.matrix(node), edge[[weight_col]] * edge_size_multiplier);
        # ext_dist <- apply(edge_ext, 1, LearnGeom::DistancePoints, c(edge$weighted_midpoint_x, edge$weighted_midpoint_y));
        # ext_max <- which.max(ext_dist);
        # edge_ext_pt <- edge_ext[ext_max,];

        svg_curve_control_multiplier = 1.1;
        if(!is.finite(perp_line[1]) || perp_line[1] == "Inf" || perp_line[1] == Inf) {
          browser()
        }
        else if(perp_line[1] == 0) {
          perp_base_intercepts <- perp_base
          perp_base_intercepts[, 2] <- 0
          perp_line_1 <- LearnGeom::CreateLinePoints(perp_base[1,], perp_base_intercepts[1,]);
          perp_line_2 <- LearnGeom::CreateLinePoints(perp_base[2,], perp_base_intercepts[2,]);
          perp_adj <- rep((radius * svg_curve_control_multiplier) * edge_size_multiplier, 2) * c(1,-1)
          perp_line_1_ints <- matrix(perp_base[1,], nrow = 2, ncol = 2, byrow = TRUE)
          perp_line_2_ints <- matrix(perp_base[2,], nrow = 2, ncol = 2, byrow = TRUE)

          perp_line_1_ints[, 2] <- perp_line_1_ints[,2] + perp_adj
          perp_line_2_ints[, 2] <- perp_line_2_ints[,2] + perp_adj
        }
        else {
          perp_line_1 <- LearnGeom::Rotate(perp_line, perp_base[1,], 90)
          perp_line_2 <- LearnGeom::Rotate(perp_line, perp_base[2,], 90);

          perp_line_1_ints <- LearnGeom::IntersectLineCircle(perp_line_1, perp_base[1,], (radius * svg_curve_control_multiplier) * edge_size_multiplier);
          perp_line_2_ints <- LearnGeom::IntersectLineCircle(perp_line_2, perp_base[2,], (radius * svg_curve_control_multiplier) * edge_size_multiplier);
        }
        perp_line_1_dists <- apply(perp_line_1_ints, 1, LearnGeom::DistancePoints, c(edge$weighted_midpoint_x, edge$weighted_midpoint_y));
        perp_line_1_dists_max <- which.max(perp_line_1_dists);
        perp_line_1_pt <- perp_line_1_ints[perp_line_1_dists_max,];

        perp_line_2_dists <- apply(perp_line_2_ints, 1, LearnGeom::DistancePoints, c(edge$weighted_midpoint_x, edge$weighted_midpoint_y));
        perp_line_2_dists_max <- which.max(perp_line_2_dists);
        perp_line_2_pt <- perp_line_2_ints[perp_line_2_dists_max,];

        new_path <- paste0(
           "M", paste(mp * round(c(pinch_on[1], pinch_on[2]), rnd), collapse = " "),
          " L", paste(mp * round(c(perp_base[1,1], perp_base[1,2]), rnd), collapse = " "),
          " C", paste(mp * round(c(perp_line_1_pt[1], perp_line_1_pt[2]), rnd), collapse = " "), ", ",
                paste(mp * round(c(perp_line_2_pt[1], perp_line_2_pt[2]), rnd), collapse = " "), ", ",
                paste(mp * round(c(perp_base[2,1], perp_base[2,2]), rnd), collapse = " "),
          " Z"
        )
        path$path <- new_path;
        paths <- c(paths, list(path));
      }
      else if (edge_end == "inward") {
        rotate_on <- c(x = edge$weighted_midpoint_x, y = edge$weighted_midpoint_y)
        pinch_on <- c(x = node$x, y = node$y);
        if(is.finite(edge_line[1])){
          perp_line <- LearnGeom::Rotate(edge_line, as.matrix(rotate_on), 90);
        }
        else {
          perp_line <- LearnGeom::CreateLinePoints(rotate_on, c(0, rotate_on[2]));
        }
        perp_base <- LearnGeom::IntersectLineCircle(perp_line, as.matrix(rotate_on), radius * edge_size_multiplier);

        svg_curve_control_multiplier = 1.1;
        if(!is.finite(perp_line[1]) || perp_line[1] == "Inf" || perp_line[1] == Inf) {
          browser()
        }
        else if(perp_line[1] == 0) {
          perp_base_intercepts <- perp_base
          perp_base_intercepts[, 2] <- 0
          perp_line_1 <- LearnGeom::CreateLinePoints(perp_base[1,], perp_base_intercepts[1,]);
          perp_line_2 <- LearnGeom::CreateLinePoints(perp_base[2,], perp_base_intercepts[2,]);
          perp_adj <- rep((radius * svg_curve_control_multiplier) * edge_size_multiplier, 2) * c(1,-1)
          perp_line_1_ints <- matrix(perp_base[1,], nrow = 2, ncol = 2, byrow = TRUE)
          perp_line_2_ints <- matrix(perp_base[2,], nrow = 2, ncol = 2, byrow = TRUE)

          perp_line_1_ints[, 2] <- perp_line_1_ints[,2] + perp_adj
          perp_line_2_ints[, 2] <- perp_line_2_ints[,2] + perp_adj
        }
        else {
          perp_line_1 <- LearnGeom::Rotate(perp_line, perp_base[1,], 90)
          perp_line_2 <- LearnGeom::Rotate(perp_line, perp_base[2,], 90);

          perp_line_1_ints <- LearnGeom::IntersectLineCircle(perp_line_1, perp_base[1,], (radius * svg_curve_control_multiplier) * edge_size_multiplier);
          perp_line_2_ints <- LearnGeom::IntersectLineCircle(perp_line_2, perp_base[2,], (radius * svg_curve_control_multiplier) * edge_size_multiplier);
        }
        perp_line_1_dists <- apply(perp_line_1_ints, 1, LearnGeom::DistancePoints, c(edge$weighted_midpoint_x, edge$weighted_midpoint_y));
        perp_line_1_dists_max <- which.max(perp_line_1_dists);
        perp_line_1_pt <- perp_line_1_ints[perp_line_1_dists_max,];

        perp_line_2_dists <- apply(perp_line_2_ints, 1, LearnGeom::DistancePoints, c(edge$weighted_midpoint_x, edge$weighted_midpoint_y));
        perp_line_2_dists_max <- which.max(perp_line_2_dists);
        perp_line_2_pt <- perp_line_2_ints[perp_line_2_dists_max,];

        if(path_as == "svg") {
          new_path <- paste0(
             "M ", paste(mp * round(c(pinch_on[1], pinch_on[2]), rnd), collapse = " "),
            " L ", paste(mp * round(c(perp_base[1,1], perp_base[1,2]), rnd), collapse = " "),
            " L ", paste(mp * round(c(perp_base[2,1], perp_base[2,2]), rnd), collapse = " "),
            " Z "
          )
          path$path <- new_path;
          paths <- c(paths, list(path));

          paths <- c(paths, list(list(
            type = "path", opacity = 0.25,
            line = list( width = 1, color = "#000000", opacity = 0.5),
            path = paste0("M", paste(mp * round(c(perp_base[1,1], perp_base[1,2]), rnd), collapse = " "),
                         " L", paste(mp * round(c(perp_base[2,1], perp_base[2,2]), rnd), collapse = " "))
          )))
        }
        else {
          new_path <- as.data.frame(matrix(
            c(mp * round(c(pinch_on[1], pinch_on[2]), rnd),
              mp * round(c(perp_base[1,1], perp_base[1,2]), rnd),
              mp * round(c(perp_base[2,1], perp_base[2,2]), rnd),
              mp * round(c(pinch_on[1], pinch_on[2]), rnd)
            ),
            ncol = 2, byrow = TRUE, dimnames = list(NULL, c("x", "y")))
          )
          path$path <- new_path;
          paths <- c(paths, list(path));
        }
      }
      else {
        stop("edge_end parameter unknown")
      }

      if(edge_arrow == TRUE &&
         (
            edge_arrows_to_show == "both" ||
           # (do.call(paste0("which.", edge_arrows_to_show), list(edge[, c("weight_1", "weight_2")])) == e)
           edges_[label == edge$label, .SD[which.max(value)]]$node_1 == edge$node_1
         )
      ) {
        x_axis <- LearnGeom::CreateLinePoints(c(-10,0), c(10,0))
        perp_line_angle <- LearnGeom::LinesAngles(perp_line, x_axis)
        if(perp_line_angle[1] == 0) {
          mid_point[1] <- mid_point[1] - 0.0001
        }
        arrow_offset <- LearnGeom::MidPoint(mid_point, as.matrix(node))
        arrow_offset_hgt <- CreatePointByDistance(as.matrix(node), arrow_offset, dist = weight * edge_arrow_width);
        arrow_offset_hgt2 <- CreatePointByDistance(as.matrix(node), arrow_offset, dist = weight * (edge_arrow_width*2));

        arrow_offset_line <- MoveLineAlong(arrow_offset[1], arrow_offset[2], perp_line[1])
        arrow_offset_hgt_line <- MoveLineAlong(arrow_offset_hgt[1], arrow_offset_hgt[2], perp_line[1])
        arrow_offset_hgt2_line <- MoveLineAlong(arrow_offset_hgt2[1], arrow_offset_hgt2[2], perp_line[1])
        edge_line_1 <- LearnGeom::CreateLinePoints(perp_base[1,], mid_point);
        edge_line_2 <- LearnGeom::CreateLinePoints(perp_base[2,], mid_point);

        # if(edge_line_1[1] == 0) {
        #   # add for offset
        #   print("ER...Tell Cody you saw me.");
        #   # arrow_int_1 <- LearnGeom::IntersectLines(edge_line_1, arrow_offset_hgt_line);
        #   # arrow_int_2 <- LearnGeom::IntersectLines(edge_line_2, arrow_offset_hgt_line);
        #   # arrow_int2_1 <- LearnGeom::IntersectLines(edge_line_1, arrow_offset_hgt2_line);
        #   # arrow_int2_2 <- LearnGeom::IntersectLines(edge_line_2, arrow_offset_hgt2_line);
        # }
        # else {
          offset_int_1 <- LearnGeom::IntersectLines(edge_line_1, arrow_offset_line);
          offset_int_2 <- LearnGeom::IntersectLines(edge_line_2, arrow_offset_line);
          arrow_int_1 <- LearnGeom::IntersectLines(edge_line_1, arrow_offset_hgt_line);
          arrow_int_2 <- LearnGeom::IntersectLines(edge_line_2, arrow_offset_hgt_line);
          arrow_int2_1 <- LearnGeom::IntersectLines(edge_line_1, arrow_offset_hgt2_line);
          arrow_int2_2 <- LearnGeom::IntersectLines(edge_line_2, arrow_offset_hgt2_line);
        # }

        arr_path <- list(line = list(width = 0));
        if ( edge_arrow_direction == 1 ) {
          new_path <- paste0(
             "M", paste(mp * round(c(arrow_int_1), rnd), collapse = " "),
            " L", paste(mp * round(c(arrow_offset), rnd), collapse = " "),
            " L", paste(mp * round(c(arrow_int_2), rnd), collapse = " "),
            " L", paste(mp * round(c(arrow_int2_2), rnd), collapse = " "),
            " L", paste(mp * round(c(arrow_offset_hgt), rnd), collapse = " "),
            " L", paste(mp * round(c(arrow_int2_1), rnd), collapse = " "),
            " Z"
          )
        }
        else {
          new_path <- paste0(
             "M", paste(mp * round(c(offset_int_1), rnd), collapse = " "),
            " L", paste(mp * round(c(arrow_offset_hgt), rnd), collapse = " "),
            " L", paste(mp * round(c(offset_int_2), rnd), collapse = " "),
            " L", paste(mp * round(c(arrow_int_2), rnd), collapse = " "),
            " L", paste(mp * round(c(arrow_offset_hgt2), rnd), collapse = " "),
            " L", paste(mp * round(c(arrow_int_1), rnd), collapse = " "),
            " Z"
          )
        }

        arr_path$path <- new_path;
        arr_path$fillcolor <- "#FFFFFF";
        arr_path$line$color <- "#FFFFFF";
        arr_path$line$width <- arrow_halo;
        arr_path$opacity <- 1
        if(!is.null(this_color_hsv)) {
          # add white
          if(edge_arrow_saturation_multiplier < 1) {
            # rgb_s_and_a_equiv_ext <- rgb_white + (rgb_just_s - rgb_white) * (this_color_hsv$alpha * edge_arrow_saturation_multiplier)
            rgb_s_and_a_equiv_ext <- lighten_color(rgb_just_s, this_color_hsv$alpha, edge_arrow_saturation_multiplier);
          }

          # add shade
          else {
            # rgb_s_and_a_equiv_ext <- (rgb_white + (rgb_just_s - rgb_white) * (this_color_hsv$alpha)) * (1 - -(1 - edge_arrow_saturation_multiplier))
            rgb_s_and_a_equiv_ext <- darken_color(rgb_just_s, edge_arrow_saturation_multiplier, this_color_hsv$alpha)
          }

          col_s_and_a_equiv_ext <- grDevices::rgb(rgb_s_and_a_equiv_ext[1,1], rgb_s_and_a_equiv_ext[2,1], rgb_s_and_a_equiv_ext[3,1], maxColorValue = 255)
          arr_path$fillcolor <- col_s_and_a_equiv_ext
        }

        paths <- c(paths, list(arr_path));
      }
    }

    else {
      # cat("Zero edge: ", edge$node_1, ", ", edge$node_2, " - ", weight, "\n");
      # print(edge);
    }
  }
  paths
}


base_plot <- function(
  width = 600, height = width,
  autosize = (is.null(width) && is.null(height)),
  max_axis_marker = 3,
  title = "ENA Plot", ...
) {
  p <-
    plotly::plot_ly(
      mode = "markers+lines", type = "scatter",
      width = width, height = height
    )
  ;

  axis_list <- list(
    type = "linear",
    showgrid = F,
    title = list(text = ""),
    showticklabels = FALSE,
    range = c(-max_axis_marker, max_axis_marker)
  );
  x_axis <- axis_list;
  x_axis$scaleanchor = "y";
  y_axis <- axis_list;

  p <- plotly::layout( p = p, title = title,
    autosize = autosize,
    showlegend = FALSE,
    xaxis = x_axis,
    yaxis = y_axis
  );

  class(p) <- c("enaplot", class(p));

  return(p);
}

#' Plot Function
#'
#' This function generates a plot based on the provided inputs.
#'
#' @param x A numeric vector or object to be plotted.
#' @param ... Additional arguments passed to the plotting method.
#' @param y An optional numeric vector for y-values.
#'
#' @return A plot object or visualization.
#' @export
'plot.ena.ordered.set' <- function(x, y, ...) {
  p <- base_plot(...)

  attr(p, "model") <- x;

  args <- list(...);
  dimensions <- colnames(as.matrix(x$points))[1:2];
  if(!is.null(args$dimensions)) {
    dimensions <- args$dimensions;
  }
  attr(p, "dimensions") <- dimensions;

  p
}

#' Plot ENA nodes
#'
#' Plots the nodes of an ENA model using plotly, with options for node size, color, labels, and self-connections.
#'
#' @param x An ENA model object or plotly object with model attached.
#' @param ... Additional arguments passed to plotting functions.
#' @param weighted_nodes Logical; if TRUE, node sizes are weighted (default: TRUE).
#' @param node_size_multiplier Numeric; scales node sizes (default: 1).
#' @param node_color Character; color(s) for nodes (default: "#3E4750").
#' @param node_sizes Numeric; custom node sizes (default: 1).
#' @param node_direction Integer; direction for node size calculation (default: 2).
#' @param self_connections Logical; if TRUE, show self-connections (default: weighted_nodes).
#' @param self_connection_color Character; color(s) for self-connections (default: "#000000").
#' @param nodes_include_self Logical; if TRUE, include self-connections in node size (default: TRUE).
#' @param node_halo Numeric; width of node halo (default: 0).
#' @param self_connection_halo Numeric; width of self-connection halo (default: 0).
#' @param node_labels Logical; if TRUE, show node labels (default: TRUE).
#' @param node_position_multiplier Numeric; scales node positions (default: 1).
#' @param node_saturation_multiplier Numeric; controls node color saturation (default: 1).
#' @param sender_direction Integer; direction for sender (default: 1).
#' @param weights Optional; edge weights to use (default: from model).
#' @param nodes Optional; node positions to use (default: from model).
#'
#' @return A plotly object with ENA nodes drawn.
#' @export
nodes <- function(
  x, ...,
  weighted_nodes = TRUE, node_size_multiplier = 1,
  node_color = "#3E4750", node_sizes = 1, node_direction = 2,
  self_connections = weighted_nodes, self_connection_color = "#000000",
  nodes_include_self = TRUE, node_halo = 0, self_connection_halo = 0,
  node_labels = TRUE, node_position_multiplier = 1, node_saturation_multiplier = 1,

  #
  sender_direction = 1,

  # parameters pulled from model `x` that can be overridden
  weights = NULL, nodes = NULL
) {
  p <- NULL;
  if(is(x, "plotly")) {
    p <- x;
    x <- attr(x = p, which = "model");
    if(is.null(weights)) {
      weights <- attr(p, "edge_weights");
    }
  }
  if(is.null(weights)) {
    weights <- colMeans(as.matrix(x$line.weights$ENA_DIRECTION$response));
  }
  if(is.null(nodes)) {
    nodes <- data.table::copy(x$rotation$nodes);
  }
  else {
    nodes <- data.table::copy(nodes);
  }
  raw_weights <- weights;
  raw_sq_weights_raw <- to_square(raw_weights);
  raw_sq_weights <- raw_weights * node_size_multiplier;

  weights <- abs(weights);

  dimensions <- attr(p, "dimensions");

  # if(is.null(scale_edges_from)) {
  #   if(!is.null(model)) {
  #     scale_edges_from <- c(0, max(as.matrix(model$line.weights$ENA_DIRECTION$response)))
  #   }
  #   else {
  #     scale_edges_from <- c(0, max(raw_weights))
  #   }
  # }

  sq_weights_raw <- to_square(weights);
  sq_weights <- sq_weights_raw * node_size_multiplier;

  if(nodes_include_self == TRUE) {
    node_sizes_raw <- (apply(sq_weights_raw, node_direction, sum));
    node_sizes <- (apply(sq_weights, node_direction, sum));
  }
  else {
    node_sizes_raw <- (apply(sq_weights_raw, node_direction, sum) - diag(sq_weights_raw));
    node_sizes <- (apply(sq_weights, node_direction, sum) - diag(sq_weights));
  }
  node_sizes[node_sizes == 0] <- 0.001

  nodes[, (seq(2, ncol(nodes))) := lapply(.SD, `*`, node_position_multiplier) , .SDcols = (seq(2, ncol(nodes)))];
  dim_indices <- which(colnames(as.matrix(nodes)) %in% dimensions);
  edge_table <- create_edge_matrix(
    NULL, weights = raw_weights, nodes = nodes, direction = sender_direction, dims = dim_indices,
    node_position_multiplier = 1
  );
  nodes_xy <- attr(edge_table, "nodes_xy");

  annotations <- list();
  node_traces <- list();
  if(weighted_nodes) {
    # Node definition
    shape = list(
      type = "circle",
      xanchor = "x",
      yanchor = "y",
      xref = "x",
      yref = "y",
      line = list (color = "#FFFFFF", width = node_halo)
    )

    # ___ Outer nodes ----
    if(length(node_color) == 1) {
      node_color <- rep(node_color, nrow(sq_weights));
    }
    if(self_connections == TRUE) {
      if(length(self_connection_color) == 1) {
        self_connection_color <- rep(self_connection_color, nrow(sq_weights));
      }
      else if(length(self_connection_color) == 2) {
        self_connection_color <- self_connection_color[((diag(raw_sq_weights_raw) < 0) * 1) + 1];
      }
    }
    for(i in seq(nrow(sq_weights))) {
      if(node_sizes[i] > 0) {
        shape$line$width = node_halo;
        node <- nodes_xy[i,];
        node_radius <- node_sizes[i] / 2;
        edge_line <- edge_fun(rbind(as.matrix(node),c(0,0)), fun = LearnGeom::CreateLinePoints);
        node_pts <- LearnGeom::IntersectLineCircle(edge_line, as.matrix(node), node_radius);
        shape$x0 <- node$x - (node_radius);
        shape$y0 <- node$y - (node_radius);
        shape$x1 <- node$x + (node_radius);
        shape$y1 <- node$y + (node_radius);
        shape$fillcolor <- paste0(node_color[i], alpha_to_hex(node_saturation_multiplier)); #hsv(0, 0, 0.1, 1); #edge_sizes[i]);

        node_traces <- c(node_traces, list(shape));

        if(node_labels == TRUE) {
          annotations <- c(annotations, list(list(
            x = shape$x1,
            y = shape$y1,
            text = node$code,
            xref = "x", yref = "y",
            showarrow = FALSE, arrowhead = 0,
            ax = 0, ay = -60
          )))
        }

        # ___ Self nodes ----
        if(self_connections) {
          shape$line$width = self_connection_halo;
          shape$x0 <- node$x - (diag(sq_weights)[i] / 2);
          shape$y0 <- node$y - (diag(sq_weights)[i] / 2);
          shape$x1 <- node$x + (diag(sq_weights)[i] / 2);
          shape$y1 <- node$y + (diag(sq_weights)[i] / 2);
          shape$fillcolor <- hsv(0, 0, 1, 1);
          node_traces <- c(node_traces, list(shape));

          shape$x0 <- node$x - (diag(sq_weights)[i] / 2);
          shape$y0 <- node$y - (diag(sq_weights)[i] / 2);
          shape$x1 <- node$x + (diag(sq_weights)[i] / 2);
          shape$y1 <- node$y + (diag(sq_weights)[i] / 2);

          this_color_hsv <- as.list(t(grDevices::rgb2hsv(grDevices::col2rgb(self_connection_color[i]))));
          names(this_color_hsv) <- c("h", "s", "v");
          this_color_hsv$alpha <- diag(sq_weights_raw)[i];
          shape$fillcolor <- do.call(grDevices::hsv, this_color_hsv); #hsv(1, 1, 1, edge_sizes[i]);
          node_traces <- c(node_traces, list(shape));
        }
      }
    }

    attr(x = p, which = "annotations") <- c(attr(x = p, which = "annotations"), annotations);
    attr(x = p, which = "shapes") <- c(attr(x = p, which = "shapes"), node_traces);
  }

  # __ Non-weighted Nodes ----
  else {
    nodes_xy$size <- 15;
    p <- plotly::add_trace(p = p, data = nodes_xy,
      x = ~x, y = ~y,
      text = ~code,
      type = "scatter",
      mode = "markers+text",
      textposition = "bottom center",
      name = "Nodes",
      marker = list(
        sizemode = "diameter",
        size = ~size,
        color = "black",
        line = list(width = 0)
      )
    );
  }

  return(p);
}

#' Plot ENA edges
#'
#' Plots the edges (connections) between nodes in an ENA model using plotly, with options for color, width, direction, and arrowheads.
#'
#' @param x An ENA model object or plotly object with model attached.
#' @param ... Additional arguments passed to plotting functions.
#' @param show_edge_center Logical; if TRUE, show edge centers (default: FALSE).
#' @param units Optional; units to include in the plot.
#' @param edge_color Character vector of edge colors (default: c("#FF0000", "#0000FF")).
#' @param edge_end Character; "outward" or "inward" (default: "outward").
#' @param desaturate_edges_by Function or string; controls edge desaturation (default: sqrt).
#' @param lighten_edges_by Function or string; controls edge lightening (default: same as desaturate_edges_by).
#' @param multiplier Numeric; scales edge sizes (default: 1).
#' @param edge_size_multiplier Numeric; scales edge width (default: multiplier).
#' @param edges_behind_nodes Logical; if TRUE, draw edges behind nodes (default: TRUE).
#' @param sender_direction Integer; direction for sender (default: 1).
#' @param scale_edges_to Numeric vector; range to scale edge weights to (default: c(0,1)).
#' @param scale_edges_from Numeric vector; range of edge weights to scale from (default: NULL).
#' @param edge_arrow Logical; if TRUE, draw arrowheads (default: TRUE).
#' @param edge_arrow_offset Numeric; offset for arrowhead placement (default: 2).
#' @param edge_arrow_width Numeric; width of arrowhead (default: 0.1).
#' @param edge_arrows_to_show Character; which arrows to show (default: "max").
#' @param edge_arrow_saturation_multiplier Numeric; controls arrowhead color saturation (default: 1.1).
#' @param fake_alpha Logical; if TRUE, simulate alpha blending (default: TRUE).
#' @param edge_halo Numeric; width of edge halo (default: 0.5).
#' @param edge_halo_multiplier Numeric; multiplier for edge halo (default: 1.1).
#' @param arrow_halo Numeric; width of arrow halo (default: 0).
#' @param edge_arrow_direction Integer; direction for arrowhead (default: 2).
#' @param weights Optional; edge weights to use (default: from model).
#' @param nodes Optional; node positions to use (default: from model).
#' @param node_position_multiplier Numeric; scales node positions (default: 1).
#'
#' @return A plotly object with ENA edges drawn.
#' @export
edges <- function(
  x, ...,
  show_edge_center = FALSE, units = NULL,
  edge_color = c("#FF0000","#0000FF"), edge_end = "outward",
  desaturate_edges_by = sqrt, lighten_edges_by = desaturate_edges_by,
  multiplier = 1,
  edge_size_multiplier = multiplier, edges_behind_nodes = TRUE, sender_direction = 1,
  scale_edges_to = c(0,1), scale_edges_from = NULL,
  edge_arrow = TRUE, edge_arrow_offset = 2, edge_arrow_width = 0.1, edge_arrows_to_show = "max",
  edge_arrow_saturation_multiplier = 1.1, fake_alpha = TRUE, edge_halo = 0.5, edge_halo_multiplier = 1.1,
  arrow_halo = 0, edge_arrow_direction = 2, node_position_multiplier = 1,

  # parameters pulled from model `x` that can be overridden
  weights = NULL, nodes = NULL
) {
  .SD <- NULL;

  p <- NULL;
  if(is(x, "plotly")) {
    p <- x;
    x <- attr(x = p, which = "model");
  }
  if(is.null(weights)) {
    weights <- colMeans(as.matrix(x$line.weights$ENA_DIRECTION$response));
  }
  else if (!is.numeric(weights)) {
    weights <- colMeans(as.matrix(weights[ENA_DIRECTION == 'response',]));
  }
  if(is.null(nodes)) {
    nodes <- data.table::copy(x$rotation$nodes);
  }
  else {
    nodes <- data.table::copy(nodes);
  }
  nodes[, (seq(2, ncol(nodes))) := lapply(.SD, `*`, node_position_multiplier) , .SDcols = (seq(2, ncol(nodes)))];
  raw_weights <- weights;
  raw_sq_weights_raw <- to_square(raw_weights);
  raw_sq_weights <- raw_weights * multiplier;

  weights <- abs(weights);

  if(is.null(scale_edges_from)) {
    if(!is.null(model)) {
      scale_edges_from <- c(0, max(as.matrix(x$line.weights$ENA_DIRECTION$response)))
    }
    else {
      scale_edges_from <- c(0, max(raw_weights))
    }
  }

  sq_weights_raw <- to_square(weights);
  sq_weights <- sq_weights_raw * multiplier;

  edge_table <- create_edge_matrix(NULL, weights = raw_weights, nodes = nodes, direction = sender_direction);
  nodes_xy <- attr(edge_table, "nodes_xy");

  edge_path_list <- list();
  nodes_xy <- attr(edge_table, "nodes_xy");
  node_combinations <- attr(edge_table, "node_combinations");

  edge_path_list <- edge_paths(
    edge_table, edge_color = edge_color, edge_end = edge_end,
    desaturate_edges_by = desaturate_edges_by, lighten_edges_by = lighten_edges_by,
    edge_size_multiplier = edge_size_multiplier,
    scale_edges_to = scale_edges_to, scale_edges_from = scale_edges_from,
    edge_arrow = edge_arrow, edge_arrow_offset = edge_arrow_offset,
    edge_arrow_width = edge_arrow_width, edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    fake_alpha = fake_alpha, edge_halo = edge_halo, edge_arrows_to_show = edge_arrows_to_show,
    arrow_halo = arrow_halo, edge_halo_multiplier = edge_halo_multiplier,
    sender_direction = sender_direction,
    edge_arrow_direction = edge_arrow_direction
  );

  if(show_edge_center) {
    p <- plotly::add_markers(p = p, data = edge_table, x = ~weighted_midpoint_x, y = ~weighted_midpoint_y, marker = list(color = "black"), name = ~label)
    p <- plotly::add_markers(p = p, data = edge_table, x = ~midpoint_x, y = ~midpoint_y, marker = list(color = "black"), name = ~label)
  }

  # p <- plotly::layout(p = p, shapes = edge_path_list);
  attr(x = p, which = "shapes") <- c(attr(x = p, which = "shapes"), edge_path_list);
  attr(x = p, which = "edge_weights") <- raw_weights;

  return (p);
}

unit_vectors <- function(points, dims = c(1,2)) {
  if(is(points, "list")) {
    meta_cols <- which(sapply(points[[1]], is, "ena.metadata"))
    dim_cols <- which(sapply(points[[1]], is, "ena.dimension"))[dims]
    dim_col_names <- colnames(points[[1]])[dim_cols]
  }
  else {
    meta_cols <- which(sapply(points, is, "ena.metadata"))
    dim_cols <- which(sapply(points, is, "ena.dimension"))[dims]
    dim_col_names <- colnames(points)[dim_cols]
  }

  merge_by = c("ENA_UNIT")
  if("_GRP" %chin% colnames(points)) {
    dim_col_names <- c("_GRP", dim_col_names);
    merge_by = c(merge_by, "_GRP");
  }

  from_rows <- points[ENA_DIRECTION == "ground", c("ENA_UNIT", dim_col_names), with = FALSE]
  to_rows <- points[ENA_DIRECTION == "response", c("ENA_UNIT", dim_col_names), with = FALSE]

  sd_cols <- paste(rep(dim_col_names,2), c(rep("from",2), rep("to",2)), sep = ".");
  if(nrow(from_rows) > 0) {
    y <- merge(from_rows, to_rows, by = merge_by, suffixes = c(".from", ".to"));
    y[, distance := {
      LearnGeom::DistancePoints(c(.SD[[sd_cols[1]]], .SD[[sd_cols[2]]]), c(.SD[[sd_cols[3]]], .SD[[sd_cols[4]]]))
    }, .SDcols = c(sd_cols), by = 1:nrow(y)];
    y
  }
  else {
    to_rows;
  }
}

#' Plot ENA units (points and vectors)
#'
#' Plots the units (points) and optionally vectors/arrows for each unit in an ENA model using plotly.
#'
#' @param x An ENA model object or plotly object with model attached.
#' @param ... Additional arguments passed to plotting functions.
#' @param units_as_vectors Logical; if TRUE, plot unit vectors/arrows (default: FALSE).
#' @param unit_vectors_towards Character; "response" or "ground" (default: "response").
#' @param mean_vector_width Numeric; width of mean vector arrow (default: 5).
#' @param mean_vector_arrow_size Numeric; size of mean vector arrowhead (default: 0.5).
#' @param unit_vector_width Numeric; width of unit vector arrows (default: 3).
#' @param unit_vector_arrow_size Numeric; size of unit vector arrowheads (default: 0.7).
#' @param show_points Logical; if TRUE, show points (default: TRUE).
#' @param point_shape Character; shape of points (default: "circle").
#' @param dynamic_response_point_size Logical; if TRUE, scale response point size by distance (default: FALSE).
#' @param show_mean Logical; if TRUE, show mean vector/point (default: TRUE).
#' @param mean_value Optional; mean value to use (default: NULL).
#' @param points_color Character; color for points (default: "#FF0000").
#' @param saturate_vectors_by Numeric or function; controls vector color saturation (default: 0.2).
#' @param point_position_multiplier Numeric; scales point positions (default: 1).
#' @param distance_scaled_by Function; rescales distances for arrow size (default: scales::rescale).
#' @param with_lines Logical; if TRUE, connect points with lines (default: FALSE).
#' @param with_ci Logical; if TRUE, show confidence intervals (default: TRUE).
#' @param points Optional; points to use (default: from model).
#'
#' @return A plotly object with ENA units (points and/or vectors) drawn.
#' @export
units <- function(
  x, ...,
  units_as_vectors = FALSE, unit_vectors_towards = "response",
  mean_vector_width = 5, mean_vector_arrow_size = 0.5,
  unit_vector_width = 3, unit_vector_arrow_size = 0.7,
  show_points = TRUE, point_shape = "circle",
  dynamic_response_point_size = FALSE,
  show_mean = TRUE,
  mean_value = NULL,
  points_color = "#FF0000",
  saturate_vectors_by = 0.2,
  point_position_multiplier = 1,
  distance_scaled_by = function(x) { scales::rescale(x, to = c(0.3, 1.5), from = c(0, 2)) },
  with_lines = FALSE,
  with_ci = TRUE,


  # parameters pulled from model `x` that can be overridden
  points = NULL
) {
  p <- NULL;
  if(is(x, "plotly")) {
    p <- x;
    x <- attr(x = p, which = "model");
  }
  if(is.null(points)) {
    points <- x$points;
  }

  point_dim_cols <- which(rENA::find_dimension_cols(points));
  points[,point_dim_cols] <- points[,c(point_dim_cols),with = FALSE] * point_position_multiplier;

  dimensions <- attr(p, "dimensions");
  dim_indices <- which(colnames(as.matrix(points)) %in% dimensions);


  unit_table <- unit_vectors(points, dims = dim_indices); # * scale_points_by;
  has_group <- "_GRP" %chin% colnames(unit_table);

  if(length(unique(points$ENA_DIRECTION)) > 1) {
    dim_cols <- paste(rep(dimensions,2), c(rep("from",2), rep("to",2)), sep = ".");
  }
  else {
    dim_cols <- dimensions;
  }
  if(is.null(mean_value)) {
    # dim_cols <- c("SVD1.from", "SVD2.from", "SVD1.to", "SVD2.to");
    mean_table <- if(has_group) {
      unit_table[, lapply(.SD, mean), .SDcols = dim_cols, by = "_GRP"];
    } else {
      unit_table[, lapply(.SD, mean), .SDcols = dim_cols];
    }
  }
  else {
    browser()
    mean_table <- data.table::data.table(
      SVD1.from = mean_value[1],
      SVD2.from = mean_value[2],
      SVD1.to = mean_value[3],
      SVD2.to = mean_value[4]
    )
  }

  # browser()

  if(length(unique(points$ENA_DIRECTION)) > 1) {
    mean_table[ , c("slope", "intercept", "center.x", "center.y") := {
      l = LearnGeom::CreateLinePoints(
        as.numeric(.SD[,c(dim_cols[1:2]), with = FALSE]),
        as.numeric(.SD[,c(dim_cols[3:4]), with = FALSE])
      );
      cen = LearnGeom::MidPoint(
        as.numeric(.SD[, c(dim_cols[1:2]), with = FALSE]),
        as.numeric(.SD[, c(dim_cols[3:4]), with = FALSE])
      );
      list(slope = l[1], intercept = l[2], center.x = cen[1], center.y = cen[2])
    }, by = 1:nrow(mean_table)]
  }

  rgb_cols <- col2rgb(points_color);
  hsv_cols <- apply(rgb_cols, 2, function(rgb_col) rgb2hsv(rgb_col[1], rgb_col[2], rgb_col[3]) );

  if(show_points) {

    if(length(unique(points$ENA_DIRECTION)) > 1) {
      units_by = if(has_group) c("ENA_UNIT", "_GRP") else c("ENA_UNIT")
      if(has_group) {
        unit_table[, c("ang_diff", "center.x", "center.y", "color") := {
          mt <- mean_table[mean_table$`_GRP` == .BY$`_GRP`,]
          list(
            LearnGeom::Angle(
              as.numeric(mt[, c(dim_cols[3:4]), with = FALSE]),
              as.numeric(mt[, c(dim_cols[1:2]), with = FALSE]),
              as.numeric(.SD[1, c(dim_cols[3:4]), with = FALSE])
            ),
            mt$center.x,
            mt$center.y,
            points_color[.BY$`_GRP`]
          )
        }, by = units_by]
      }
      else {
        unit_table[, c("ang_diff", "center.x", "center.y","color") := {
          list(
            LearnGeom::Angle(as.numeric(mean_table[1,3:4, with = FALSE]), as.numeric(mean_table[1,1:2, with = FALSE]), as.numeric(.SD[1,3:4, with = FALSE])),
            mean_table$center.x,
            mean_table$center.y,
            points_color
          )
        }, by = "ENA_UNIT"];
      }

      mean_table[, c("length") := {
        LearnGeom::DistancePoints(
          as.numeric(.SD[, c("center.x", "center.y"), with = FALSE]),
          as.numeric(.SD[, c(dim_cols[3:4]), with = FALSE])
        )
      }, by = 1:nrow(mean_table)]
    }
    else {
      if(has_group) {
        warning("Colors should be updated.")
      }
      else {
        unit_table$color = points_color;
      }
    }

    if(is.function(saturate_vectors_by)) {
      # saturate_vectors_by(unit_table$ang_diff, abs(unit_table$distance - mean_length))
      browser()
      units_by = if(has_group) c("ENA_UNIT", "_GRP") else c("ENA_UNIT")
      unit_table$ang_diff_scaled <- 1
    }
    else {
      unit_table$ang_diff_scaled <- saturate_vectors_by
    }

    unit_table[, color_adjusted := {
      adjusted_color <- lighten_color(.SD$color, alpha = .SD$ang_diff_scaled)
      grDevices::rgb(adjusted_color[1], adjusted_color[2], adjusted_color[3], maxColorValue = 255)
    }, by = 1:nrow(unit_table)]

    unit_points <- NULL;
    if(units_as_vectors == TRUE) {
      data.table::setorder(unit_table, -"distance")

      if(!is.null(point_shape)) {
        p <- plotly::add_markers(
          p = p,
          data = unit_table, x = ~center.x, y = ~center.y,
          marker = list(color = ~color_adjusted, shape = point_shape)
        );
      }

      arrow_size <- if(dynamic_response_point_size) {
        unit_table$distance_scaled <- distance_scaled_by(unit_table$distance);
        unit_table$distance_scaled
      }
      else {
        unit_vector_arrow_size
      }

      point_towards <- NULL
      if ( unit_vectors_towards == "ground" ) {
        point_towards <- unit_table[, c(dim_cols[1:2]), with = FALSE]
      }
      else {
        point_towards <- unit_table[, c(dim_cols[3:4]), with = FALSE]
      }
      p <- plotly::add_annotations(
        p = p,
        data = unit_table,
        # x = ~SVD1.to, y = ~SVD2.to,
        x = point_towards[[1]], y = point_towards[[2]],
        arrowsize = arrow_size,
        arrowwidth = unit_vector_width,
        xref = "x", yref = "y",
        axref = "x", ayref = "y",
        text = "", showarrow = TRUE,
        ax = ~center.x, ay = ~center.y,
        arrowcolor = ~color_adjusted
      );
    }

    # Not arrows
    else {
      if ( unit_vectors_towards == "ground" ) {
        unit_points <- unit_table[, c(dim_cols[1:2]), with = FALSE]
      }
      else {
        unit_points <- unit_table[, c(which(colnames(unit_table) %in% dim_cols)), with = FALSE]
      }
      p <- plotly::add_markers(
        p = p,
        # data = unit_table, x = ~SVD1.to, y = ~SVD2.to,
        data = unit_table, x = unit_points[[1]], y = unit_points[[2]],
        marker = list(color = ~color, shape = point_shape),
        text = ~ENA_UNIT,
        hovertemplate = paste('(%{x},%{y})', '<b>%{text}</b>')
      )
      if(with_lines == TRUE) {
        if ( unit_vectors_towards == "ground" ) {
          point_towards_x <- c("center.x", dim_cols[1])
          point_towards_y <- c("center.y", dim_cols[2])
        }
        else {
          point_towards_x <- c("center.x", dim_cols[3])
          point_towards_y <- c("center.y", dim_cols[4])
        }
        for( i in 1:nrow(unit_table) ) {
          p <- plotly::add_trace(
            p = p,
            type = "scatter",
            mode = "lines",
            data = data.frame(
              X1 = unlist(unit_table[i, c(point_towards_x), with = FALSE]),
              X2 = unlist(unit_table[i, c(point_towards_y), with = FALSE]),
              row.names = NULL
            ),
            x = ~X1, y = ~X2,
            line = list(
              color = unit_table[i,]$color_adjusted,
              width = unit_vector_width
            )
          )
        }
      }
    }
  }
  if(show_mean){
    mean_table$color <- points_color;
    error = list(
      x = list(color = points_color[1], visible=T, type="data"),
      y = list(color = points_color[1], visible=T, type="data")
    );
    int.values = NULL;
    if(with_ci == TRUE) {
      if(has_group) {
        int.values <- unit_table[, lapply(.SD, function(x) { t.test(x, conf.level = 0.95)$conf.int }), .SDcols = dim_cols, by = "_GRP"];

        error$x$array = unlist(int.values[c(1,3), c(4)]);
        error$y$array = unlist(int.values[c(2,4), c(5)]);
        error$x$color = "black"; #points_color;
        error$y$color = "black"; #points_color;
      } else {
        int.values <- matrix(
          c(
            as.vector(stats::t.test(points[,point_dim_cols[1],with=FALSE], conf.level = 0.95)$conf.int),
            as.vector(stats::t.test(points[,point_dim_cols[2],with=FALSE], conf.level = 0.95)$conf.int)
          ),
          ncol=2, byrow = TRUE
        );
        #error$x$array = int.values[,1] Carl thinks these two lines are not right
        #error$y$array = int.values[,2]
        error$x$array = (int.values[1, 2]-int.values[1, 1])/2 # have to compute radius
        error$y$array = (int.values[2, 2]-int.values[2, 1])/2 # if you use two numbers, it may just use the first number.

      }
    }

    if(units_as_vectors == TRUE) {
      p <- plotly::add_markers(
        p = p,
        data = mean_table,
        x = mean_table$center.x,
        y = mean_table$center.y,
        error_x = error$x, error_y = error$y,
        marker = list(color = ~color, symbol = "square", size = 10)
      );
    }
    else {
      # mean_cols <- dim_cols[seq.int(ifelse(unit_vectors_towards == "ground", 1, 3), length.out = 2)];
      mean_cols <- colnames(mean_table)[1:2]; #colnames(unit_points);
      p <- plotly::add_markers(
        p = p,
        data = mean_table,
        # x = mean_table$center.x,
        # y = mean_table$center.y,
        x = mean_table[[mean_cols[1]]],
        y = mean_table[[mean_cols[2]]],
        error_x = error$x, error_y = error$y,
        marker = list(color = ~color, symbol = "square", size = 10)
      );
    }
  }

  return(p);
}

#' Print ENA plot
#'
#' Prints an ENA plotly object, applying any stored shapes and annotations.
#'
#' @param x An ENA plot object (class "enaplot").
#' @param ... Additional arguments (currently ignored).
#'
#' @return Invisibly returns the plotly object after printing.
#' @export
'print.enaplot' <- function(x, ...) {
  class(x) <- c("plotly", "htmlwidget");

  shapes <- attr(x = x, which = "shapes");
  annotations <- attr(x = x, which = "annotations");
  x <- plotly::layout(p = x, shapes = shapes, annotations = annotations);

  print(x);
}

##' Calculate correlations between SVD and centroid differences in ONA model
##'
##' Computes Pearson and Spearman correlations between pairwise differences of points and centroids along specified dimensions in an Epistemic Network Analysis (ENA) model.
##'
##' @param enaset An ENA set object containing points and model centroids (typically output from ENA model fitting).
##' @param pts Optional data frame of points to use. If NULL, points are selected from `enaset$points` matching the specified direction.
##' @param cts Optional data frame or matrix of centroids to use. If NULL, centroids are taken from `enaset$model$centroids`.
##' @param dims Integer vector specifying which dimensions to use for the correlation (default: c(1,2)).
##' @param direction Character string specifying which direction(s) to include from `enaset$points$ENA_DIRECTION` (default: "response").
##'
##' @return A data frame with Pearson and Spearman correlations for each specified dimension.
##' @export
correlations <- function(enaset, pts = NULL, cts = NULL, dims = c(1:2), direction = "response") {
  if(is.null(pts) || is.null(cts)) {
    # Not every ordered set carries an ENA_DIRECTION column -- a model with one
    # row per unit has no directions to separate. Check for the column by name
    # rather than reading it: filtering on an absent column would drop every
    # row, and rENA's `$.ena.points` raises an error rather than returning NULL.
    if("ENA_DIRECTION" %in% colnames(enaset$points)) {
      rows <- enaset$points$ENA_DIRECTION %in% direction
    } else {
      rows <- rep(TRUE, nrow(enaset$points))
    }
    pts <- enaset$points[rows, ]
    cts <- enaset$model$centroids[, ]
  }

  pComb = combn(nrow(pts), 2)
  point1 = pComb[1,]
  point2 = pComb[2,]

  points = as.matrix(pts)
  centroids = as.matrix(cts)
  svdDiff = matrix(points[point1, dims] - points[point2, dims], ncol=length(dims), nrow=length(point1))
  optDiff = matrix(centroids[point1, dims] - centroids[point2, dims], ncol=length(dims), nrow=length(point1))

  # svdDiff/optDiff have one column per requested dimension, so index them by
  # position rather than by the caller's dimension number. Indexing by `dims`
  # only works for the default c(1,2); any other pair reads past the end.
  corrs = as.data.frame(mapply(function(method) {
    sapply(seq_along(dims), function(dim) {
      cor(as.numeric(svdDiff[,dim]), as.numeric(optDiff[,dim]), method=method)
    });
  }, c("pearson","spearman")))

  return(corrs);
}

###
#' ona writeup
#'
#' @description provide writeup of the provided model
#'
#' @param enaset ENAset to view methods of
#' @param tool c("rENA","webENA")
#' @param tool.version as.character(packageVersion(tool))
#' @param comparison character string representing the comparison used, c(NULL, "parametric", "non-parametric"). Default NULL
#' @param comparison.groups Groups that were used for the comparison
#' @param sig.dig Integer for the number of digits to round to
#' @param output_dir Where to save the output file
#' @param type c("file","stream") File will save to a file in output_dir, Stream returns the contents directly
#' @param theory Logical indicating whether to include theory in the writeup
#' @param methods Logical indicating whether to include methods in the writeup
#' @param params additional parameters for rmarkdown::render
#' @param output_file character
#' @param output_format character
#' 
#' @importFrom utils packageVersion
#'
#' @export
#'
#' @return String representing the methods used to generate the model
ona.writeup <- function(
    enaset,
    tool = "ona", tool.version = as.character(packageVersion(tool)),
    comparison = NULL, comparison.groups = NULL, sig.dig = 2,
    output_dir = getwd(), type = c("file","stream"), theory = T, methods = T,
    params = NULL, output_file = NULL, output_format = NULL
) {
  if(is.null(enaset$`_function.params`$weight.by))
    enaset$`_function.params`$weight.by <- enaset$`_function.params`$args$weight.by

  type = match.arg(type, choices = c("file","stream"), several.ok = FALSE)

  if(type == "file") {
    output_format = "word_document"
  }

  file = rmarkdown::render(system.file("rmd","methods_ona.rmd", package="ona"), output_dir = output_dir,
                           knit_root_dir = output_dir, intermediates_dir = output_dir, quiet = TRUE,
                           params = params, output_file = output_file
                           ,output_format = output_format
                           # ,output_format = ifelse(type == "file", rENA::methods_report, rENA::methods_report_stream)
  )

  if(type == "file")
    file
  else if (type == "stream" && endsWith(file, ".plain"))
    readChar(file, file.info(file)$size)
}

#' @title methods_report
#' @description Output format for methods report in R Markdown (Word/docx)
#' @param toc Logical; include a table of contents (default: FALSE).
#' @param toc_depth Integer; depth of table of contents (default: 3).
#' @param fig_width Numeric; width (in inches) of figures (default: 5).
#' @param fig_height Numeric; height (in inches) of figures (default: 4).
#' @param keep_md Logical; keep the intermediate markdown file (default: FALSE).
#' @param md_extensions Character vector of additional markdown extensions to use (default: NULL).
#' @param pandoc_args Character vector of additional Pandoc command line arguments (default: NULL).
#'
#' @return An R Markdown output format object for use with `rmarkdown::render` (Word/docx output).
#' @export
methods_report <- function(toc = FALSE,
                           toc_depth = 3,
                           fig_width = 5,
                           fig_height = 4,
                           keep_md = FALSE,
                           md_extensions = NULL,
                           pandoc_args = NULL) {

  # knitr options and hooks
  knitr <- rmarkdown::knitr_options(
    opts_chunk = list(dev = 'png',
                      dpi = 96,
                      fig.width = fig_width,
                      fig.height = fig_height)
  )

  # build pandoc args
  args <- c("--standalone")

  # table of contents
  args <- c(args, rmarkdown::pandoc_toc_args(toc, toc_depth))

  # pandoc args
  args <- c(args, pandoc_args)

  preserved_chunks <- character()

  # pre_processor <- function(metadata, input_file, runtime, knit_meta,
  #                           files_dir, output_dir) {
  #   preserved_chunks <<- extract_preserve_chunks(input_file, knitr::extract_raw_output)
  #   NULL
  # }

  # post_processor <- function(metadata, input_file, output_file, clean, verbose) {
  #   output_str <- readLines(output_file, encoding = 'UTF-8')
  #   output_res <- knitr::restore_raw_output(output_str, preserved_chunks)
  #   if (!identical(output_str, output_res))
  #     writeLines(enc2utf8(output_res), output_file, useBytes = TRUE)
  #   output_file
  # }

  # return output format
  rmarkdown::output_format(
    knitr = knitr,
    pandoc = rmarkdown::pandoc_options(to = "docx",
                                       from = rmarkdown::from_rmarkdown(extensions = md_extensions),
                                       args = args),
    keep_md = keep_md
    # ,pre_processor = pre_processor,
    # post_processor = post_processor
  )
}

#' @title methods_report_stream
#' @description Output format for methods report in R Markdown (plain text stream)
#' @param toc Logical; include a table of contents (default: FALSE).
#' @param toc_depth Integer; depth of table of contents (default: 3).
#' @param fig_width Numeric; width (in inches) of figures (default: 5).
#' @param fig_height Numeric; height (in inches) of figures (default: 4).
#' @param keep_md Logical; keep the intermediate markdown file (default: FALSE).
#' @param md_extensions Character vector of additional markdown extensions to use (default: NULL).
#' @param pandoc_args Character vector of additional Pandoc command line arguments (default: NULL).
#'
#' @return An R Markdown output format object for use with `rmarkdown::render` (plain text output).
#' @export
methods_report_stream <- function(toc = FALSE,
                                  toc_depth = 3,
                                  fig_width = 5,
                                  fig_height = 4,
                                  keep_md = FALSE,
                                  md_extensions = NULL,
                                  pandoc_args = NULL) {

  # knitr options and hooks
  knitr <- rmarkdown::knitr_options(
    opts_chunk = list(dev = 'png',
                      dpi = 96,
                      fig.width = fig_width,
                      fig.height = fig_height)
  )

  # build pandoc args
  args <- c("--standalone")

  # table of contents
  args <- c(args, rmarkdown::pandoc_toc_args(toc, toc_depth))

  # pandoc args
  args <- c(args, pandoc_args)

  preserved_chunks <- character()

  # return output format
  rmarkdown::output_format(
    knitr = knitr,
    pandoc = rmarkdown::pandoc_options(to = "plain",
                                       from = rmarkdown::from_rmarkdown(extensions = md_extensions),
                                       args = args),
    keep_md = keep_md
  )
}

.onLoad <- function(libname, pkgname) {
  utils::globalVariables(c(
    "%chin%",".BY",".I",".SD",":=","ENA_DIRECTION","as.data.table","code","col2rgb",
    "color_adjusted","combn","data.table","distance","hsv","is","label","node_1","node_2",
    "prcomp","rgb2hsv","set","t.test","value","var","variable"
  ))
}

.onAttach <- function(libname, pkgname) {
  packageStartupMessage("For the latest features and updates, install from https://cran.qe-libs.org");
  invisible();
}

# R-level wrappers replacing the former C++ exports in src/optim.cpp.
# These preserve backward compatibility for any code calling these functions
# directly as ona::accum_stanza_window(), ona::vector_to_ut(), etc.

#' Accumulate Stanza Window (directed, full-square output)
#'
#' For each row, accumulates co-occurrences over a backward stanza window.
#' Preserves directionality: output column (i*n_codes + j) represents
#' code i (context/ground) co-occurring with code j (current row/response).
#' Returns a matrix with \code{n_codes²} columns (full square adjacency).
#'
#' Delegates to \code{\link[libqe]{accumulate_stanza}(ordered = TRUE)}.
#'
#' @param df         A data frame of code columns (one row per utterance)
#' @param windowSize Number of prior rows in the window (default 1; use
#'   \code{Inf} for all prior rows in the stanza)
#' @param binary     Logical; binarise non-zero co-occurrences (default TRUE)
#' @return Numeric matrix with \code{nrow(df)} rows and \code{ncol(df)²} columns
#' @importFrom libqe accumulate_stanza
#' @export
accum_stanza_window <- function(df, windowSize = 1, binary = TRUE) {
  INT_MAX <- .Machine$integer.max
  wb <- if (is.infinite(windowSize) || windowSize >= INT_MAX) INT_MAX
        else as.integer(windowSize)
  as.data.frame(
    libqe::accumulate_stanza(as.matrix(df), window_back = wb,
                             window_forward = 0L, binary = binary,
                             ordered = TRUE)
  )
}

#' Upper Triangle from Vector (numeric)
#'
#' Compute pairwise products v[i] * v[j] for all i < j.
#' Thin wrapper around \code{\link[rENA]{vector_to_ut}}.
#'
#' @param v Numeric vector or single-row matrix of code values
#' @return Numeric row vector of pairwise products (length choose_two(length(v)))
#' @export
vector_to_ut <- function(v) {
  rENA::vector_to_ut(v)
}

#' Least-squares node positions for directed ENA
#'
#' Thin wrapper around \code{\link[rENA]{directed_node_positions}}.
#'
#' @param line_weights Numeric matrix (units x connections)
#' @param points       Numeric matrix of rotated points (units x dims)
#' @param numDims      Number of dimensions
#' @return List with nodes, centroids, weights, points
#' @export
directed_node_positions <- function(line_weights, points, numDims) {
  rENA::directed_node_positions(line_weights, points, numDims)
}

#' Directed node positions with ground+response rows combined
#'
#' Thin wrapper around
#' \code{\link[rENA]{directed_node_positions_with_ground_response_added}}.
#'
#' @param line_weights Numeric matrix (units x connections)
#' @param points       Numeric matrix of rotated points (units x dims)
#' @param numDims      Number of dimensions
#' @return List with nodes, centroids, weights, points
#' @export
directed_node_positions_with_ground_response_added <- function(line_weights,
                                                                points,
                                                                numDims) {
  rENA::directed_node_positions_with_ground_response_added(line_weights,
                                                           points, numDims)
}

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
}

#' Normalize and optimize a directed (ordered) network accumulation
#'
#' Thin wrapper around \code{\link[rENA]{model}} that applies the ONA-specific
#' defaults: SVD rotation, zero-network exclusion from centering, and
#' center-to-origin post-processing.
#'
#' @param x     An \code{ena.ordered.set} produced by
#'   \code{\link[rENA]{accumulate}(ordered = TRUE)} or
#'   \code{tma::accumulate_contexts}.
#' @param ...   Additional arguments forwarded to \code{rENA::model}.
#' @param rotate.using Character. \code{"svd"} (default) for standard principal
#'   components rotation; \code{"mean"} for means rotation.
#' @param rotation.params Additional parameters for the rotation function
#'   (passed as \code{rotate_params} to \code{rENA::model}).
#' @param rotation.set Optional pre-computed \code{ENARotationSet}. When
#'   supplied the rotation step is skipped and this set is used directly.
#' @param center_to_origin Logical. When \code{TRUE} (default) the projected
#'   space is translated so the mean of all points is at the origin.
#'
#' @return An \code{ena.ordered.set} with calculated points and node positions.
#' @export
model <- function(
  x, ...,
  rotate.using     = "svd",
  rotation.params  = NULL,
  rotation.set     = NULL,
  center_to_origin = TRUE
) {
  rotate_fun <- if (rotate.using == "mean") {
    rENA::ena.rotate.by.mean
  } else {
    rENA::ena.rotate.by.svd
  }

  rotate_params <- if (!is.null(rotation.params)) rotation.params else list()

  if (!is.null(rotation.set)) {
    # Custom rotation: normalise + centre, then apply the supplied rotation set
    # directly (skipping the rotation step) before projecting and optimising.
    x <- rENA::sphere_norm(x)
    x <- rENA::center(x, exclude_zero_networks = TRUE)

    if (!is(rotation.set, "ena.rotation.set")) {
      stop("Supplied rotation.set is not an instance of ENARotationSet")
    }
    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

    x <- rENA::project(x)
    x <- rENA::optimize(x)
  } else {
    x <- rENA::model(
      x, ...,
      rotate_fun           = rotate_fun,
      rotate_params        = rotate_params,
      exclude_zero_networks = TRUE,
      center_to_origin     = center_to_origin
    )
  }

  return(x)
}

#' 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))
}

###
#' ONA rotation by mean difference
#'
#' Computes a dimensional reduction (rotation) such that the first dimension
#' passes through the means of two groups. Subsequent dimensions use SVD on
#' the deflated data. Thin wrapper around \code{\link[rENA]{ena.rotate.by.mean}};
#' kept for backward compatibility with existing ONA scripts.
#'
#' @param x      A matrix or \code{ena.points} data table of ENA points
#'   (typically centered and normalized). Metadata columns are stripped
#'   automatically via \code{as.matrix}.
#' @param groups A list of length 2 (or more), each element being a length-2
#'   list \code{list(a, b)} where \code{a} and \code{b} are logical vectors or
#'   character unit IDs identifying the two groups.
#' @param codes  Character vector of code names corresponding to columns of
#'   \code{x}; used to label the rows of the returned rotation matrix.
#'
#' @return A list with \code{rotation}, \code{codes}, \code{eigenvalues}, and
#'   \code{node.positions = NULL}.
#' @export
###
ona.rotate.by.mean <- function(x, groups, codes) {
  # Construct a minimal enaset-compatible list so we can delegate to
  # rENA::ena.rotate.by.mean which now calls libqe::means_rotation internally.
  data_matrix <- as.matrix(x)   # as.matrix.ena.matrix strips metadata columns
  enaset <- list(
    model                  = list(points.for.projection = data_matrix),
    points.normed.centered = NULL,
    connection.counts      = list(ENA_UNIT = x$ENA_UNIT),
    line.weights           = data_matrix,
    rotation               = list(codes = codes)
  )
  rENA::ena.rotate.by.mean(enaset, groups = groups)
}

1		##' Calculate correlations between SVD and centroid differences in ONA model
2		##'
3		##' Computes Pearson and Spearman correlations between pairwise differences of points and centroids along specified dimensions in an Epistemic Network Analysis (ENA) model.
4		##'
5		##' @param enaset An ENA set object containing points and model centroids (typically output from ENA model fitting).
6		##' @param pts Optional data frame of points to use. If NULL, points are selected from `enaset$points` matching the specified direction.
7		##' @param cts Optional data frame or matrix of centroids to use. If NULL, centroids are taken from `enaset$model$centroids`.
8		##' @param dims Integer vector specifying which dimensions to use for the correlation (default: c(1,2)).
9		##' @param direction Character string specifying which direction(s) to include from `enaset$points$ENA_DIRECTION` (default: "response").
10		##'
11		##' @return A data frame with Pearson and Spearman correlations for each specified dimension.
12		##' @export
13		correlations <- function(enaset, pts = NULL, cts = NULL, dims = c(1:2), direction = "response") {
14	!	if(is.null(pts) \|\| is.null(cts)) {
15		# Not every ordered set carries an ENA_DIRECTION column -- a model with one
16		# row per unit has no directions to separate. Check for the column by name
17		# rather than reading it: filtering on an absent column would drop every
18		# row, and rENA's `$.ena.points` raises an error rather than returning NULL.
19	!	if("ENA_DIRECTION" %in% colnames(enaset$points)) {
20	!	rows <- enaset$points$ENA_DIRECTION %in% direction
21		} else {
22	!	rows <- rep(TRUE, nrow(enaset$points))
23		}
24	!	pts <- enaset$points[rows, ]
25	!	cts <- enaset$model$centroids[, ]
26		}
27
28	!	pComb = combn(nrow(pts), 2)
29	!	point1 = pComb[1,]
30	!	point2 = pComb[2,]
31
32	!	points = as.matrix(pts)
33	!	centroids = as.matrix(cts)
34	!	svdDiff = matrix(points[point1, dims] - points[point2, dims], ncol=length(dims), nrow=length(point1))
35	!	optDiff = matrix(centroids[point1, dims] - centroids[point2, dims], ncol=length(dims), nrow=length(point1))
36
37		# svdDiff/optDiff have one column per requested dimension, so index them by
38		# position rather than by the caller's dimension number. Indexing by `dims`
39		# only works for the default c(1,2); any other pair reads past the end.
40	!	corrs = as.data.frame(mapply(function(method) {
41	!	sapply(seq_along(dims), function(dim) {
42	!	cor(as.numeric(svdDiff[,dim]), as.numeric(optDiff[,dim]), method=method)
43		});
44	!	}, c("pearson","spearman")))
45
46	!	return(corrs);
47		}

1		###
2		#' ona writeup
3		#'
4		#' @description provide writeup of the provided model
5		#'
6		#' @param enaset ENAset to view methods of
7		#' @param tool c("rENA","webENA")
8		#' @param tool.version as.character(packageVersion(tool))
9		#' @param comparison character string representing the comparison used, c(NULL, "parametric", "non-parametric"). Default NULL
10		#' @param comparison.groups Groups that were used for the comparison
11		#' @param sig.dig Integer for the number of digits to round to
12		#' @param output_dir Where to save the output file
13		#' @param type c("file","stream") File will save to a file in output_dir, Stream returns the contents directly
14		#' @param theory Logical indicating whether to include theory in the writeup
15		#' @param methods Logical indicating whether to include methods in the writeup
16		#' @param params additional parameters for rmarkdown::render
17		#' @param output_file character
18		#' @param output_format character
19		#'
20		#' @importFrom utils packageVersion
21		#'
22		#' @export
23		#'
24		#' @return String representing the methods used to generate the model
25		ona.writeup <- function(
26		enaset,
27		tool = "ona", tool.version = as.character(packageVersion(tool)),
28		comparison = NULL, comparison.groups = NULL, sig.dig = 2,
29		output_dir = getwd(), type = c("file","stream"), theory = T, methods = T,
30		params = NULL, output_file = NULL, output_format = NULL
31		) {
32	!	if(is.null(enaset$`_function.params`$weight.by))
33	!	enaset$`_function.params`$weight.by <- enaset$`_function.params`$args$weight.by
34
35	!	type = match.arg(type, choices = c("file","stream"), several.ok = FALSE)
36
37	!	if(type == "file") {
38	!	output_format = "word_document"
39		}
40
41	!	file = rmarkdown::render(system.file("rmd","methods_ona.rmd", package="ona"), output_dir = output_dir,
42	!	knit_root_dir = output_dir, intermediates_dir = output_dir, quiet = TRUE,
43	!	params = params, output_file = output_file
44	!	,output_format = output_format
45		# ,output_format = ifelse(type == "file", rENA::methods_report, rENA::methods_report_stream)
46		)
47
48	!	if(type == "file")
49	!	file
50	!	else if (type == "stream" && endsWith(file, ".plain"))
51	!	readChar(file, file.info(file)$size)
52		}
53
54		#' @title methods_report
55		#' @description Output format for methods report in R Markdown (Word/docx)
56		#' @param toc Logical; include a table of contents (default: FALSE).
57		#' @param toc_depth Integer; depth of table of contents (default: 3).
58		#' @param fig_width Numeric; width (in inches) of figures (default: 5).
59		#' @param fig_height Numeric; height (in inches) of figures (default: 4).
60		#' @param keep_md Logical; keep the intermediate markdown file (default: FALSE).
61		#' @param md_extensions Character vector of additional markdown extensions to use (default: NULL).
62		#' @param pandoc_args Character vector of additional Pandoc command line arguments (default: NULL).
63		#'
64		#' @return An R Markdown output format object for use with `rmarkdown::render` (Word/docx output).
65		#' @export
66		methods_report <- function(toc = FALSE,
67		toc_depth = 3,
68		fig_width = 5,
69		fig_height = 4,
70		keep_md = FALSE,
71		md_extensions = NULL,
72		pandoc_args = NULL) {
73
74		# knitr options and hooks
75	!	knitr <- rmarkdown::knitr_options(
76	!	opts_chunk = list(dev = 'png',
77	!	dpi = 96,
78	!	fig.width = fig_width,
79	!	fig.height = fig_height)
80		)
81
82		# build pandoc args
83	!	args <- c("--standalone")
84
85		# table of contents
86	!	args <- c(args, rmarkdown::pandoc_toc_args(toc, toc_depth))
87
88		# pandoc args
89	!	args <- c(args, pandoc_args)
90
91	!	preserved_chunks <- character()
92
93		# pre_processor <- function(metadata, input_file, runtime, knit_meta,
94		# files_dir, output_dir) {
95		# preserved_chunks <<- extract_preserve_chunks(input_file, knitr::extract_raw_output)
96		# NULL
97		# }
98
99		# post_processor <- function(metadata, input_file, output_file, clean, verbose) {
100		# output_str <- readLines(output_file, encoding = 'UTF-8')
101		# output_res <- knitr::restore_raw_output(output_str, preserved_chunks)
102		# if (!identical(output_str, output_res))
103		# writeLines(enc2utf8(output_res), output_file, useBytes = TRUE)
104		# output_file
105		# }
106
107		# return output format
108	!	rmarkdown::output_format(
109	!	knitr = knitr,
110	!	pandoc = rmarkdown::pandoc_options(to = "docx",
111	!	from = rmarkdown::from_rmarkdown(extensions = md_extensions),
112	!	args = args),
113	!	keep_md = keep_md
114		# ,pre_processor = pre_processor,
115		# post_processor = post_processor
116		)
117		}
118
119		#' @title methods_report_stream
120		#' @description Output format for methods report in R Markdown (plain text stream)
121		#' @param toc Logical; include a table of contents (default: FALSE).
122		#' @param toc_depth Integer; depth of table of contents (default: 3).
123		#' @param fig_width Numeric; width (in inches) of figures (default: 5).
124		#' @param fig_height Numeric; height (in inches) of figures (default: 4).
125		#' @param keep_md Logical; keep the intermediate markdown file (default: FALSE).
126		#' @param md_extensions Character vector of additional markdown extensions to use (default: NULL).
127		#' @param pandoc_args Character vector of additional Pandoc command line arguments (default: NULL).
128		#'
129		#' @return An R Markdown output format object for use with `rmarkdown::render` (plain text output).
130		#' @export
131		methods_report_stream <- function(toc = FALSE,
132		toc_depth = 3,
133		fig_width = 5,
134		fig_height = 4,
135		keep_md = FALSE,
136		md_extensions = NULL,
137		pandoc_args = NULL) {
138
139		# knitr options and hooks
140	!	knitr <- rmarkdown::knitr_options(
141	!	opts_chunk = list(dev = 'png',
142	!	dpi = 96,
143	!	fig.width = fig_width,
144	!	fig.height = fig_height)
145		)
146
147		# build pandoc args
148	!	args <- c("--standalone")
149
150		# table of contents
151	!	args <- c(args, rmarkdown::pandoc_toc_args(toc, toc_depth))
152
153		# pandoc args
154	!	args <- c(args, pandoc_args)
155
156	!	preserved_chunks <- character()
157
158		# return output format
159	!	rmarkdown::output_format(
160	!	knitr = knitr,
161	!	pandoc = rmarkdown::pandoc_options(to = "plain",
162	!	from = rmarkdown::from_rmarkdown(extensions = md_extensions),
163	!	args = args),
164	!	keep_md = keep_md
165		)
166		}

1		.onLoad <- function(libname, pkgname) {
2	!	utils::globalVariables(c(
3	!	"%chin%",".BY",".I",".SD",":=","ENA_DIRECTION","as.data.table","code","col2rgb",
4	!	"color_adjusted","combn","data.table","distance","hsv","is","label","node_1","node_2",
5	!	"prcomp","rgb2hsv","set","t.test","value","var","variable"
6		))
7		}
8
9		.onAttach <- function(libname, pkgname) {
10	2x	packageStartupMessage("For the latest features and updates, install from https://cran.qe-libs.org");
11	2x	invisible();
12		}

1		# R-level wrappers replacing the former C++ exports in src/optim.cpp.
2		# These preserve backward compatibility for any code calling these functions
3		# directly as ona::accum_stanza_window(), ona::vector_to_ut(), etc.
4
5		#' Accumulate Stanza Window (directed, full-square output)
6		#'
7		#' For each row, accumulates co-occurrences over a backward stanza window.
8		#' Preserves directionality: output column (i*n_codes + j) represents
9		#' code i (context/ground) co-occurring with code j (current row/response).
10		#' Returns a matrix with \code{n_codes²} columns (full square adjacency).
11		#'
12		#' Delegates to \code{\link[libqe]{accumulate_stanza}(ordered = TRUE)}.
13		#'
14		#' @param df A data frame of code columns (one row per utterance)
15		#' @param windowSize Number of prior rows in the window (default 1; use
16		#' \code{Inf} for all prior rows in the stanza)
17		#' @param binary Logical; binarise non-zero co-occurrences (default TRUE)
18		#' @return Numeric matrix with \code{nrow(df)} rows and \code{ncol(df)²} columns
19		#' @importFrom libqe accumulate_stanza
20		#' @export
21		accum_stanza_window <- function(df, windowSize = 1, binary = TRUE) {
22	!	INT_MAX <- .Machine$integer.max
23	!	wb <- if (is.infinite(windowSize) \|\| windowSize >= INT_MAX) INT_MAX
24	!	else as.integer(windowSize)
25	!	as.data.frame(
26	!	libqe::accumulate_stanza(as.matrix(df), window_back = wb,
27	!	window_forward = 0L, binary = binary,
28	!	ordered = TRUE)
29		)
30		}
31
32		#' Upper Triangle from Vector (numeric)
33		#'
34		#' Compute pairwise products v[i] * v[j] for all i < j.
35		#' Thin wrapper around \code{\link[rENA]{vector_to_ut}}.
36		#'
37		#' @param v Numeric vector or single-row matrix of code values
38		#' @return Numeric row vector of pairwise products (length choose_two(length(v)))
39		#' @export
40		vector_to_ut <- function(v) {
41	!	rENA::vector_to_ut(v)
42		}
43
44		#' Least-squares node positions for directed ENA
45		#'
46		#' Thin wrapper around \code{\link[rENA]{directed_node_positions}}.
47		#'
48		#' @param line_weights Numeric matrix (units x connections)
49		#' @param points Numeric matrix of rotated points (units x dims)
50		#' @param numDims Number of dimensions
51		#' @return List with nodes, centroids, weights, points
52		#' @export
53		directed_node_positions <- function(line_weights, points, numDims) {
54	!	rENA::directed_node_positions(line_weights, points, numDims)
55		}
56
57		#' Directed node positions with ground+response rows combined
58		#'
59		#' Thin wrapper around
60		#' \code{\link[rENA]{directed_node_positions_with_ground_response_added}}.
61		#'
62		#' @param line_weights Numeric matrix (units x connections)
63		#' @param points Numeric matrix of rotated points (units x dims)
64		#' @param numDims Number of dimensions
65		#' @return List with nodes, centroids, weights, points
66		#' @export
67		directed_node_positions_with_ground_response_added <- function(line_weights,
68		points,
69		numDims) {
70	!	rENA::directed_node_positions_with_ground_response_added(line_weights,
71	!	points, numDims)
72		}

1		to_square <- function(x) {
2	!	n <- sqrt(length(x))
3
4	!	dim(x) <- c(n, n);
5	!	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	!	if(is.factor(x)) {
16	!	x = as.character(x)
17		}
18	!	class(x) = c("ena.co.occurrence", class(x))
19	!	x
20		}

1		#' Normalize and optimize a directed (ordered) network accumulation
2		#'
3		#' Thin wrapper around \code{\link[rENA]{model}} that applies the ONA-specific
4		#' defaults: SVD rotation, zero-network exclusion from centering, and
5		#' center-to-origin post-processing.
6		#'
7		#' @param x An \code{ena.ordered.set} produced by
8		#' \code{\link[rENA]{accumulate}(ordered = TRUE)} or
9		#' \code{tma::accumulate_contexts}.
10		#' @param ... Additional arguments forwarded to \code{rENA::model}.
11		#' @param rotate.using Character. \code{"svd"} (default) for standard principal
12		#' components rotation; \code{"mean"} for means rotation.
13		#' @param rotation.params Additional parameters for the rotation function
14		#' (passed as \code{rotate_params} to \code{rENA::model}).
15		#' @param rotation.set Optional pre-computed \code{ENARotationSet}. When
16		#' supplied the rotation step is skipped and this set is used directly.
17		#' @param center_to_origin Logical. When \code{TRUE} (default) the projected
18		#' space is translated so the mean of all points is at the origin.
19		#'
20		#' @return An \code{ena.ordered.set} with calculated points and node positions.
21		#' @export
22		model <- function(
23		x, ...,
24		rotate.using = "svd",
25		rotation.params = NULL,
26		rotation.set = NULL,
27		center_to_origin = TRUE
28		) {
29	1x	rotate_fun <- if (rotate.using == "mean") {
30	!	rENA::ena.rotate.by.mean
31		} else {
32	1x	rENA::ena.rotate.by.svd
33		}
34
35	1x	rotate_params <- if (!is.null(rotation.params)) rotation.params else list()
36
37	1x	if (!is.null(rotation.set)) {
38		# Custom rotation: normalise + centre, then apply the supplied rotation set
39		# directly (skipping the rotation step) before projecting and optimising.
40	!	x <- rENA::sphere_norm(x)
41	!	x <- rENA::center(x, exclude_zero_networks = TRUE)
42
43	!	if (!is(rotation.set, "ena.rotation.set")) {
44	!	stop("Supplied rotation.set is not an instance of ENARotationSet")
45		}
46	!	x$rotation.matrix <- rotation.set$rotation.matrix
47	!	x$rotation$rotation.matrix <- rotation.set$rotation.matrix
48	!	x$rotation$nodes <- rotation.set$nodes
49	!	x$rotation$eigenvalues <- rotation.set$eigenvalues
50
51	!	x <- rENA::project(x)
52	!	x <- rENA::optimize(x)
53		} else {
54	1x	x <- rENA::model(
55	1x	x, ...,
56	1x	rotate_fun = rotate_fun,
57	1x	rotate_params = rotate_params,
58	1x	exclude_zero_networks = TRUE,
59	1x	center_to_origin = center_to_origin
60		)
61		}
62
63	1x	return(x)
64		}

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	!	points = as.matrix(set$points)[filter_rows, ];
13	!	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	!	if (length(unique(set$points[filter_rows,]$ENA_DIRECTION)) == 2) {
17	!	positions = directed_node_positions_with_ground_response_added(weights, points, ncol(points));
18		}
19		else {
20	!	positions = directed_node_positions(weights, points, ncol(points));
21		}
22
23	!	node.positions = positions$nodes;
24	!	rownames(node.positions) = set$enadata$codes;
25
26	!	return(list("node.positions" = node.positions, "centroids" = positions$centroids))
27		}

1		###
2		#' ONA rotation by mean difference
3		#'
4		#' Computes a dimensional reduction (rotation) such that the first dimension
5		#' passes through the means of two groups. Subsequent dimensions use SVD on
6		#' the deflated data. Thin wrapper around \code{\link[rENA]{ena.rotate.by.mean}};
7		#' kept for backward compatibility with existing ONA scripts.
8		#'
9		#' @param x A matrix or \code{ena.points} data table of ENA points
10		#' (typically centered and normalized). Metadata columns are stripped
11		#' automatically via \code{as.matrix}.
12		#' @param groups A list of length 2 (or more), each element being a length-2
13		#' list \code{list(a, b)} where \code{a} and \code{b} are logical vectors or
14		#' character unit IDs identifying the two groups.
15		#' @param codes Character vector of code names corresponding to columns of
16		#' \code{x}; used to label the rows of the returned rotation matrix.
17		#'
18		#' @return A list with \code{rotation}, \code{codes}, \code{eigenvalues}, and
19		#' \code{node.positions = NULL}.
20		#' @export
21		###
22		ona.rotate.by.mean <- function(x, groups, codes) {
23		# Construct a minimal enaset-compatible list so we can delegate to
24		# rENA::ena.rotate.by.mean which now calls libqe::means_rotation internally.
25	!	data_matrix <- as.matrix(x) # as.matrix.ena.matrix strips metadata columns
26	!	enaset <- list(
27	!	model = list(points.for.projection = data_matrix),
28	!	points.normed.centered = NULL,
29	!	connection.counts = list(ENA_UNIT = x$ENA_UNIT),
30	!	line.weights = data_matrix,
31	!	rotation = list(codes = codes)
32		)
33	!	rENA::ena.rotate.by.mean(enaset, groups = groups)
34		}