Nexus: Optimal Tile Matrix Refactor and Documentation

After a chat with @jtkrogel and subsequently with @prckent, I took a look at the function `optimal_tilematrix()` in `structure.py`, and I think it could use some improvement and documentation.

As I do not work in condensed matter very often, I'm not entirely sure what this function does, but as someone who writes a lot of Python, I've noticed that it is not written in a Pythonic fashion and I think there is much room for improvement.

### The Code

To start, here is the entire section of code that is in `structure.py` that defines the variables and components required for `optimal_tilematrix()` to work (if I haven't missed anything):
```python
opt_tm_matrices = obj()
opt_tm_wig_indices = obj()

def trivial_filter(T):
return True
#end def trival_filter

class MaskFilter(DevBase):
def set(self,mask,dim=3):
omask = np.array(mask)
mask = np.array(mask,dtype=bool)
if mask.size==dim:
mvec = mask.ravel()
mask = np.empty((dim,dim),dtype=bool)
i=0
for mi in mvec:
j=0
for mj in mvec:
mask[i,j] = mi==mj
j+=1
#end for
i+=1
#end for
elif mask.shape!=(dim,dim):
error('shape of mask array must be {0},{0}\nshape received: {1},{2}\nmask array received: {3}'.format(dim,mask.shape[0],mask.shape[1],omask),'optimal_tilematrix')
#end if
self.mask = not mask
#end def set

def __call__(self,T):
return (T[self.mask]==0).all()
#end def __call__
#end class MaskFilter
mask_filter = MaskFilter()

def optimal_tilematrix(axes,volfac,dn=1,tol=1e-3,filter=trivial_filter,mask=None,nc=5,Tref=None):
if mask is not None:
mask_filter.set(mask)
filter = mask_filter
#end if
dim = 3
if isinstance(axes,Structure):
axes = axes.axes
else:
axes = np.array(axes,dtype=float)
#end if
if not isinstance(volfac,int):
volfac = int(np.around(volfac))
#end if
volume = np.abs(det(axes))*volfac
axinv = inv(axes)
cube = volume**(1./3)*np.identity(dim)
if Tref is None:
Tref = np.array(np.around(dot(cube,axinv)),dtype=int)
else:
Tref = np.asarray(Tref)
#end if
# calculate and store all tiling matrix variations
if dn not in opt_tm_matrices:
mats = []
rng = tuple(range(-dn,dn+1))
for n1 in rng:
for n2 in rng:
for n3 in rng:
for n4 in rng:
for n5 in rng:
for n6 in rng:
for n7 in rng:
for n8 in rng:
for n9 in rng:
mats.append((n1,n2,n3,n4,n5,n6,n7,n8,n9))
#end for
#end for
#end for
#end for
#end for
#end for
#end for
#end for
#end for
mats = np.array(mats,dtype=int)
mats.shape = (2*dn+1)**(dim*dim),dim,dim
opt_tm_matrices[dn] = mats
else:
mats = opt_tm_matrices[dn]
#end if
# calculate and store all wigner image indices
if nc not in opt_tm_wig_indices:
inds = []
rng = tuple(range(-nc,nc+1))
for k in rng:
for j in rng:
for i in rng:
if i!=0 or j!=0 or k!=0:
inds.append((i,j,k))
#end if
#end for
#end for
#end for
inds = np.array(inds,dtype=int)
opt_tm_wig_indices[nc] = inds
else:
inds = opt_tm_wig_indices[nc]
#end if
# track counts of tiling matrices
ntilings = len(mats)
nequiv_volume = 0
nfilter = 0
nequiv_inscribe = 0
nequiv_wigner = 0
nequiv_cubicity = 0
nequiv_shape = 0
# try a faster search for cells w/ target volume
det_inds_p = [
[(0,0),(1,1),(2,2)],
[(0,1),(1,2),(2,0)],
[(0,2),(1,0),(2,1)]
]
det_inds_m = [
[(0,0),(1,2),(2,1)],
[(0,1),(1,0),(2,2)],
[(0,2),(1,1),(2,0)]
]
volfacs = np.zeros((len(mats),),dtype=int)
for (i1,j1),(i2,j2),(i3,j3) in det_inds_p:
volfacs += (Tref[i1,j1]+mats[:,i1,j1])*(Tref[i2,j2]+mats[:,i2,j2])*(Tref[i3,j3]+mats[:,i3,j3])
#end for
for (i1,j1),(i2,j2),(i3,j3) in det_inds_m:
volfacs -= (Tref[i1,j1]+mats[:,i1,j1])*(Tref[i2,j2]+mats[:,i2,j2])*(Tref[i3,j3]+mats[:,i3,j3])
#end for
Tmats = mats[np.abs(volfacs)==volfac]
nequiv_volume = len(Tmats)
# find the set of cells with maximal inscribing radius
inscribe_tilings = []
rmax = -1e99
for mat in Tmats:
T = Tref + mat
if filter(T):
nfilter+=1
Taxes = dot(T,axes)
rc1 = norm(cross(Taxes[0],Taxes[1]))
rc2 = norm(cross(Taxes[1],Taxes[2]))
rc3 = norm(cross(Taxes[2],Taxes[0]))
r = 0.5*volume/max(rc1,rc2,rc3) # inscribing radius
if r>rmax or np.abs(r-rmax)rwmax or np.abs(rw-rwmax)0).sum()-(o<0).sum()
if s>smax or np.abs(s-smax)0:
cells = diagonal
elif len(symmetric)>0:
cells = symmetric
elif len(antisymmetric)>0:
cells = antisymmetric
s = -1
elif len(other)>0:
cells = other
else:
cells = []
#end if
skew_min = 1e99
if len(cells)>0:
for rw,T,Taxes in cells:
Td = np.diag(np.diag(T))
skew = np.abs(T.T-s*T-(1-s)*Td).sum()
if skew