I am interested in the following problem, and am curious if Mathematica could be used to solve it.

I start with a radius function $\rho:D\to \mathbf{R}_{\ge 0}$ on a region $D\subseteq \mathbf{R}^n$. Fix $r_2>r_1$ and a finite subset $S_0\subset D$ of points.

Problem: I am interested randomly in generating a fairly large set $S\supset S_0$ of dots $x_i$ in $D$ such that $$|x_i-x_j|\ \not\in\ [\rho(x_i)r_1,\rho(x_i)r_2]$$ for all $i\ne j$.

There are lots of example uses of Mathematica generating a random set of points, but I am not sure how to alter these solutions to include this condition on the distances between the points.

To be extremely precise for the sake of argument, take for instance $D=[0,1]^2\subseteq \mathbf{R}^2$ with $r_1=1$, $r_2=2$ and radius function $\rho(x,y)=x^2+xy+y^2.$

  • 1
    $\begingroup$ A specific example with $n=2$ or $n=3$ is essential to include. What you've tried (even if you didn't get too far) with that example in Mathematica would also be helpful. $\endgroup$
    – JimB
    Commented Apr 8, 2022 at 21:12
  • $\begingroup$ Shouldn't your density function integrate to 1 over $[0,1]^2$ ? Integrate[x^2 + x y + y^2, {x, 0, 1}, {y, 0, 1}] = 11/12. $\endgroup$
    – JimB
    Commented Apr 8, 2022 at 21:30
  • $\begingroup$ Apologies- it was a bad name, $\rho$ only appears in the radius condition rather than the distribution using which the points are sampled. I have edited my question to call $\rho$ a radius function instead. $\endgroup$
    – Pulcinella
    Commented Apr 8, 2022 at 21:32
  • $\begingroup$ So I'm a little slow: if $\rho$ is not a density function, what it the distribution from which the points are sampled? $\endgroup$
    – JimB
    Commented Apr 8, 2022 at 21:34
  • 1
    $\begingroup$ While it's not computationally efficient, you could do "sample rejection". Generate points from the whole interval, then check them to see if they violate your constraint. If they do, throw them away, else keep them. $\endgroup$
    – bill s
    Commented Apr 9, 2022 at 15:29

2 Answers 2


This (I hope after discussion with OP), is closer to what is desired. This does not scale well (not >100 points). I just use euclidean norm but this could be changed.

I only post it in case it motivates others.

rand[n_, a_, b_] := Module[{r = RandomReal[{0, 1}, {n, 2}], rg, cl},
  rg = RelationGraph[
    Not[a < EuclideanDistance[#1, #2] < b && 
       EuclideanDistance[#1, #2] != 0] &, r];
  cl = First@FindClique[rg];
  Graphics[{Red, Point[r], Black, PointSize[0.03], Point[#]} & /@ cl, 
   Axes -> True, PlotRange -> {{0, 1}, {0, 1}}, Frame -> True]]

The red points are the original points and the black over red points are those that comply with constraint.

Illustrating this:

 With[{plot = rand[100, a, a + b]}, 
    Graphics[{Circle[p, a], Circle[p, a + b]}]], {p, {0.5, 0.5}, 
    Locator, Appearance -> "+"}]], {a, {0.05, 0.1, 0.2}, 
  PopupMenu}, {b, {0.1, 0.2, 0.3}, PopupMenu}]

enter image description here


Here is a brute-force approach (which is essentially a poor copy of @ubpdqn 's answer):


(* Set parameters *)
f[r_, {x_, y_}] := r (x^2 + x y + y^2)
r1 = 0.1;
r2 = 0.3;

(* Number of points to select *)
n = 50;

(* Array to store sequential plots *)
plt = ConstantArray[0, n];

(* Initial point *)
pts = {RandomVariate[UniformDistribution[{0, 1}], 2]};
noPts = {};
m = 100;  (* Fineness of grid showing area where no points can be selected *)
(* Create a list of points on a fine grid that can no longer be selected *)
Do[Do[status = True;
  If[f[r1, {i/m, j/m}] <= Norm[{i/m, j/m} - pts[[1]]] <= f[r2, {i/m, j/m}], status = False];
  If[status == False, noPts = Join[noPts, {{i/m, j/m}}]], {i, m}], {j, m}];
plt[[1]] = 
  ListPlot[{pts, noPts, pts}, PlotRange -> {{0, 1}, {0, 1}}, AspectRatio -> 1, 
   PlotStyle -> {{Blue, PointSize[0.01]}, {Gray, PointSize[0.005]}, {Red, PointSize[0.02]}}, Frame -> True];

(* Selete points 2 through n *)
Do[good = False; (* Initial status of potential sample point *)
 While[! good,
  {x0, y0} = RandomVariate[UniformDistribution[{0, 1}], 2];
  status = True;  (* 
  Assume status of potential sample point is true until proven otherwise *)
  (* Check on relationship with all other points *)
  Do[If[f[r1, {x0, y0}] <= Norm[{x0, y0} - pts[[j]]] <= 
     f[r2, {x0, y0}], status = False], {j, 1, h - 1}];
  If[status, good = True; pts = Join[pts, {{x0, y0}}]]];
 (* Highlight locations where the next point cannot occur *)
 noPts = {};
 m = 100;
 Do[Do[status = True;
   Do[If[f[r1, {i/m, j/m}] <= Norm[{i/m, j/m} - pts[[k]]] <= f[r2, {i/m, j/m}], status = False], 
   {k, h}];
   If[status == False, noPts = Join[noPts, {{i/m, j/m}}]], {i, m}], {j, m}];
 (* Create plot of currently sampled points and the area that can no longer be sampled *)
 plt[[h]] = ListPlot[{pts, noPts, {pts[[h]]}}, PlotRange -> {{0, 1}, {0, 1}}, 
   AspectRatio -> 1, 
   PlotStyle -> {{Blue, PointSize[0.01]}, {Gray, PointSize[0.005]}, {Red, PointSize[0.02]}}, Frame -> True],
 {h, 2, n}]

Sequential sampling of points

The red dot represents the currently sampled point and the blue dots are all of the previously sampled points. The smaller gray dots represent the area than can no longer be sampled due to the restrictions imposed.

  • $\begingroup$ like it…have some other thoughts but too busy right now to test…mine not efficient but got ‘ball rolling’ :) $\endgroup$
    – ubpdqn
    Commented Apr 11, 2022 at 6:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.