Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

The Mathematica 10 documentation was updated for FindInstance adding support for regions.

In my use case I'm trying to sample points in a set of disks:

region = DiscretizeRegion@RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}]
FindInstance[{x, y} ∈ region, {x, y}, Reals, 2] // N

However the above code fails and generates the following error:

FindInstance::elemc: "Unable to resolve the domain or region membership condition {x,y}∈"

What's going wrong here?

share|improve this question
    
Thanks @RunnyKine how can we escalate it to get a response? –  user5601 Aug 22 at 19:39
    
It would be nice if more people reported it to them and linked to these questions here like I did. –  RunnyKine Aug 22 at 19:40
    
FindInstance is not for generating random numbers ... If you already have the disks, it's trivial to find an instance of a point which lies within them. Take the centre of any disk. It doesn't even matter that you have several disks, you can just use a single one. If you have any expectation that these numbers will be uncorrelated and uniformly distributed (within the disks) then do not use FindInstance. –  Szabolcs Aug 22 at 20:09
    
@Szabolcs, while I agree with your comment above, that's not the problem here. The problem is that these numerical functions that have been advertised to work with regions don't work with MeshRegions or BoundaryMeshRegions e.g. see here. –  RunnyKine Aug 22 at 20:32
    
@RunnyKine it has nothing to do with MeshRegion, try it and see (remove DiscretizeRegion). Its weird, and I'm reporting it. –  rcollyer Aug 22 at 20:53

6 Answers 6

It would be nice if UniformDistribution worked on arbitrary regions, then we could simply do RandomVariate[UniformDistribution[region]]. Someone at Wolfram should get on that.

In the meantime, it seems we have to write our own sampling routines. @m_goldberg's answer is very nice (vote it up!) and uses rejection sampling, which works for arbitrary regions. However it will become slow if the measure of the region is small compared to that of its bounding box, for example if the disks are small and far apart. On the other hand, since you have a MeshRegion, it is actually possible to generate uniformly distributed samples without any rejection. We will pick a mesh cell randomly with probability proportional to its measure, then pick a point uniformly within that cell.

Sorry, I haven't bothered to wrap the code up into a neat module.

SeedRandom[0];
region = DiscretizeRegion@RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}]

enter image description here

d = RegionDimension[region];
coords = MeshCoordinates[region];
cells = MeshCells[region, d];
cellMeasures = PropertyValue[{region, d}, MeshCellMeasure];

randomCell[] := RandomChoice[cellMeasures -> cells]
randomBarycentric[] := Differences@Flatten@{0, Sort@RandomReal[1, d], 1}
randomRegionPoint[] := randomBarycentric[].coords[[First@randomCell[]]]

ListPlot[Table[randomRegionPoint[], {10000}], AspectRatio -> Automatic]

enter image description here

share|improve this answer
    
This is very nice upgrade version of goldberg's answer –  Junho Lee Aug 23 at 23:02
3  
You can collect the already available cellMeasures like this: PropertyValue[{region, 2}, MeshCellMeasure] –  RunnyKine Aug 23 at 23:26
    
@RunnyKine: Thanks, added! –  Rahul Aug 23 at 23:59
1  
@RahulNarain very nice...learned a lot...+1 –  ubpdqn Aug 24 at 8:34

There are already good answers, but I'm going to improve the performance, generalize to any region in any dimensions and make the function more convenient. The main idea is to use DirichletDistribution (the uniform distribution on a simplex, e.g. triangle or tetrahedron). This idea was implemented by PlatoManiac and me in the related question obtaining random element of a set given by multiple inequalities (there is also Metropolis algorithm, but it is not suitable here).

The code is relatively short:

