%pylab inline
from uncertainties import ufloat

import tess

nm = 8
ptsr = [tuple(array((l,m,n))+.5) + ((l+m+n) % 2,)
       for l in range(nm) for m in range(nm) for n in range(nm)]
ptsr = np.array(ptsr)
pts = ptsr[:, :3]
r = [.5 if r < 1 else .6 for r in ptsr[:, 3]]

cells1 = tess.Container(pts, nm, periodic=False)
len(cells1)

cells2 = tess.Container(pts, nm, radii=r, periodic=True)
len(cells2)

# a useful function for later
def ufloatms(arr):
    "Return a ufloat based on the mean and standard deviation of an array"
    from uncertainties import ufloat
    import numpy as np
    mn, sd = np.mean(arr), np.std(arr)
    return ufloat(mn, sd)

for cells in (cells1, cells2):
    alln = []
    for v in cells:
        alln.append(len(v.neighbors()))

    print(ufloatms(alln))
    for n in sorted(set(alln)):
        print('{:3d}:{:3d}'.format(n, alln.count(n)))

for cells in (cells1, cells2):
    allVs = []
    for v in cells:
        allVs.append(v.volume())

    print(ufloatms(allVs))
    vlst = ['{:8.3g}'.format(V) for V in allVs]
    for n in sorted(set(vlst)):
        print('{:8s}:{:3d}'.format(n, vlst.count(n)))

from mpl_toolkits.mplot3d import Axes3D

for cells in (cells1, cells2):
    (fig, (ax1, ax2)) = plt.subplots(1,2,figsize=(10,4), subplot_kw=dict(projection='3d'))

    p1, p2 = array((1.5,1.5,1.5)), array((1.5,1.5,2.5))
    for v in cells:
        if np.allclose(p1, v.pos()):
            ax = ax1
        elif np.allclose(p2, v.pos()):
            ax = ax2
        else:
            continue

        x1,y1,z1 = v.pos()
        xyz2 = [cells[m].pos() for m in v.neighbors()]

        for x2,y2,z2 in sorted(xyz2):
            ls = 'k-'
            if z1 - z2 > z1*1e-3: ls = 'r-'
            elif z1 - z2 < -z1*1e-3: ls = 'b-'
            ax.plot([x1,x2], [y1,y2], [z1,z2], ls)

c, c2 = cells[:2]

c.pos(), c2.pos()

np.around(c.centroid(), 5)

c.volume()

c.surface_area()

c.number_of_faces()

c.max_radius_squared()

c.total_edge_distance()

c.face_areas()

np.around(c.normals(), 8)

expected_coords = [-1.0, -1/sqrt(2.), 0, 1/sqrt(2.0), 1.0]
found_coords = sorted(set(np.around(np.array(c2.normals()).flatten(), 12)))
np.allclose(found_coords, expected_coords)

cells.order(l=6, local=False, weighted=True)