<jupyter_start><jupyter_text>Hex grid===================*<NAME>**November 12, 2015*<jupyter_code>apothem = 1.0 # inner radius
offset = 0.0 # set to 0.0 for flat-topped ("rhex"), 0.5 for pointy-topped (standard KENO "hex")
num_sides = 6
offset = math.fmod(num_sides * 3 + 4 * offset, 4.0) / 4.0
hextype = "flattop" if offset == 0.5 else "pointytop"
# Face normals and displacement vectors
normals = np.empty((num_sides, 2))
displacements = np.empty((num_sides, 2))
for n in range(num_sides):
    theta = 2 * math.pi * (n + offset) / num_sides
    normals[n] = [math.cos(theta), math.sin(theta)]
    displacements[n] = normals[n] * apothem
normals[:3, :]<jupyter_output><empty_output><jupyter_text>With hexes, the unit basis vectors are along the normals through the zeroth and first faces. It's easier to scale the basis vectors by twice the apothem (distance from the center to a face), so that the distance between adjacent hex center points is always 1 in the hex coordinate system:$$ \vec{u} = 2a \cos \frac{\pi \delta}{3} \vec{i} + 2a \sin \frac{\pi \delta}{3} \vec{j} $$and $$ \vec{v} = 2a \cos \frac{\pi (1 + \delta)}{3} \vec{i} + 2a \sin \frac{\pi (1 + \delta)}{3} \vec{j} $$Here, $\delta$ is the "net" offset (different from the user-input offset): $\tfrac{1}{2}$ for flat-topped hexes (in KENO, "rhexagon"), and 0 for pointy-topped hexes (in KENO, "hexagon").<jupyter_code>bases = 2 * displacements[0:2,:]
bases
outer_radius = apothem / math.cos(math.pi / 6)
poly_orient = math.pi * offset
1/math.cos(math.pi / 6)
fig, ax = plt.subplots()

# Plot hex
patch = RegularPolygon((0., 0.), num_sides, radius=outer_radius, orientation=poly_orient,
                       facecolor='none', edgecolor='k')
ax.set_xlim(-2.2 * apothem, 2.2 * apothem)
ax.set_ylim(-1.5 * apothem, 2.5 * apothem)
ax.set_aspect('equal')
ax.add_patch(patch)
ax.plot(0, 0, 'ko')

# Plot faces
ax.plot(displacements[:,0], displacements[:,1], 'ro')
for n in range(num_sides):
    ax.text(1.25 * displacements[n,0], 1.25 * displacements[n,1], str(n),
           ha='center', va='center')
    
# Plot bases
for (i, lab) in enumerate("uv"):
    basis = bases[i,:]
    ax.arrow(0, 0, basis[0], basis[1], head_width=.1 * apothem, head_length=.15 * apothem, fc='k')
    ax.text(1.1 * basis[0], 1.1 * basis[1], "${}$".format(lab))
# Additional basis function
if 1:
    basis = 2 * displacements[2,:]
    lab = 'w'
    c = (.75,.75,.75)
    ax.arrow(0, 0, basis[0], basis[1], head_width=.1 * apothem, head_length=.15 * apothem, fc=c, color=c)
    ax.text(1.1 * basis[0], 1.1 * basis[1], "${}$".format(lab), color=c)
    