RegionDistribution /: 
 Random`DistributionVector[RegionDistribution[reg_MeshRegion], n_Integer, prec_?Positive] :=
  Module[{d = RegionDimension@reg, cells, measures, s, m},
   cells = Developer`ToPackedArray@MeshPrimitives[reg, d][[All, 1]];
   s = RandomVariate[DirichletDistribution@ConstantArray[1, d + 1], n];
   measures = PropertyValue[{reg, d}, MeshCellMeasure];
   m = RandomVariate[#, n] &@EmpiricalDistribution[measures -> Range@Length@cells];
   #[[All, 1]] (1 - Total[s, {2}]) + Total[#[[All, 2 ;;]] s, {2}] &@
   cells[[m]]]

Examples

Random disks (2D in 2D)

SeedRandom[0];
region = DiscretizeRegion@RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}];
pts = RandomVariate[RegionDistribution[region], 10000]; // AbsoluteTiming
ListPlot[pts, AspectRatio -> Automatic]

{0.004473, Null}

enter image description here

Precise test

pts = RandomVariate[RegionDistribution[region], 200000000]; // AbsoluteTiming

{85.835022, Null}

Histogram3D[pts, 50, "PDF", BoxRatios -> {Automatic, Automatic, 1.5}]

enter image description here

It is fast for $2\cdot10^8$ points and the distribution is really flat!

Intervals (1D in 1D)

region = DiscretizeRegion[Interval[{0, 1}, {2, 4}]];
pts = RandomVariate[RegionDistribution[region], 100000]; // AbsoluteTiming
Histogram[Flatten@pts]

{0.062430, Null}

enter image description here

Random circles (1D in 2D)

region = DiscretizeRegion@RegionUnion[Circle /@ RandomReal[10, {100, 2}]];
pts = RandomVariate[RegionDistribution[region], 10000]; // AbsoluteTiming
ListPlot[pts, AspectRatio -> Automatic]

{0.006216, Null}

enter image description here

Balls (3D in 3D)

region = DiscretizeRegion@RegionUnion[Ball[{0, 0, 0}], Ball[{1.5, 0, 0}], Ball[{3, 0, 0}]];
pts = RandomVariate[RegionDistribution[region], 10000]; // AbsoluteTiming
ListPointPlot3D[pts, BoxRatios -> Automatic]

{0.082202, Null}

enter image description here

Surface cow disctribution (2D in 3D)

region = DiscretizeGraphics@ExampleData[{"Geometry3D", "Cow"}];
pts = RandomVariate[RegionDistribution[region], 2000]; // AbsoluteTiming
ListPointPlot3D[pts, BoxRatios -> Automatic]

{0.026357, Null}

enter image description here

Line in space (1D in 3D)

region = DiscretizeGraphics@ParametricPlot3D[{Sin[2 t], Cos[3 t], Cos[5 t]}, {t, 0, 2 π}];
pts = RandomVariate[RegionDistribution[region], 1000]; // AbsoluteTiming
ListPointPlot3D[pts, BoxRatios -> Automatic]

{0.005056, Null}

enter image description here

share|improve this answer
1  
Very nice (+1). –  kale Oct 9 at 13:53
    
@ybeltukov Very nice answer. Could you just explain briefly why you use DirichletDistribution? –  s.s.o Oct 9 at 14:33
1  
@s.s.o It is the uniform distribution on a simplex (e.g. triangle, tetrahedron). –  ybeltukov Oct 9 at 15:42
    
Ah, I should have known my randomBarycentric function was (the special case of) a well-known probability distribution. –  Rahul Oct 9 at 19:14
SeedRandom[0];

region = RegionUnion @ Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}];

DiscretizeRegion @ region

enter image description here

points = RandomReal[{-1, 5}, {10000, 2}];

circles = List @@ (region /. Disk -> Circle);

Graphics[{AbsolutePointSize[2],
  Transpose[{RegionMember[region, points] /.
     {False -> White, True -> Gray}, Point /@ points}], circles}, Frame -> True]

enter image description here

share|improve this answer
    
