<jupyter_start><jupyter_text>Quasi-uniform grids on 3D sphere<NAME>  11.2021Usage:```python 'ipynb as module' loaderimport ipynb_loader import ipynb itselfimport quasi_uniform_3d as qu3dN = 10 Generating pointscube_points = qu3d.generate_cube_uniform_grid(N) sphere_points = qu3d.generate_sphere_uniform_grid(N) sphere_ponts_2 = qu3d.generate_spherical_uniform_grid(N) Distance statistics: maximal distance between closest points, minimal distance,   and histogram of distances between closest points (point-wise)density, tightness, min_dists = get_grid_density(sphere_points)```See also examples in [quasi-uniform-3d-demo.ipynb](quasi-uniform-3d-demo.ipynb)1. Quusi-uniform gridD. <NAME>, Uniform spherical grids via equal area projection from the cube to the sphereJournal of Computationanl and Applied Mathimatics, 236:6 (2011), 1033-1041.2. Spherical coordinatesSee also: https://medium.com/@oscarsc/four-ways-to-create-a-mesh-for-a-sphere-d7956b825db4<jupyter_code>import numpy as np
# Rosca-Plonka 2011 
# Uniform spherical grids via equal area projection from the cube to the sphere

def get_XY_coordinates(points):
    """Map points (x, y) of a square [-a, a]^2 to a curved square. <NAME>"""
    # Size "a" of the square shall be chosen prior to calling this function
    # TODO: process (0, 0) point as it causes warning
    x, y = points.T
    
    ind_x_gte_y = (np.abs(x) >= np.abs(y))
    ind_y_gt_x = np.logical_not(ind_x_gte_y)
    
    X = np.zeros_like(x)
    Y = np.zeros_like(x)
    
    # constants, common parts and shortcuts
    c2_14_over_beta = 2**0.25 / np.sqrt(np.pi / 6)
    sqrt_2 = np.sqrt(2)
    
    y_on_x_1 = y[ind_x_gte_y] / x[ind_x_gte_y]
    cos_a = np.cos(y_on_x_1 * np.pi / 12)
    sin_a = np.sin(y_on_x_1 * np.pi / 12)
    
    x_on_y_2 = x[ind_y_gt_x] / y[ind_y_gt_x]
    cos_b = np.cos(x_on_y_2 * np.pi / 12)
    sin_b = np.sin(x_on_y_2 * np.pi / 12)
    
    X[ind_x_gte_y] = c2_14_over_beta * x[ind_x_gte_y] * (sqrt_2 * cos_a - 1) / np.sqrt(sqrt_2 - cos_a)
    Y[ind_x_gte_y] = c2_14_over_beta * x[ind_x_gte_y] * (sqrt_2 * sin_a) / np.sqrt(sqrt_2 - cos_a)
    
    X[ind_y_gt_x] = c2_14_over_beta * y[ind_y_gt_x] * (sqrt_2 * sin_b) / np.sqrt(sqrt_2 - cos_b)
    Y[ind_y_gt_x] = c2_14_over_beta * y[ind_y_gt_x] * (sqrt_2 * cos_b - 1) / np.sqrt(sqrt_2 - cos_b)

    return X, Y
    
def get_xyz_from_XY_curved(X, Y):
    """Map [X, Y] points of a curved square [X, Y, 1] to a sphere upper hat [x, y, z]"""
    c = 1 - (X**2 + Y**2) / 4
    x = np.sqrt(c) * X
    y = np.sqrt(c) * Y
    z = 2 * c - 1
    
    return x, y, z # [zip(x, y, z)]

def get_xyz_from_XY_flat(X, Y):
    """Map [X, Y] points of a square [X, Y, 1] to a cube upper hat [x, y, z]"""
    x = X
    y = Y
    z = np.ones_like(X)
    
    return x, y, z # [zip(x, y, z)]


# # Case 1: grid with overlapping

# # points = np.random.rand(2, 100)
# def generate_uniform_2d_overlap_grid(N):
#     """Generate N^2 points (take even N) on a [-1, 1]^2 square of XY-plane"""
#     p1, p2 = np.meshgrid(np.linspace(-1, 1, N), np.linspace(-1, 1, N))
# #     print(f'size of p1 = {p1.shape}, p2 = {p2.shape}')
#     points = np.vstack([np.ravel(p1), np.ravel(p2)]).T
# #     print(f'size of points = {points.shape}')
    
#     return points

# def generate_sphere_uniform_overlap_grid_by_xyz(x, y, z):
#     """From a uniform grid on sphere hat make fully-covered uniform sphere (with some overlapping)"""
#     uniform_dir_points = np.vstack([
#         np.stack([x, y, z]).T,
#         np.stack([-x, -y, -z]).T, # symmetric cover
#         np.stack([z, x, y]).T,
#         np.stack([-z, x, y]).T,
#         np.stack([x, z, y]).T,
#         np.stack([x, -z, y]).T
#     ])
#     return uniform_dir_points

