[Revised match factor calculation]

    .. todo:: Complete geom_symm_match docstring

def geom_symm_match(g, atwts, ax, theta, do_refl):
    """ [Revised match factor calculation]

    .. todo:: Complete geom_symm_match docstring

    """

    # Imports
    import numpy as np
    from scipy import linalg as spla

    # Convert g and atwts to n-D vectors
    g = make_nd_vec(g, nd=None, t=np.float64, norm=False)
    atwts = make_nd_vec(atwts, nd=None, t=np.float64, norm=False)

    # Ensure proper dimensionality
    if not g.shape[0] == 3 * atwts.shape[0]:
        raise ValueError("Size of 'g' is not 3*size of 'atwts'")
    ## end if

    # Calculate transformed geometry
    gx = symm_op(g, ax, theta, do_refl)

    # Push g to a column vector
    g = g.reshape((g.shape[0],1))

    # Augment g and gx with imaginary atomic weights
    ex_wts = atwts.repeat(3,axis=0).T.reshape((atwts.shape[0]*3,1)) * 1.j
    g = np.add(g, ex_wts)
    gx = np.add(gx, ex_wts)

##        # Define calc as the outer product of the augmented vectors
##        calc = np.dot(g.reshape((g.shape[0],1)), \
##                            np.reciprocal(gx.reshape((1,gx.shape[0]))))
##
##        # Calculate the complex magnitude of each element and take log10,
##        #  then abs again
##        calc = np.abs(np.log10(np.abs(calc)))

    # Expand g and gx as column vectors of coordinates
    calc_g = g.reshape((g.shape[0] // 3, 3))
    calc_gx = gx.reshape((gx.shape[0] // 3, 3))
##
##        # Expand each into a square matrix of identical column vectors
##        calc_g = calc_g.repeat(g.shape[0], axis=1)
##        calc_gx = gx.repeat(gx.shape[0], axis=1)

    # Calc is the absolute distance between the calc-ed values,
    #  scaled by the maximum of the individual atom distances or unity.

    # Calculate the unscaled distances
    calc = [[spla.norm(np.subtract(calc_g[i,:], calc_gx[j,:])) \
                            for j in range(calc_gx.shape[0])] \
                            for i in range(calc_g.shape[0])]

    # Calculate the scale factors
    scale_g = np.array([spla.norm(calc_g[i,:]) for i in \
                range(calc_g.shape[0])]).reshape((calc_g.shape[0],1)) \
                .repeat(calc_g.shape[0], axis=1)
    scale_gx = np.array([spla.norm(calc_gx[j,:]) for j in \
                range(calc_g.shape[0])]).reshape((1,calc_gx.shape[0])) \
                .repeat(calc_gx.shape[0], axis=0)
    scale = np.maximum(np.maximum(scale_g, scale_gx),
                np.ones_like(scale_g, dtype=np.float64))

    # Scale calc
    calc = np.divide(calc, scale)

    # Take the minimum of each row
    mins = np.min(calc, axis=1)

    # Take the maximum of the minima for the final factor
    fac = np.max(mins)

    # Using the atomic weights for checking matching can result in 'fac'
    #  being greater than unity. Return the minimum of fac and unity.
    fac = min(fac, 1.0)
    return fac