thanks, but (region = RegionUnion@ Table[Disk[RandomReal[{0, 100}, {2}], RandomReal[1]], {10000}]; RegionMember[region, RandomReal[{0, 100}, {100000, 2}]]) crashes the kernel –  user5601 Aug 22 at 20:46
2  
@user5601 Your question is "How to generate random points in a region". It's not about FindInstance (see several comments), and it's not about the endless possibilities of "crashing the kernel". I answered your question using your own definition of region. –  eldo Aug 22 at 21:05
    
Right! I should have specified that I need many points! –  user5601 Aug 25 at 15:08
    
region = RegionUnion @ Table[Disk[RandomReal[4, {2}], DiscretizeRegion is very slow, so if the number of disks exceeds 10000 then it freezes... –  user5601 Aug 25 at 15:09

Here is a work-around that will let you extract as many randoms points from a region as you wish.

SeedRandom[0]; 
region = 
  DiscretizeRegion @ RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}]

region

randomFromRegion[region_, trials_] :=
  Module[{bounds, randPts},
    bounds = RegionBounds@region;
    randPts = 
      Transpose @ {RandomReal[bounds[[1]], trials], RandomReal[bounds[[2]], trials]};
    Select[randPts, RegionMember[region]]]

Let's see what we get from 5000 trials.

pts = randomFromRegion[region, 5000];
Length @ pts
2550

We extracted 2550 points that lie within the region. Here is how they are distributed.

Graphics @ Point @ pts

points

share|improve this answer
    
RegionBounds make rectangle region from circle. From the beginning, You ought to have made rectangle region. –  Junho Lee Aug 22 at 21:20
1  
Sorry, but I don't understand your objection. –  m_goldberg Aug 22 at 21:30
    
@JunhoLee. If you were referring to my assertion concerning approximating area, I agree it was erroneous and have removed it. –  m_goldberg Aug 22 at 22:42

FindInstance have to apply to a expression like an equation or an inequality. DiscretizeRegion does not make an inequality. You should apply like the following example.

region = RegionUnion@Table[Disk[RandomReal[4, 2], RandomReal[1]],
    {10}
    ];
DiscretizeRegion[region]

Blockquote

RegionQ[region]

True

FindInstance[RegionMember[region, {x, y}], {x, y}, Reals, 2] // N

{{x -> 2.4597, y -> 4.09129}, {x -> 1.31552, y -> 2.00384}}

share|improve this answer
    
But DiscretizeRegion is a region. Because RegionQ[] returns True –  RunnyKine Aug 22 at 20:48
    
Seems like it is too slow: (region = RegionUnion@ Table[Disk[RandomReal[{0, 100}, {2}], RandomReal[1]], {10000}]; FindInstance[RegionMember[region, {x, y}], {x, y}, Reals, 2] // N) –  user5601 Aug 22 at 20:48
    
@user5601 Please check technique applying method of FindInstance. I might guess you find the solution like this posting example. –  Junho Lee Aug 22 at 20:50
    
@RunnyKine you right. DiscretizeRegion is not a quantified system of equations and inequalities. –  Junho Lee Aug 22 at 20:55
    
@user5601 RegionMember[region, {x, y}] make inequalities. too slow is natural.. –  Junho Lee Aug 22 at 20:56

you can look at this question here

This method also may work fine:

p1 = {x, y} /. FindInstance[x^2 + y^2 < 1, {x, y}, Reals, 1000];
p2 = {y, x} /. FindInstance[x^2 + y^2 < 1, {x, y}, Reals, 1000];
point = RandomChoice[Join[p1, p2], 1000];
Graphics[{Circle[], {Red, Point[point]}}]
share|improve this answer
1  
FindInstance is not a random number generator. There's no guarantee that it will give you uncorrelated random numbers with a uniform distribution. What the OP is asking for a misuse of FindInstance that is going to lead to wrong results. –  Szabolcs Aug 22 at 20:12
    
I know FindInstance is not a random number generator but I think they will be uniformly distributed even if the pointes are correlated. then I used RandomChoice to get random pointes. I don't know if this way is correct or not. –  Algohi Aug 23 at 17:02

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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