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