1"""
   2- Authors: Peter Mawhorter, Tiffany Lin, & Presha Goel
   3- Consulted:
   4- Date: 2022-4-15
   5- Purpose: Code to support visualizing decision graphs and explorations.
   6
   7Defines functions for graph layout and drawing for
   8`exploration.core.DecisionGraph` objects. See the `explorationViewer`
   9module for more info on how these are used. This module computes layout
  10positions, but actually displaying the graphs is done via HTML.
  11
  12TODO: Anchor-free localization implementation?
  13"""
  14
  15from typing import (
  16    Dict, Tuple, Literal, TypeAlias, Sequence, Optional, Set, Union,
  17    List, cast
  18)
  19
  20import math
  21import copy
  22
  23import networkx as nx
  24
  25from . import base
  26from . import core
  27from . import analysis
  28
  29BlockPosition: 'TypeAlias' = Tuple[int, int, int]
  30"""
  31A type alias: block positions indicate the x/y coordinates of the
  32north-west corner of a node, as well as its side length in grid units
  33(all nodes are assumed to be square).
  34"""
  35
  36
  37BlockLayout: 'TypeAlias' = Dict[base.DecisionID, BlockPosition]
  38"""
  39A type alias: block layouts map each decision in a particular graph to a
  40block position which indicates both position and size in a unit grid.
  41"""
  42
  43
  44def roomSize(connections: int) -> int:
  45    """
  46    For a room with the given number of connections, returns the side
  47    length of the smallest square which can accommodate that many
  48    connections. Note that outgoing/incoming reciprocal pairs to/from a
  49    single destination should only count as one connection, because they
  50    don't need more than one space on the room perimeter. Even with zero
  51    connections, we still return 1 as the room size.
  52    """
  53    if connections == 0:
  54        return 1
  55    return 1 + (connections - 1) // 4
  56
  57
  58def expandBlocks(layout: BlockLayout) -> None:
  59    """
  60    Modifies the given block layout by adding extra space between each
  61    positioned node: it triples the coordinates of each node, and then
  62    shifts them south and east by 1 unit each, by maintaining the nodes
  63    at their original sizes, TODO...
  64    """
  65    # TODO
  66
  67
  68#def blockLayoutFor(region: core.DecisionGraph) -> BlockLayout:
  69#    """
  70#    Computes a unit-grid position and size for each room in an
  71#    `exploration.core.DecisionGraph`, laying out the rooms as
  72#    non-overlapping square blocks. In many cases, connections will be
  73#    stretched across empty space, but no explicit space is reserved for
  74#    connections.
  75#    """
  76#    # TODO
  77
  78GraphLayoutMethod: 'TypeAlias' = Literal[
  79    "stacked",
  80    "square",
  81    "line",
  82    "arc"
  83]
  84"""
  85The options for layouts of a decision graph. They are:
  86
  87- 'stacked': Assigns *all* nodes to position (0, 0). Use this if you want
  88    to generate an empty layout that you plan to modify. Doesn't require
  89    any attributes.
  90- 'square': Takes the square root of the number of decisions, then places
  91    them in order into a square with that side length (rounded up). This
  92    is a very simple but also terrible algorithm. Doesn't require any
  93    attributes.
  94- 'line': Lays out the decisions in a straight line. Doesn't require any
  95    attributes.
  96- 'arc': Lays out decisions on an arc. Doesn't require any attributes.
  97"""
  98
  99def assignPositions(
 100    decisions: Sequence[base.DecisionID],
 101    attributes: Optional[Dict[base.DecisionID, dict]] = None,
 102    method: GraphLayoutMethod = "square"
 103) -> base.Layout:
 104    """
 105    Given a sequence of decision IDs, plus optionally a dictionary
 106    mapping those IDs to attribute dictionaries, computes a layout for
 107    the decisions according to the specified method, returning a
 108    dictionary mapping each decision ID to its position in the layout.
 109
 110    Different layout methods may required different attributes to be
 111    available.
 112    """
 113    if method == "stacked":
 114        return {d: (0, 0) for d in decisions}  #  all nodes at (0, 0)
 115    elif method == "square":
 116        return assignSquarePositions(decisions)
 117    elif method == "line":
 118        return assignLinePositions(decisions)
 119    elif method == "arc":
 120        return assignArcPositions(decisions)
 121    else:
 122        raise ValueError(f"Invalid layout method {method!r}.")
 123
 124
 125def assignSquarePositions(
 126    decisions: Sequence[base.DecisionID]
 127) -> base.Layout:
 128    """
 129    Creates and returns a dictionary of positions for the given sequence
 130    of decisions using the 'square' layout: it arranges them into a big
 131    square.
 132    """
 133    result = {}
 134    # Figure out side length of the square that will fit them all
 135    side = math.ceil((len(decisions)**0.5))
 136    # Put 'em in a square
 137    for i, d in enumerate(decisions):
 138        result[d] = (i % side, i // side)
 139    return result
 140
 141
 142def assignLinePositions(
 143    decisions: Sequence[base.DecisionID]
 144) -> base.Layout:
 145    """
 146    Creates and returns a dictionary of positions for the given sequence
 147    of decisions using the 'line' layout: it arranges them into a
 148    straight horizontal line.
 149    """
 150    result = {}
 151    # Put 'em in a line
 152    for i, d in enumerate(decisions):
 153        result[d] = (float(i), 0.0)
 154    return result
 155
 156
 157def assignArcPositions(
 158    decisions: Sequence[base.DecisionID],
 159    radiusFactor: float = 1.5
 160) -> base.Layout:
 161    """
 162    Creates and returns a dictionary of positions for the given sequence
 163    of decisions using the 'arc' layout: it arranges them along an arc
 164    curving towards -y relative to a straight line. The start and end
 165    of the arc are at (0, 0) and (N, 0) respectively, where N is the
 166    number of decisions to be laid out. The `radiusFactor` must be at
 167    least 1; when it's 1 the radius is N/2 and the shape is a full
 168    semi-circle, as `radiusFactor` increases the radius of the circle
 169    increases proportionally making the arc gentler.
 170    """
 171    n = len(decisions)
 172    cx = n / 2  # x-position of center halfway down line
 173    r = radiusFactor * cx  # radius of the circle
 174    # y-position so that that radius intersects both (0, 0) and (n, 0)
 175    cy = (r * r - cx * cx)**0.5
 176    # angle from straight-down to either end of arc (note y and x
 177    # intentionally reversed here)
 178    th = math.atan2(cx, cy)
 179    thSlice = th * 2 / n
 180    thStart = 1.5 * math.pi - th
 181    result = {}
 182    # Traverse the arc
 183    for i, d in enumerate(decisions):
 184        angle = thStart + thSlice * i
 185        result[d] = (
 186            cx + r * math.cos(angle),
 187            cy + r * math.sin(angle)
 188        )
 189    return result
 190
 191
 192def setFinalPositions(
 193    exploration: core.DiscreteExploration,
 194    method: GraphLayoutMethod = "square"
 195) -> None:
 196    """
 197    Adds a "final" layout to the given exploration which contains a
 198    dictionary mapping decision IDs to `Position`s. Every decision that
 199    ever existed over the course of the exploration is assigned a
 200    position.
 201    """
 202    exploration.layouts["final"] = assignPositions(
 203        exploration.allDecisions(),
 204        method=method
 205    )
 206
 207
 208def setPathPositions(
 209    exploration: core.DiscreteExploration,
 210    arc: bool = True
 211) -> None:
 212    """
 213    Adds a "path" layout to the given exploration which contains a
 214    dictionary mapping decision IDs to `Position`s. This includes every
 215    visited decision and all of their neighbors, but does NOT include
 216    decisions which were never visited and are not neighbors of a visited
 217    decision.
 218
 219    Positions nodes along an arc below the X axis, each slightly more
 220    than 1 unit apart, with unvisited neighbors descending below the node
 221    they're a neighbor of in a parabola. This layout attempts to set
 222    things up so that edges don't end up parallel and overlapping. You
 223    can set `arc` to False to instead lay things out on a line, with
 224    unvisited neighbors descending in a line below the node they connect
 225    to. The linear layout would have lots of ambiguities if you tried to
 226    relax it using a force simulation, whereas the arc layout would
 227    have few ambiguities.
 228    """
 229    # Lay out visited decisions in a line:
 230    onPath = exploration.allVisitedDecisions()
 231    if arc:
 232        result = assignArcPositions(onPath, 5)  # a gentle arc
 233    else:
 234        result = assignLinePositions(onPath)  # a straight line
 235    # Get the final graph to add neighbors from
 236    finalGraph = exploration.getSituation().graph
 237    # Track already-accounted-for neighbors
 238    seen: Set[base.DecisionID] = set()
 239    # Add neighbors to our layout
 240    for decision in onPath:
 241        # copy x of this decision
 242        x = result[decision][0]
 243        dx = 0.05
 244        # Track y coordinates going down below the arc
 245        y = result[decision][1] - 1.0
 246        try:
 247            neighbors = finalGraph.destinationsFrom(decision)
 248        except core.MissingDecisionError:
 249            continue
 250            # decision on path may have been merged or deleted by end of
 251            # exploration. We don't want to get neighbors in the step the
 252            # node was visited, because it's likely many of those will
 253            # have been explored by the time we get to the final graph,
 254            # and we also don't want to include any merged/deleted
 255            # neighbors in our graph, despite the fact we're including
 256            # merged or deleted path nodes.
 257        # TODO: Sort these by step added?
 258        for dID in neighbors.values():
 259            # We only include all neighbors which are not elsewhere on
 260            # the path. It's possible some of these may be confirmed, but
 261            # that's fine, they were not active at any step.
 262            if dID not in onPath:
 263                result[dID] = (x, y)
 264                # not just straight down if we're arcing
 265                if arc:
 266                    dx += 0.03  # dx changes so nodes aren't in a line
 267                    x += dx  # next node will be slightly farther to the right
 268                y -= 1.0  # next node will be lower down
 269
 270    exploration.layouts["path"] = result
 271
 272
 273def setCustomPositions(
 274    exploration: core.DiscreteExploration,
 275    fromLayout: str = "final"
 276) -> None:
 277    """
 278    Sets up a 'custom' layout to be edited. Copies from the "final"
 279    layout by default, or you can specify which layout to copy from. If
 280    the specified layout doesn't exist, a square layout of all decisions
 281    is used.
 282    """
 283    copyFrom = exploration.layouts.get(fromLayout)
 284    if copyFrom is None:
 285        copied = assignPositions(
 286            exploration.allDecisions(),
 287            method="square"
 288        )
 289    else:
 290        copied = copy.deepcopy(copyFrom)
 291    exploration.layouts["custom"] = copied
 292
 293
 294#---------------#
 295# Baryeccentric #
 296#---------------#
 297
 298Number: 'TypeAlias' = Union[int, float]
 299"""
 300For arguments that can be either an integer or a float.
 301"""
 302
 303
 304# Testing support functions
 305def pz(n: Number) -> Number:
 306    """
 307    Converts -0.0 to 0.0 and leaves all other numbers alone.
 308    """
 309    if n == -0.0:
 310        return 0.0
 311    else:
 312        return n
 313
 314
 315def rt(t: base.LayoutPosition) -> base.LayoutPosition:
 316    """
 317    Rounds off both parts of a 2-element tuple to 6 decimal places, and
 318    then also converts -0.0 to +0.0 in either position.
 319    """
 320    return (
 321        pz(round(t[0], 6)),
 322        pz(round(t[1], 6))
 323    )
 324
 325
 326def rtl(tl: Sequence[base.LayoutPosition]) -> List[base.LayoutPosition]:
 327    """
 328    Applies `rt` to a sequence of positions, returning a list.
 329    """
 330    return [rt(t) for t in tl]
 331
 332
 333def distance(a: base.LayoutPosition, b: base.LayoutPosition) -> float:
 334    """
 335    Calculates the distance between two points, using the distance
 336    formula in 2 dimensions. For example:
 337
 338    >>> distance((0, 0), (0, 3))
 339    3.0
 340    >>> distance((0, 0), (3, 0))
 341    3.0
 342    >>> distance((0, 0), (3, 4))
 343    5.0
 344    """
 345    x1, y1 = a
 346    x2, y2 = b
 347
 348    return (((x2 - x1) ** 2) + ((y2 - y1) ** 2)) ** 0.5
 349
 350
 351def mid(
 352    a: base.LayoutPosition,
 353    b: base.LayoutPosition
 354) -> base.LayoutPosition:
 355    """
 356    Returns the midpoint between two points. For example:
 357
 358    >>> rt(mid((0, 0), (1, 0)))
 359    (0.5, 0.0)
 360    >>> rt(mid((0, 0), (3, 8)))
 361    (1.5, 4.0)
 362    >>> rt(mid((3, -3), (-3, 3)))
 363    (0.0, 0.0)
 364    """
 365    x1, y1 = a
 366    x2, y2 = b
 367
 368    return ((x1 + x2) / 2, (y1 + y2) / 2)
 369
 370
 371def vAdd(
 372    a: base.LayoutPosition,
 373    b: base.LayoutPosition
 374) -> base.LayoutPosition:
 375    """
 376    Returns the vector addition result for two layout positions.
 377    For example:
 378
 379    >>> vAdd((0.0, 0.0), (1.0, 1.0))
 380    (1.0, 1.0)
 381    >>> vAdd((1.0, 1.0), (0.0, 0.0))
 382    (1.0, 1.0)
 383    >>> vAdd((1.0, 1.0), (2.0, -3.0))
 384    (3.0, -2.0)
 385    >>> vAdd((1.0, 1.0), (1.0, 1.0))
 386    (2.0, 2.0)
 387    """
 388    x1, y1 = a
 389    x2, y2 = b
 390
 391    return (x1 + x2, y1 + y2)
 392
 393
 394def vSub(
 395    a: base.LayoutPosition,
 396    b: base.LayoutPosition
 397) -> base.LayoutPosition:
 398    """
 399    Returns the vector between points a and b (from b to a, which is
 400    also a - b in vector math). For example:
 401
 402    >>> vSub((1.0, 1.0), (0.0, 0.0))
 403    (1.0, 1.0)
 404    >>> vSub((2.0, -3.0), (1.0, 1.0))
 405    (1.0, -4.0)
 406    >>> vSub((1.0, 1.0), (1.0, 1.0))
 407    (0.0, 0.0)
 408    """
 409    x1, y1 = a
 410    x2, y2 = b
 411
 412    return (x1 - x2, y1 - y2)
 413
 414def norm(v: base.LayoutPosition) -> base.LayoutPosition:
 415    """
 416    Normalizes the given vector, returning a vector in the same direction
 417    whose length is 1. Returns the zero-vector if given the zero-vector,
 418    which is the only case where the length of the result is not 1.
 419
 420    For example:
 421
 422    >>> norm((0, 0))
 423    (0.0, 0.0)
 424    >>> norm((2, 0))
 425    (1.0, 0.0)
 426    >>> norm((102, 0))
 427    (1.0, 0.0)
 428    >>> norm((0, -3))
 429    (0.0, -1.0)
 430    """
 431    length = distance((0, 0), v)
 432    if length == 0:
 433        return (0.0, 0.0)
 434    else:
 435        return (v[0] / length, v[1] / length)
 436
 437
 438def isATriangle(x: Number, y: Number, z: Number) -> bool:
 439    """
 440    Checks whether three side lengths can form a triangle. For example:
 441
 442    >>> isATriangle(3, 4, 5)
 443    True
 444    >>> isATriangle(1, 1, 1)
 445    True
 446    >>> isATriangle(100, 99, 1)
 447    False
 448    >>> isATriangle(100, 99, 2)
 449    True
 450    >>> isATriangle(2, 100, 99)
 451    True
 452    >>> isATriangle(99, 2, 100)
 453    True
 454    >>> isATriangle(3, 2, 10)
 455    False
 456    >>> isATriangle(5, 1, 1)
 457    False
 458    >>> isATriangle(9, 18.01, 9)
 459    False
 460    """
 461    return ((x + y > z) and (y + z > x) and (z + x > y))
 462
 463
 464def scaleBy(
 465    vector: base.LayoutPosition,
 466    scale: Number
 467) -> base.LayoutPosition:
 468    """
 469    Scales the given vector by the specified scale.
 470    Examples:
 471
 472    >>> rt(scaleBy((1.0, 0.0), 3))
 473    (3.0, 0.0)
 474    >>> rt(scaleBy((3.0, 4.0), 10))
 475    (30.0, 40.0)
 476    >>> rt(scaleBy((6.0, 8.0), 5))
 477    (30.0, 40.0)
 478    >>> rt(scaleBy((0.0, 2.0), -2))
 479    (0.0, -4.0)
 480    >>> rt(scaleBy((0.0, 0.0), 1000))
 481    (0.0, 0.0)
 482    """
 483    x, y = vector
 484    return (x * scale, y * scale)
 485
 486
 487def scaleTo(
 488    vector: base.LayoutPosition,
 489    length: Number
 490) -> base.LayoutPosition:
 491    """
 492    Scales the given vector to the specified length. Note that if the
 493    vector is (0, 0), it will remain (0, 0) so the result won't actually
 494    have the specified length in that one case. Examples:
 495
 496    >>> rt(scaleTo((1, 0), 3))
 497    (3.0, 0.0)
 498    >>> rt(scaleTo((3, 4), 10))
 499    (6.0, 8.0)
 500    >>> rt(scaleTo((6, 8), 5))
 501    (3.0, 4.0)
 502    >>> rt(scaleTo((0, 2), -2))
 503    (0.0, -2.0)
 504    >>> rt(scaleTo((0, 0), 1000))
 505    (0.0, 0.0)
 506    """
 507    lengthNow = distance((0, 0), vector)
 508    if lengthNow == 0:
 509        return (0.0, 0.0)
 510    else:
 511        x, y = vector
 512        return (
 513            (x / lengthNow) * length,
 514            (y / lengthNow) * length
 515        )
 516
 517
 518def circleIntersections(
 519    a: base.LayoutPosition,
 520    b: base.LayoutPosition,
 521    aRadius: Number,
 522    bRadius: Number
 523) -> List[base.LayoutPosition]:
 524    """
 525    Calculates the intersection point(s) between two circles centered at
 526    points `a` and `b` with the given radii. Returns a list of 0, 1, or 2
 527    positions depending on the relationship between the circles. Note
 528    that if two circles are the same circle, it should in theory return a
 529    list of infinite positions; in that case we return a list with four
 530    positions that are the places where horizontal and vertical lines
 531    through the shared center intersect the shared circle. Examples:
 532
 533    >>> rtl(circleIntersections((0, 0), (2, 0), 1, 1))  # single point
 534    [(1.0, 0.0)]
 535    >>> rtl(circleIntersections((0, 0), (0, 2), 1, 1))  # single point
 536    [(0.0, 1.0)]
 537    >>> rtl(circleIntersections((0, 0), (6, 8), 5, 5))  # two 3/4/5 triangles
 538    [(3.0, 4.0)]
 539    >>> rtl(circleIntersections((0, 0), (2, 0), 1.5, 1.5))  # two points
 540    [(1.0, -1.118034), (1.0, 1.118034)]
 541    >>> rtl(circleIntersections((0, 0), (0, 0), 2, 3))  # no points
 542    []
 543    >>> rtl(circleIntersections((0, 0), (2, 0), 0.5, 0.5))  # no points
 544    []
 545    >>> rtl(circleIntersections((-3, 4), (3, 4), 5, 5))  # two 3/4/5 triangles
 546    [(0.0, 0.0), (0.0, 8.0)]
 547    >>> rtl(circleIntersections((-4, -3), (4, -3), 5, 5))  # two 3/4/5 triangles
 548    [(0.0, -6.0), (0.0, 0.0)]
 549    >>> rtl(circleIntersections((0.0, 0.0), (0.0, 0.0), 5, 5))  # infinity
 550    [(5.0, 0.0), (0.0, -5.0), (-5.0, 0.0), (0.0, 5.0)]
 551    """
 552    x1, y1 = a
 553    x2, y2 = b
 554
 555    d = distance(a, b)
 556
 557    # If the circles are too far apart or if the distance between their
 558    # centers is so small compared to the difference between their radii
 559    # that one is entirely inside the other, there are no points of
 560    # intersection.
 561    if d > aRadius + bRadius or d < abs(aRadius - bRadius):
 562        return []
 563    elif aRadius + bRadius == d:  # one point of intersection
 564        vec = scaleTo((x2 - x1, y2 - y1), aRadius)
 565        return [ (x1 + vec[0], y1 + vec[1]) ]
 566    elif x1 == x2 and y1 == y2 and aRadius == bRadius:  # same circle
 567        return [  # clockwise from 3 o'clock if +y is up and +x is right
 568            (x1 + aRadius, y1),
 569            (x1, y1 - aRadius),
 570            (x1 - aRadius, y1),
 571            (x1, y1 + aRadius)
 572        ]
 573
 574    # Otherwise we have 2 points of intersection
 575    # TODO: Explain this math a bit...
 576    aa = (aRadius**2 - bRadius**2 + d**2) / (2 * d)
 577    h = (aRadius**2 - aa**2)**0.5
 578
 579    x0 = x1 + aa * (x2 - x1) / d
 580    y0 = y1 + aa * (y2 - y1) / d
 581
 582    x3 = x0 + h * (y2 - y1) / d
 583    y3 = y0 - h * (x2 - x1) / d
 584
 585    x4 = x0 - h * (y2 - y1) / d
 586    y4 = y0 + h * (x2 - x1) / d
 587
 588    return [(x3, y3), (x4, y4)]
 589
 590
 591def bestFitIntersection(
 592    a: base.LayoutPosition,
 593    b: base.LayoutPosition,
 594    aRad: Number,
 595    bRad: Number
 596):
 597    """
 598    Given two circles which may or may not intersect (specified as
 599    centers `a` and `b` plus respective radii), returns a point that's 
 600    on the line through their centers that's on the shortest segment of
 601    that line which connects one circle to the other (this may or may
 602    not be between the two centers if one circle encircles the other).
 603
 604    The point is placed along that segment so that its distance from
 605    circle a divided by its distance from circle b is proportional to
 606    the radius of circle a divided by the radius of circle b (it ends up
 607    closer to the smaller circle).
 608
 609    If the two circles have the same center, we return a point that has
 610    the same y-coordinate as that center, and of the two equally valid
 611    points of that nature, we return the one with the greater
 612    x-coordinate.
 613
 614    Some examples:
 615    >>> rt(bestFitIntersection((0, 0), (100, 0), 180, 40))
 616    (147.272727, 0.0)
 617    >>> rt(bestFitIntersection((0, 0), (50, 87), 60, 78))
 618    (21.73913, 37.826087)
 619    >>> rt(bestFitIntersection((100, 0), (50, 87), 70, 78))
 620    (76.351351, 41.148649)
 621
 622    >>> rt(bestFitIntersection((0, 0), (8, 6), 5, 5))  # circles touch
 623    (4.0, 3.0)
 624    >>> rt(bestFitIntersection((0, 0), (12, 9), 10, 5))  # circles touch
 625    (8.0, 6.0)
 626    >>> rt(bestFitIntersection((-20, -20), (-30, 20), 10, 10))  # r1 == r2
 627    (-25.0, 0.0)
 628    >>> rt(bestFitIntersection((-30, 20), (-20, -20), 10, 10))  # other order
 629    (-25.0, 0.0)
 630    >>> rt(bestFitIntersection((0, 0), (0, 0), 12, 24))  # same center
 631    (16.0, 0.0)
 632    >>> # we arbitrarily pick a point horizontal from the center
 633    >>> # note that (-16.0, 0.0) is equally valid but we pick the option
 634    >>> # with the higher x-coordinate
 635    >>> rt(bestFitIntersection((0, 0), (0, 0), 24, 12))  # works same other way
 636    (16.0, 0.0)
 637    >>> rt(bestFitIntersection((0, 0), (0, 0), 10, 10))  # same circle
 638    (10.0, 0.0)
 639    >>> rt(bestFitIntersection((0, 0), (0, 0), 0, 0))  # zero-radius same center
 640    (0.0, 0.0)
 641    >>> rt(bestFitIntersection((0, 0), (2, 0), 0, 0))  # zero-radius diff center
 642    (1.0, 0.0)
 643    >>> rt(bestFitIntersection((2, 0), (0, 0), 0, 0))  # other direction
 644    (1.0, 0.0)
 645    >>> rt(bestFitIntersection((0, 0), (2, 0), 1, 0))  # single zero-radius
 646    (2.0, 0.0)
 647    """
 648    import sys
 649    dist = distance(a, b)
 650    vect = norm(vSub(b, a))
 651    if vect == (0.0, 0.0):  # same center
 652        vect = (1.0, 0.0)  # horizontal
 653
 654    # Find all four points at which the line between a and b intersects
 655    # the a and b circles:
 656    aVec = scaleBy(vect, aRad)
 657    aLow = vSub(a, aVec)
 658    aHigh = vAdd(a, aVec)
 659    bVec = scaleBy(vect, bRad)
 660    bLow = vSub(b, bVec)
 661    bHigh = vAdd(b, bVec)
 662
 663    # Now find which pair of one point from A and one from B is closest
 664    # There's probably a more mathy way to do this, but this is simple to
 665    # reason about.
 666    closest = None
 667    bestDist = None
 668    for (p1, p2) in [
 669        (aHigh, bHigh),
 670        (aHigh, bLow),
 671        (aLow, bHigh),
 672        (aLow, bLow),
 673    ]:
 674        pointSep = distance(p1, p2)
 675        # Note strict < here biases towards earlier pairs in the order
 676        # above, such that 'high' points beat low ones on ties
 677        if closest is None or pointSep < bestDist:
 678            closest = (p1, p2)
 679            bestDist = pointSep
 680
 681    assert closest is not None
 682
 683    # Now find a point between the two closest points-on-circle where the
 684    # proportion between distances to each matches the ratio of the radii
 685    # of the circles:
 686    onA, onB = closest
 687    between = vSub(onB, onA)
 688    if between == (0.0, 0.0):  # same point, so return it
 689        return onA
 690    dirBetween = norm(between)
 691    distBetween = distance(onA, onB)
 692    if aRad + bRad == 0:  # both zero-radius; return average of the two
 693        return ((onA[0] + onB[0]) / 2, (onA[1] + onB[1]) / 2)
 694    howFarAlong = aRad / (aRad + bRad)
 695    return vAdd(onA, scaleBy(dirBetween, howFarAlong * distBetween))
 696
 697
 698def baryeccentricPosition(
 699    a: base.LayoutPosition,
 700    b: base.LayoutPosition,
 701    c: base.LayoutPosition,
 702    distA: Number,
 703    distB: Number,
 704    distC: Number
 705):
 706    """
 707    Returns a "baryeccentric" position given three reference points and
 708    three numbers indicating distances to each of them. If the distances
 709    are in agreement and together specify a particular point within (or
 710    outside of) the reference triangle, we return that point. If the
 711    two or more of the distances are too short to touch each other, we
 712    compromise at a position most consistent with them, and if the
 713    distances are too long we also compromise.
 714
 715    For best results, you should ensure that the reference points make a
 716    triangle rather than a line or point.
 717
 718    We find a compromise by treating each reference point + distance as a
 719    circle. We first compute the intersection points between each pair of
 720    circles, resulting in 0-4 intersection points per pair (see
 721    `circleIntersection`). For pairs with no intersection, we use
 722    `bestFitIntersection` to come up with a single "intersection" point.
 723    Now for pairs with 2+ intersection points, we pick the single
 724    intersection point whose distance to the third point is most
 725    consistent with the measured third distance. This leaves us with 3
 726    intersection points: one for each pair of reference points. We
 727    average these three points to come up with the final result.
 728
 729    TODO: consider the perfectly-overlapping circles case a bit more...
 730
 731    Some examples:
 732
 733    >>> baryeccentricPosition((0, 0), (6, 8), (6, 0), 5, 5, 5)
 734    (3.0, 4.0)
 735    >>> baryeccentricPosition((0, 0), (-6, 8), (-6, 0), 5, 5, 5)
 736    (-3.0, 4.0)
 737    >>> baryeccentricPosition((0, 0), (-6, -8), (-6, 0), 5, 5, 5)
 738    (-3.0, -4.0)
 739    >>> baryeccentricPosition((0, 0), (3.0, 4.0), (3.0, 0), 5, 0, 4)
 740    (3.0, 4.0)
 741    >>> baryeccentricPosition((0, 0), (3.0, 4.0), (3.0, 0), 0, 5, 3)
 742    (0.0, 0.0)
 743    >>> baryeccentricPosition((0, 0), (3.0, 4.0), (3.0, 0), 3, 4, 0)
 744    (3.0, 0.0)
 745    >>> rt(baryeccentricPosition((-8, 6), (8, 6), (0, -10), 10, 10, 10))
 746    (0.0, 0.0)
 747    >>> rt(baryeccentricPosition((-8, 6), (8, 6), (0, -12), 10, 10, 0))
 748    (0.0, -8.0)
 749    >>> rt(baryeccentricPosition((-8, -6), (0, 12), (8, -6), 10, 0, 10))
 750    (0.0, 8.0)
 751    >>> rt(baryeccentricPosition((0, 12), (-8, -6), (8, -6), 0, 10, 10))
 752    (0.0, 8.0)
 753    >>> rt(baryeccentricPosition((-4, 3), (4, 3), (0, -5), 5, 5, 0))
 754    (0.0, -3.333333)
 755    >>> rt(baryeccentricPosition((-1, 0), (1, 0), (0, -1), 1, 1, 1))
 756    (0.0, 0.0)
 757    >>> rt(baryeccentricPosition(
 758    ...     (-25.3, 45.8), (12.4, -24.3), (35.9, 58.2),
 759    ...     61.2, 35.5, 28.4
 760    ... ))
 761    (27.693092, 20.240286)
 762    >>> rt(baryeccentricPosition(
 763    ...     (-25.3, 45.8), (12.4, -24.3), (35.9, 58.2),
 764    ...     102.5, 12.8, 89.4
 765    ... ))
 766    (28.437607, -32.62218)
 767
 768    Edge case examples:
 769
 770    >>> baryeccentricPosition((0, 0), (0, 0), (0, 0), 5, 5, 5)
 771    (5.0, 0.0)
 772    >>> baryeccentricPosition((0, 0), (0, 0), (0, 0), 0, 0, 0)
 773    (0.0, 0.0)
 774    """
 775    # TODO: Should we print a warning if the points aren't a triangle?
 776
 777    # First, find intersection point(s) for each pair
 778    abPoints = circleIntersections(a, b, distA, distB)
 779    acPoints = circleIntersections(a, c, distA, distC)
 780    bcPoints = circleIntersections(b, c, distB, distC)
 781
 782    # if circles don't touch, add an estimated point
 783    if len(abPoints) == 0:
 784        abPoints = [bestFitIntersection(a, b, distA, distB)]
 785    if len(acPoints) == 0:
 786        acPoints = [bestFitIntersection(a, c, distA, distC)]
 787    if len(bcPoints) == 0:
 788        bcPoints = [bestFitIntersection(b, c, distB, distC)]
 789
 790    # If circles touch a multiple places, narrow that down to one by
 791    # figuring out which is most consistent with the third distance
 792    if len(abPoints) == 1:
 793        abPoint = abPoints[0]
 794    else:  # must be > 1 point per above
 795        assert len(abPoints) > 1
 796        abPoint = None
 797        bestError = None
 798        for p in abPoints:
 799            thirdDist = distance(p, c)
 800            error = abs(thirdDist - distC)
 801            if abPoint is None or error < cast(float, bestError):
 802                abPoint = p
 803                bestError = error
 804
 805    if len(acPoints) == 1:
 806        acPoint = acPoints[0]
 807    else:  # must be > 1 point per above
 808        assert len(acPoints) > 1
 809        acPoint = None
 810        bestError = None
 811        for p in acPoints:
 812            thirdDist = distance(p, b)
 813            error = abs(thirdDist - distB)
 814            if bestError is None or error < bestError:
 815                acPoint = p
 816                bestError = error
 817
 818    if len(bcPoints) == 1:
 819        bcPoint = bcPoints[0]
 820    else:  # must be > 1 point per above
 821        assert len(bcPoints) > 1
 822        bcPoint = None
 823        bestError = None
 824        for p in bcPoints:
 825            thirdDist = distance(p, a)
 826            error = abs(thirdDist - distA)
 827            if bestError is None or error < bestError:
 828                bcPoint = p
 829                bestError = error
 830
 831    assert abPoint is not None
 832    assert acPoint is not None
 833    assert bcPoint is not None
 834
 835    # At this point, ab/ac/bc point variables should be assigned properly
 836    return (
 837        (abPoint[0] + acPoint[0] + bcPoint[0]) / 3,
 838        (abPoint[1] + acPoint[1] + bcPoint[1]) / 3,
 839    )
 840
 841
 842def baryeccentricLayout(
 843    exploration: core.DiscreteExploration,
 844    specifiedNodes: Optional[base.Layout] = None
 845) -> base.Layout:
 846    """
 847    Computes a baryeccentric coordinate layout for all decisions in the
 848    final step of the given exploration, using the specified positions
 849    of a few nodes given in `specifiedNodes`. `specifiedNodes` should
 850    specify positions for at least 3 decisions, and those positions must
 851    form a triangle, not a line or point. If `specifiedNodes` does not
 852    contain enough decisions (or if it's not provided), decisions will
 853    be added to it as follows:
 854
 855    - If it's empty, add the node with the lowest id at position (0, 0).
 856    - If it's got only one decision or we just added one node, add the
 857        node that's furthest from that node in terms of hop distance.
 858        We'll position this second node at the same y-coordinate as the
 859        first, but with an x-coordinate equal to the hop distance between
 860        it and the first node. If multiple nodes are tied for furthest,
 861        add the one with the lowest id.
 862    - If it's got only two decisions or we just added one or two, add the
 863        node whose sum of hop distances to the two already selected is
 864        largest. We position this third node such that the hop distances
 865        to each of the already-placed nodes are respected and it forms a
 866        triangle, or if that's not possible due to those distances being
 867        too short, we position it partway between them proportional to
 868        those two distances with an artificial offset perpendicular to
 869        the line between the two other points. Ties are broken towards
 870        nodes with a shorter max hop distance to either of the two
 871        already-placed nodes, and then towards lower node IDs.
 872
 873    If the number of nodes in the entire graph is 1 or 2, we return a
 874    layout positioning the first node at (0, 0) and (if it exists) the
 875    second node at (1, 0).
 876
 877    Some examples:
 878
 879    # TODO
 880    >> baryeccentricLayout(TODO)
 881    """
 882    hops = analysis.shortestHopPaths(
 883        exploration[-1].graph,
 884        lambda src, transition, dst, graph: (
 885            'journey' not in graph.transitionTags(src, transition)
 886        )
 887    )
 888    # Now we can use `analysis.hopDistnace` given `hops` plus two
 889    # decision IDs to get the hop distance between any two decisions.
 890
 891    # Create empty layout by default:
 892    if specifiedNodes is None:
 893        specifiedNodes = {}
 894
 895    # Select at least 3 specific nodes
 896    if len(specifiedNodes) < 3:
 897        finalGraph = exploration[-1].graph
 898        allDecisions = sorted(finalGraph)
 899
 900        # Bail out if we have fewer than 3 total decisions
 901        if len(allDecisions) < 3:
 902            result = {}
 903            if len(allDecisions) > 0:
 904                result[allDecisions[0]] = (0.0, 0.0)
 905                if len(allDecisions) > 1:
 906                    result[allDecisions[1]] = (1.0, 0.0)
 907            return result
 908
 909        # Add a decision at (0, 0) if we didn't have any specified
 910        if len(specifiedNodes) < 1:
 911            # Find largest weakly connected component:
 912            bigCC = max(nx.weakly_connected_components(finalGraph), key=len)
 913            # Use an arbitrary node from that component
 914            specifiedNodes[list(bigCC)[0]] = (0.0, 0.0)
 915
 916        assert len(specifiedNodes) >= 1
 917
 918        # If 1 specified or just added, add furthest-away decision
 919        if len(specifiedNodes) < 2:
 920            first = list(specifiedNodes)[0]
 921            best = None
 922            bestDist = None
 923            # Find furthest connected node
 924            for dID in allDecisions:
 925                if dID == first:
 926                    # Skip node that we already assigned
 927                    continue
 928                dist = analysis.hopDistance(hops, dID, first)
 929                # Note > here breaks ties towards lower IDs
 930                if dist is not None and (bestDist is None or dist > bestDist):
 931                    best = dID
 932                    bestDist = dist
 933            # if no nodes are connected, we've got a big problem, but
 934            # we'll push on by selecting the node with the second-lowest
 935            # node ID.
 936            if best is None:
 937                # Find first un-specified node ID:
 938                second = None
 939                for second in allDecisions:
 940                    dFirst = analysis.hopDistance(hops, second, first)
 941                    if (
 942                         second not in specifiedNodes
 943                     and dFirst is not None
 944                    ):
 945                        # second will remain at this value after loop
 946                        break
 947                else:  # if we never hit break
 948                    for second in allDecisions:
 949                        if second not in specifiedNodes:
 950                            # second will remain at this value after loop
 951                            break
 952                assert second is not None
 953                # Just put it at (1, 0) since hops aren't informative
 954                specifiedNodes[second] = (1.0, 0.0)
 955            else:
 956                assert best != first
 957                assert bestDist is not None
 958                firstPos = specifiedNodes[first]
 959                # Same y-value as first one, with x-dist as hop dist
 960                specifiedNodes[best] = (firstPos[0] + bestDist, firstPos[1])
 961
 962        assert len(specifiedNodes) >= 2
 963
 964        # If only two specified (and or one or two just added) we look
 965        # for the node with best combined distance to those two,
 966        # breaking ties towards smaller max-distance to either and then
 967        # towards smaller ID values.
 968        if len(specifiedNodes) < 3:
 969            first, second = list(specifiedNodes)[:2]
 970            best = None
 971            bestCombined = None
 972            bestLonger = None
 973            bestDists = None
 974            for dID in allDecisions:
 975                if dID in specifiedNodes:
 976                    # Skip already-placed nodes
 977                    continue
 978                distA = analysis.hopDistance(hops, dID, first)
 979                distB = analysis.hopDistance(hops, dID, second)
 980                if distA is None or distB is None:
 981                    # Note: *shouldn't* be possible for only one to be
 982                    # None, but we don't take chances here
 983                    continue
 984                combined = distA + distB
 985                longer = max(distA, distB)
 986                if (
 987                    # first one
 988                    bestCombined is None
 989                    # better combined distance (further away)
 990                 or combined > bestCombined
 991                    # tied combined and better max distance (more evenly
 992                    # placed between at *shorter* max dist)
 993                 or (
 994                        combined == bestCombined
 995                    and longer < cast(float, bestLonger)
 996                        # Note strict < here breaks ties towards lower IDs
 997                    )
 998                ):
 999                    best = dID
1000                    bestCombined = combined
1001                    bestLonger = longer
1002                    bestDists = (distA, distB)
1003
1004            firstPos = specifiedNodes[first]
1005            secondPos = specifiedNodes[second]
1006
1007            abDist = analysis.hopDistance(hops, first, second)
1008           # They were chosen based on being connected...
1009            assert abDist is not None
1010            assert bestDists is not None
1011            # Could happen if only two nodes are connected, for example
1012            if best is None or sum(bestDists) < abDist:
1013                # Just put it artificially between them
1014                vect = (
1015                    secondPos[0] - firstPos[0],
1016                    secondPos[1] - firstPos[1]
1017                )
1018                # perpendicular vector
1019                ortho = (vect[1], -vect[0])
1020                # Just use first decision that's not already specified
1021                if best is not None:
1022                    third = best
1023                else:
1024                    third = None
1025                    for third in allDecisions:
1026                        thirdHopsA = analysis.hopDistance(hops, first, third)
1027                        thirdHopsB = analysis.hopDistance(hops, second, third)
1028                        if (
1029                             third not in specifiedNodes
1030                         and thirdHopsA is not None
1031                         and thirdHopsB is not None
1032                        ):
1033                            # third will remain on this node
1034                            break
1035                    else:  # if we never hit the break
1036                        for third in allDecisions:
1037                            if third not in specifiedNodes:
1038                                # third will remain on this node
1039                                break
1040                    assert third is not None
1041                assert third != first
1042                assert third != second
1043                # Offset orthogonally by half the distance between
1044                specifiedNodes[third] = (
1045                    firstPos[0] + vect[0]/2 + ortho[0]/2,
1046                    firstPos[1] + vect[1]/2 + ortho[0]/2
1047                )
1048            else:
1049                # Position the best candidate to form a triangle where
1050                # distances are proportional; we know distances are long
1051                # enough to make a triangle
1052                distA = analysis.hopDistance(hops, dID, first)
1053                candidates = circleIntersections(
1054                    firstPos,
1055                    secondPos,
1056                    *bestDists
1057                )
1058                if len(candidates) == 0:
1059                    assert distA is not None
1060                    assert distB is not None
1061                    where = bestFitIntersection(
1062                        firstPos,
1063                        secondPos,
1064                        distA,
1065                        distB
1066                    )
1067                else:
1068                    where = candidates[0]
1069                assert best != first
1070                assert best != second
1071                specifiedNodes[best] = where
1072
1073            assert len(specifiedNodes) >= 3
1074
1075    # TODO: Don't just use first 3 here...
1076    # Grab first 3 decision IDs from layout
1077    a, b, c = list(specifiedNodes.keys())[:3]
1078    # Get their positions
1079    aPos = specifiedNodes[a]
1080    bPos = specifiedNodes[b]
1081    cPos = specifiedNodes[c]
1082    # create initial result using just specified positions
1083    result = {
1084        a: aPos,
1085        b: bPos,
1086        c: cPos
1087    }
1088    # Now we need to compute positions of each other node...
1089    # We use `exploration.allDecisions` as the set of nodes we want to
1090    # establish positions for, even though some of them may have been
1091    # deleted by the end and thus may not appear in our hops data.
1092    toLayOut = exploration.allDecisions()
1093    # value for default positions
1094    default = 1.0
1095    for decision in toLayOut:
1096        aHops = analysis.hopDistance(hops, a, decision)
1097        bHops = analysis.hopDistance(hops, b, decision)
1098        cHops = analysis.hopDistance(hops, c, decision)
1099
1100        # if hops is none for one, it should be none for all
1101        if aHops is None or bHops is None or cHops is None:
1102            # Put it at a default position on a parabola
1103            # TODO: Better default here?
1104            result[decision] = (default, default**1.1)
1105            default += 0.1
1106        else:
1107            assert aHops is not None
1108            assert bHops is not None
1109            assert cHops is not None
1110
1111            # Place according to baryeccentric position
1112            result[decision] = baryeccentricPosition(
1113                aPos,
1114                bPos,
1115                cPos,
1116                aHops,
1117                bHops,
1118                cHops
1119            )
1120
1121    # Return result at end...
1122    return result
1123
1124def setBaryeccentricPositions(
1125    exploration: core.DiscreteExploration,
1126    method: GraphLayoutMethod = "square"
1127) -> None:
1128    """
1129    Adds a "baryeccentric" layout to the given exploration that uses the
1130    `baryeccentricLayout` function to determine node positions. Uses
1131    an empty set of specified nodes so that they'll be determined
1132    automatically.
1133    """
1134    exploration.layouts["baryeccentric"] = baryeccentricLayout(
1135        exploration,
1136        {}
1137    )