# def generate_sphere_uniform_overlap_grid(N):
#     """Generate quasi uniform grid on unit sphere with 6 N^2 points (with overlapping)"""
#     b = np.sqrt(np.pi / 6)
#     scaled_flat_points = b * generate_uniform_2d_overlap_grid(N)
#     X, Y = get_XY_coordinates(scaled_flat_points)
#     x, y, z = get_xyz_from_XY_curved(X, Y) # sphere "hat"
#     uniform_dir_points = generate_sphere_uniform_overlap_grid_by_xyz(x, y, z)
    
#     return uniform_dir_points
# Case 2: non-overlapping grid

def generate_uniform_2d_grid_type_1(N):
    """Generate N (N - 2) points (take even N >= 4) without top and bottom stripes on a [-1, 1]^2 square of XY-plane"""
    # first index - "x", second - "y"
    p1, p2 = np.meshgrid(np.linspace(-1, 1, N), np.linspace(-1, 1, N), indexing='ij')
    # exclude lowest and highest point lines (in y)
    p1_cut = p1[:, 1:-1]
    p2_cut = p2[:, 1:-1]
    
#     print(f'size of p1_cut = {p1_cut.shape}, p2_cut = {p2_cut.shape}')
    points = np.vstack([np.ravel(p1_cut), np.ravel(p2_cut)]).T
#     print(f'size of points = {points.shape}')
    return points

def generate_uniform_2d_grid_type_2(N):
    """Generate N (N - 2) + 4 points (take even N >= 4) without top and bottom stripes but with corneres of XY-plane"""
    # first index - "x", second - "y"
    stripe_points = generate_uniform_2d_grid_type_1(N)
    
    # add 4 corner points
    points = np.vstack([ 
        stripe_points,
        np.array([[-1, -1], [-1, 1], [1, -1], [1, 1]])
        ])
#     print(f'size of resulting array = {points.shape}')
    return points

def sphere_patching(type_1_flat_points, type_2_flat_points):
    """Glue a sphere from non-overlapping flat patches"""
    b = np.sqrt(np.pi / 6)
    
    scaled_flat_points_type_1 = b * type_1_flat_points
    X_type_1, Y_type_1 = get_XY_coordinates(scaled_flat_points_type_1)
    x_type_1, y_type_1, z_type_1 = get_xyz_from_XY_curved(X_type_1, Y_type_1)

    scaled_flat_points_type_2 = b * type_2_flat_points
    X_type_2, Y_type_2 = get_XY_coordinates(scaled_flat_points_type_2)
    x_type_2, y_type_2, z_type_2 = get_xyz_from_XY_curved(X_type_2, Y_type_2)
    
    uniform_dir_points = np.vstack([
        np.stack([x_type_2, y_type_2, z_type_2]).T, # top
        np.stack([y_type_1, z_type_1, x_type_1]).T, # front
        np.stack([z_type_1, x_type_1, y_type_1]).T, # left
        np.stack([-z_type_1, -x_type_1, -y_type_1]).T, # right
        np.stack([-y_type_1, -z_type_1, -x_type_1]).T, # back
        np.stack([-x_type_2, -y_type_2, -z_type_2]).T # bottom
    ])
    return uniform_dir_points

def cube_patching(type_1_points, type_2_points):
    """Glue a sphere from non-overlapping patches. <NAME>"""
    x_type_1, y_type_1, z_type_1 = get_xyz_from_XY_flat(*type_1_points.T)
    x_type_2, y_type_2, z_type_2 = get_xyz_from_XY_flat(*type_2_points.T)
    
    uniform_dir_points = np.vstack([
        np.stack([x_type_2, y_type_2, z_type_2]).T, # top
        np.stack([y_type_1, z_type_1, x_type_1]).T, # front
        np.stack([z_type_1, x_type_1, y_type_1]).T, # left
        np.stack([-z_type_1, -x_type_1, -y_type_1]).T, # right
        np.stack([-y_type_1, -z_type_1, -x_type_1]).T, # back
        np.stack([-x_type_2, -y_type_2, -z_type_2]).T # bottom
    ])
    return uniform_dir_points

def generate_cube_uniform_grid(N):
    """Generate uniform grid on unit cube with 6 N (N - 2) 6 + 8 points (without overlapping). <NAME>"""
    flat_points_type_1 = generate_uniform_2d_grid_type_1(N)
    flat_points_type_2 = generate_uniform_2d_grid_type_2(N)
    uniform_cube_points = cube_patching(flat_points_type_1, flat_points_type_2)
    
    return uniform_cube_points
    
def generate_sphere_uniform_grid(N):
    """Generate quasi uniform grid on unit sphere with 6 N (N - 2) + 8 points (without overlapping). <NAME>"""
    flat_points_type_1 = generate_uniform_2d_grid_type_1(N)
    flat_points_type_2 = generate_uniform_2d_grid_type_2(N)
    
    uniform_sphere_points = sphere_patching(flat_points_type_1, flat_points_type_2)
    
    return uniform_sphere_points<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13309352517985612, 'fraction_numerical': 0.0341726618705036, 'mean_word_length': 3.4007914075749013, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 14, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, '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': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '470267', 'license': 'CC0-1.0', 'n_tokens_mistral': 3314, 'n_tokens_neox': 3174, 'n_words': 700, 'path': 'quasi_uniform_3d.ipynb', 'repo_name': 'atremba/reach_set', 'size': '7829'}