ax.grid()
fig.savefig("../doc/hex-{}-bases.pdf".format(hextype))<jupyter_output><empty_output><jupyter_text>(Note correspondence between faces of adjacent hexes: `(n + 3) % 6`. We will need to check surface senses, however, to determine whether they're orthogonal or not (because external senses will change). Erf.) The basis functions allow us to convert from the hexagonal coordinate system $u \vec{u} + v \vec{v}$ to the Cartesian coordinate system $x \vec{i} + y \vec{j}$. We simply evaluate the basis vectors given the rotation offset $\delta$.For flat-topped hexes, the "hex coordinates" $(u,v)$ relate to the cartesian coordinates $(x, y)$ by the transform$$  \begin{bmatrix}     x \\     y  \end{bmatrix}    =   2a  \begin{bmatrix}    \tfrac{\sqrt{3}}{2} & \tfrac{1}{2} \\    0 & 1   \end{bmatrix}    \begin{bmatrix}     u \\     v   \end{bmatrix} $$For pointy-topped hexes ($\delta = \tfrac{1}{2}$), the transformation is $$  \begin{bmatrix}     x \\     y  \end{bmatrix}    =   2a  \begin{bmatrix}    1 & 0 \\    \tfrac{1}{2} & \tfrac{\sqrt{3}}{2}   \end{bmatrix}    \begin{bmatrix}     u \\     v   \end{bmatrix} $$<jupyter_code>def retick(ax, dx=2):
    (lo, hi) = ax.get_xlim()
    ax.set_xticks(np.arange(math.ceil(lo), math.floor(hi + dx), dx))
    (lo, hi) = ax.get_ylim()
    ax.set_yticks(np.arange(math.ceil(lo), math.floor(hi + dx), dx))
fig, ax = plt.subplots()

patches = []
(umax, vmax) = (3, 5)
for (u, v) in itertools.product(range(umax), range(vmax)):
    uv = np.array([u, v])
    xy = uv.dot(bases) # multiply by the column vector (u,v)
    patches.append(RegularPolygon(xy, num_sides, outer_radius, poly_orient))
    ax.text(xy[0], xy[1], "$({},{})$".format(u, v), ha='center', va='center')

ax.set_xlim(-1.2 * outer_radius, xy[0] + 1.2 * outer_radius)
ax.set_ylim(-1.2 * outer_radius, xy[1] + 1.2 * outer_radius)
ax.set_aspect('equal')
retick(ax)
ax.add_collection(PatchCollection(patches, facecolor='none', edgecolor='k'))
ax.grid()
fig.savefig("../doc/hex-{}-coords-{:d}x{:d}.pdf".format(hextype, umax, vmax))<jupyter_output><empty_output><jupyter_text>Convert from Cartesian to $(u,v)$ To convert from Cartesian to grid coordinates, we invert the basis functions: $$  \begin{bmatrix}     u \\     v  \end{bmatrix}    =   B^{-1}   \begin{bmatrix}     x \\     y   \end{bmatrix}    $$<jupyter_code>inv_bases = np.linalg.inv(bases)
inv_bases
def xy_to_uv(x, y):
    xy = np.array([x,y])
    return xy.dot(inv_bases)<jupyter_output><empty_output><jupyter_text>This gives us decimal coordinates in the $(u,v)$ basis system<jupyter_code>print(xy_to_uv(2.0, 0))
print(xy_to_uv(1.0, 2.0))
print(xy_to_uv(8.0, 10.0))
x = np.linspace(-2 * apothem, 2 * apothem, 65)
y = x
(X, Y) = np.meshgrid(x, y)
XY  = np.stack([X,Y])
UV = np.empty((2, y.size, x.size))
for (j, i) in itertools.product(range(y.size), range(x.size)):
    UV[:,j,i] = XY[:,j,i].dot(inv_bases)
U = UV[0,:,:]
V = UV[1,:,:]
norm = Normalize(vmin=-1, vmax=1.)
cmap = 'rb_linear'

def plot_hex_field(ax, val):
    pc = ax.pcolormesh(X, Y, val, norm=norm, cmap=cmap)
    patches = []
    points = []
    for (u, v) in itertools.product(range(-1,2), range(-1,2)):
        uv = np.array([u, v])
        xy = uv.dot(bases) # multiply by the column vector (u,v)
        patches.append(RegularPolygon(xy, num_sides, outer_radius, poly_orient))
        patch = RegularPolygon((0., 0.), num_sides, radius=outer_radius, orientation=poly_orient,
                           facecolor='none', edgecolor='k')
        points.append(xy)
    ax.add_collection(PatchCollection(patches, facecolor='none', edgecolor='k'))
    points = np.array(points)
    ax.plot(points[:,0], points[:,1], 'ko')
    ax.set_aspect('equal')
    ax.set_xlim(x[0], x[-1])
    ax.set_ylim(y[0], y[-1])
    retick(ax, 1.0)
    return pc
(fig, (axu, axv)) = plt.subplots(ncols=2, figsize=(8, 3))
pcu = plot_hex_field(axu, U)
axu.set_title('$u$ value')
pcv = plot_hex_field(axv, V)
axv.set_title('$v$ value')
fig.colorbar(pcu)<jupyter_output><empty_output><jupyter_text>These pseudocolor plots show that our $(x, y) \to (u, v)$ conversion correctly project the hex axes; isovalues are along the $\vec u$ and $\vec v$ normals.  Third basis functionA solution to rounding these correctly to the integer hex coordinates is to rephrase the coordinate system to add a third "unit vector", drawn in the earlier graph as the grayed-out $w$ along face index 2. It's a simple linear combination of the other two axes, defined so that $u + v + w = 0$:<jupyter_code>W = -U - V
fig, ax = plt.subplots()
pc = plot_hex_field(ax, W)
axv.set_title('$w$ value')
fig.colorbar(pc)
max_uvw = np.maximum.reduce([U,V,W,-U,-V,-W,])

fig, ax = plt.subplots()
pc = plot_hex_field(ax, max_uvw)
ax.contour(X, Y, max_uvw, [0.5, 1.0], colors='g')
fig.colorbar(pc);<jupyter_output><empty_output><jupyter_text>Next, we round the $(u,v,w)$ coordinates to the nearest whole number:<jupyter_code>UVW = np.stack([U, V, W])

rounded = np.round(UVW)
fig, axes = plt.subplots(ncols=3, figsize=(12, 2.5))
for (i, ax) in enumerate(axes):
    pc = plot_hex_field(ax, rounded[i,:,:])
fig.colorbar(pc)<jupyter_output><empty_output><jupyter_text>To be consistent, the rounded $(u,v,w)$ values should add to zero just as the unrounded did. We can see that on the corners of the hex, this is not the case:<jupyter_code>summed = rounded[0] + rounded[1] + rounded[2]

fig, ax = plt.subplots()
pc = plot_hex_field(ax, summed)
fig.colorbar(pc)
diffs   = np.abs(rounded - UVW)
(du, dv, dw) = (diffs[i] for i in range(3))
u_worst = (du > dw) & (du > dv)
v_worst = ~u_worst & (dv > dw)
w_worst = ~u_worst & ~v_worst
(ru, rv, rw) = (rounded[i] for i in range(3))
ru[u_worst] = -rv[u_worst] - rw[u_worst]
rv[v_worst] = -ru[v_worst] - rw[v_worst]
rw[w_worst] = -ru[w_worst] - rv[w_worst]
summed = ru + rv + rw

fig, ax = plt.subplots()
pc = plot_hex_field(ax, summed)
fig.colorbar(pc)<jupyter_output><empty_output><jupyter_text>This shows that the rounded values satisfy the $u + v + w = 0$ relationship. Finally, we look at the rounded $(u,v,w)$ values themselves to see if they give us the correct coordinates:<jupyter_code>fig, axes = plt.subplots(ncols=3, figsize=(12, 2.5))
for (label, data, ax) in zip("uvw", (ru, rv, rw), axes):
    pc = plot_hex_field(ax, data)
    ax.set_title(label)
fig.colorbar(pc)<jupyter_output><empty_output><jupyter_text>Simpler 'rounding' function This is based on the "Branchless Method" described in http://www.redblobgames.com/grids/hexagons/more-pixel-to-hex.htmlneighbors> All the pixel to hex conversions I’ve seen use branches or a lookup table. I was mystified when <NAME> sent me pixel to hex conversion code that uses floor() five times, but no branches. First, divide x and y by size * sqrt(3); then find q, r with:>>     temp = floor(x + sqrt(3) * y + 1)    q = floor((floor(2*x+1) + temp) / 3);    r = floor((temp + floor(-x + sqrt(3) * y + 1))/3);Extracting the embedded expressions, it's clear that the floored expressions are coordinates along the $u$, $v$, and $w$ basis vectors:<jupyter_code>uvw_bases = normals[:3] / apothem
uvw_bases
A = np.floor(uvw_bases[0,0] * X + uvw_bases[0,1] * Y)
fig, ax = plt.subplots()
plot_hex_field(ax, A)
B = np.floor(uvw_bases[1,0] * X + uvw_bases[1,1] * Y)
fig, ax = plt.subplots()
plot_hex_field(ax, B)
C = np.floor(uvw_bases[2,0] * X + uvw_bases[2,1] * Y)
fig, ax = plt.subplots()
plot_hex_field(ax, C)<jupyter_output><empty_output><jupyter_text>Adding these is effectively like 'OR'ing the enclosed regions; dividing by 2 or 3 allows the non-integer regions to expand to fill parts of the hexes. Different combinations select different regions; we want regions that have isovalues inside the hexagon.<jupyter_code>fig, axes = plt.subplots(ncols=3, nrows=2, figsize=(10,6))
plot_hex_field(axes[0,0], (A + B) / 3)
plot_hex_field(axes[0,1], (B + C) / 3)
plot_hex_field(axes[0,2], (A + C) / 3)
plot_hex_field(axes[1,0], (A - B) / 3)
plot_hex_field(axes[1,1], (B - C) / 3)
plot_hex_field(axes[1,2], (C - A) / 3)<jupyter_output><empty_output><jupyter_text>We can translate these rhombi as needed; in this case we want to move the $(C - A)$ choice up to select the third corner of the hex in analog to $(A+B)$ and $(B + C)$...<jupyter_code>fig, axes = plt.subplots(ncols=3, nrows=1, figsize=(10,6))
plot_hex_field(axes[0], (A - B + 1) / 3)
plot_hex_field(axes[1], (B - C + 1) / 3)
plot_hex_field(axes[2], (A - C + 1) / 3)
fig, axes = plt.subplots(ncols=3, figsize=(10,3))
plot_hex_field(axes[0], (A + B))
plot_hex_field(axes[1], (A + B) / 3)
plot_hex_field(axes[2], np.ceil((A + B) / 3))
fig, axes = plt.subplots(ncols=4, figsize=(10,3))
plot_hex_field(axes[0], (B + C))
plot_hex_field(axes[1], (B + C) / 3)
plot_hex_field(axes[2], np.ceil((B + C) / 3))
fig, axes = plt.subplots(ncols=3, figsize=(10,3))
plot_hex_field(axes[0], (1 + A - C))
plot_hex_field(axes[1], (1 + A - C) / 3)
plot_hex_field(axes[2], np.floor((1 + A - C) / 3))<jupyter_output><empty_output><jupyter_text>This is much more straightforward than that cube-and-rounding system of doing things.<jupyter_code>uvw_bases = normals[:3]
def xy_to_uv(x, y):
    # Dot basis functions with xy coordinates
    (a, b, c) = np.array(uvw_bases.dot([x, y]) / apothem, dtype=int)
    u = -((-(a + b)) // 3)
    v = (1 + a - c) // 3
    return (u, v)
print(xy_to_uv(2.0, 0))
print(xy_to_uv(1.0, 2.0))
print(xy_to_uv(8.0, 10.0))
uvw_bases
uvw_bases.dot([0, 1])<jupyter_output><empty_output><jupyter_text>ShexagonalKeno's 'hexagonal' and 'triangular' are the same; they generate a rhombus-looking array with pointy tops. The 'shexagonal' option staggers the grid to look like a square layout.<jupyter_code>fig, ax = plt.subplots()

patches = []
umax, vmax = (3, 5)
(urangemax, vrangemax) = (umax, vmax)
if offset == 0.0: # pointy top
    urangemax += (vmax - 1) // 2
else:
    vrangemax += (umax - 1) // 2

for (u, v) in itertools.product(range(urangemax), range(vrangemax)):
    uv = np.array([u, v])
    xy = uv.dot(bases)
    fc = 'none'
    if offset == 0.0:
        uprime = u + v // 2 - (vmax - 1) // 2
        vprime = v
    else:
        vprime = v + u // 2 - (umax - 1) // 2
        uprime = u
    if (0 <= uprime < umax) and (0 <= vprime < vmax):
        fc = (0,0,0,.25)
    patches.append(RegularPolygon(xy, num_sides, outer_radius, poly_orient, fc=fc, edgecolor='k'))
    ax.text(xy[0], xy[1], "{},{}".format(uprime, vprime), ha='center', va='center')

ax.set_xlim(-1.2 * outer_radius, xy[0] + 1.2 * outer_radius)
ax.set_ylim(-1.2 * outer_radius, xy[1] + 1.2 * outer_radius)
ax.set_aspect('equal')
retick(ax)
ax.add_collection(PatchCollection(patches, match_original=True))
ax.grid()
fig.savefig("../doc/hex-{}-rect-{:d}x{:d}.pdf".format(hextype, umax, vmax))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17659878259562636, 'fraction_numerical': 0.02570075899902307, 'mean_word_length': 4.158139534883721, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 65, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '623345', 'license': 'Apache-2.0', 'n_tokens_mistral': 5545, 'n_tokens_neox': 5302, 'n_words': 1368, 'path': 'doc/.ipynbcheckpoints/hex-grid-checkpoint.ipynb', 'repo_name': 'celeritas-project/orange-port', 'size': '13325'}