Sample given points stochastically in a 3D space with minimum nearest-neighbor distance and maximum density

I have n points in a 3D space. I want to stochastically sample a subset of points with all nearest-neighbor distances larger than r. The size of the subset m is unknown, but I want the sampled points to be as dense as possible.

There are similar questions, but they are all about generating points, rather than sampling from given points.
Generate random points in 3D space with minimum nearest-neighbor distance

Generate 3-d random points with minimum distance between each of them?

Say I have 300 random 3D points,

import numpy as np
n = 300
points = np.random.uniform(0, 10, size=(n, 3))

What is the fastest way to get a subset of m points with minimum nearest-neighbor distance r = 1 while maximizing m?

Testing @David Eisenstat 's answer:

from ortools.linear_solver import pywraplp

solver = pywraplp.Solver.CreateSolver("SCIP")
variables = [solver.BoolVar("x[{}]".format(i)) for i in range(n)]
solver.Maximize(sum(variables))
for j, q in enumerate(points):
    for i, p in enumerate(points[:j]):
        if np.linalg.norm(p - q) <= 1:
            solver.Add(variables[i] + variables[j] <= 1)
solver.EnableOutput()
solver.Solve()
indices = [i for (i, variable) in enumerate(variables) if variable.SolutionValue()]
subset = points[indices]

Plot the subset:

x, y, z = subset[:,0], subset[:,1], subset[:,2]
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, zdir='z', c= 'red')
plt.show()

Output:

Check minimum pairwise distance:

from scipy.spatial.distance import cdist
dist = cdist(subset, subset, metric='euclidean')
np.fill_diagonal(dist, 100) # Fill diagonal with a large number
print(np.min(dist))

Output:

1.0044196918966806

from Sample given points stochastically in a 3D space with minimum nearest-neighbor distance and maximum density