9
$\begingroup$

Another great video from 3Blue1Brown describes a neat proof (at about 5:45) that this has to be 25%. I wanted to do the experiment in Mathematica as a way of showing how exploring the problem computationally might give you a clue about the proof.

While I'm happy I got there in the end, the timing is really bad, about 10 seconds for 3000 triangles. Any hints?

RegionMember[#, 
   origin] & /@ (Table[{Cos[RandomReal[2 Pi]], 
       Sin[RandomReal[2 Pi]]}, 9000] // Partition[#, 3] & // 
    Map[Polygon, #] &) // Count[#, True] &
$\endgroup$
  • $\begingroup$ Your solution suggests that the probability that the origin is inside your triangles is 1/4. That is interesting, because your triangles are NOT random triangles on the unit circle. Your points have the form {Cos[a], Sin[b]} where a and b are different random numbers. Maybe you are getting on track an extension of what you want to demonstrate. BTW, can you give me the link to the proof? Your link lead me to some great mathematical video's, but I was unable to find the video you are refering to. $\endgroup$ – Fred Simons Dec 19 '17 at 10:06
  • $\begingroup$ For generating the random points, you could just use RandomPoint[Circle[], 9000] or RandomPoint[Disk[], 9000] (depending on what you need). $\endgroup$ – Sjoerd Smit Dec 19 '17 at 19:06
  • $\begingroup$ the proof is outlined in the linked video at about 5:45 $\endgroup$ – Joe Dec 23 '17 at 21:23
  • $\begingroup$ Repeating the proof on the video, we can find a generalization. Consider a probability on the plane that is symmetric with respect to the origin and consider a random triangle ABC. Draw the lines AA' and BB' such that the origin is the center, so A' is opposite to A and B' is opposite to B. Then each of the four triangles CAB, CA'B, CAB' and CA'B' have equal probability and only one of them contains the origin. Therefore, the probability that such a triangle contains the origin is 1/4. That is what has been tested in the question. $\endgroup$ – Fred Simons Jan 7 '18 at 13:32
15
$\begingroup$

Some issues with your code are that the list of random data you generate is not a PackedArray; moreover creating random numbers is more efficient if you generate all at once like angles = RandomReal[2 Pi, {n, 3}]. Even if the data were packed, it would be unpacked when mapping Polygon to it. Moreover, testing against triangles with RegionMember is like "shooting sparrows with a cannon". In general: Map on big data is not so fast, even if it is considerably better than using For loops. Best efficiency is obtained with CompiledFunctions with the option RuntimeAttributes -> {Listable}.

The following CompiledFunction computes the signed areas spanned by the origin and each pair of corners of the input triangle and returns 1 if all the signs of these signed areas are equal (that's only possible if the origin lies within the triangle).

Quiet[Block[{P, PP},
  PP = Table[P[[i, j]], {i, 1, 3}, {j, 1, 2}] /. Part -> Compile`GetElement;
  With[{
    det3 = Det[{PP[[1]], PP[[2]]}],
    det1 = Det[{PP[[2]], PP[[3]]}],
    det2 = Det[{PP[[3]], PP[[1]]}]},
   testTriangle = Compile[{{P, _Real, 2}},
     Block[{s1, s2, s3},
      s1 = Sign[det1];
      s2 = Sign[det2];
      s3 = Sign[det3];
      Boole[s1 == s2 && s2 == s3]],
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"
     ]
   ]
  ]]

On my machine, this can test about 25000000 triangles per second. That's about 10000 times faster than the original code. Note that generating the random triangles is more expensive than checking them.

n = 25000000;
triangles = With[{angles = RandomReal[2 Pi, {n, 3, 1}]},
   Join[Cos[angles], Sin[angles], 3]
   ];
N[Total[testTriangle[triangles]]/n] // AbsoluteTiming

{1.01982, 0.249998}

So, the experimental probability is pretty close to 25 %.

$\endgroup$
6
$\begingroup$

This is not efficient but I post it as another way and also for the visualization.

A query function:

qf = Compile[{{p, _Real}, {q, _Real}}, Module[{v = {p, q}, j, k, b},
    {j, k} = Sort@v;
    If[j > Pi, 0, Boole[Pi < k < j + Pi]]], CompilationTarget -> "C", 
   Parallelization -> True, RuntimeOptions -> "Speed"];

Simulation:

p[v_] := Module[{s = {#2 - #1, #3 - #1} & @@ (Sort@v)},
   qf @@ s];
nu[n_] := With[{rnd = RandomReal[2 Pi, {n, 3}]}, N[Total[p /@ rnd]/n]]

Visualization functions:

nug[v_] := 
 With[{res = 
    p@v}, {Graphics[{Circle[], PointSize[0.04], Point[{0, 0}], 
     RGBColor@RotateRight[{1, 0, 0}, res], 
     Polygon[{Cos@#, Sin@#} & /@ v]}], res}]
nugf[n_] := 
 With[{y = Transpose[nug /@ RandomReal[2 Pi, {n^2, 3}]]}, 
  Column[{N[Total@y[[2]]/n^2], 
    GraphicsGrid[Partition[y[[1]], n], Spacings -> {0, 0}, 
     ImageSize -> 800]}, Alignment -> Center]]

Some simulations:

Increasing sample size enter image description here

After a long time (minutes) on my old machine:

enter image description here

Some visualizations:

  • sample of 25 triangles enter image description here

  • sample of 100 triangles enter image description here

  • sample of 1600 triangles enter image description here

$\endgroup$
6
$\begingroup$

This is a further extension of @Henrik's answer, a little bit too long for a comment.

When playing with Henric's great solution, I observed that the two settings Parallelization->True and Parallelization->False did not make the difference I hoped for. My parallel kernel configuration shows that automatically 4 local kernels will be launched, whereas the difference between two two settings was only a factor 1.5. So I tried to do the parallelization myself.

The next function is the function testTriangle of Henric, but with Parallelization set to False.

origininside=
  With[{Part=Compile`GetElement},
    Compile[{{tr, _Real,2}},
      Module[{s1, s2, s3},
        s1=Sign[tr[[2,1]]tr[[3,2]]-tr[[2,2]]tr[[3,1]]];
        s2=Sign[tr[[3,1]]tr[[1,2]]-tr[[3,2]]tr[[1,1]]];
        s3=Sign[tr[[1, 1]]tr[[2,2]]-tr[[1,2]]tr[[2,1]]];
        Boole[s1==s2==s3]],
      CompilationTarget->"C",
      Parallelization->False,
      RuntimeAttributes->{Listable}]
    ];

On my quadcore computer, I launch 8 kernels:

LaunchKernels[8];

I seem to remember that a C-compiled function is not immediately available in the subkernels, so I expected that now I had to do the compilation of the function origininside on each of the subkernels. Maybe that idea dates from a long time ago, maybe I am wrong, but I was pleasantly surprised to see that, having compiled the function in the main kernel, it is immediately avaible on the subkernels as well.

Therefore, with the following command I create and test 25,000,000 triangles:

nn=2500;
  N@Total[ParallelTable[
    Total[origininside @  With[{angles=RandomReal[2 Pi,{nn,3,1}]},
      Join[Cos[angles],Sin[angles],3]] ], {10000}]]/10000/nn // AbsoluteTiming

{0.770372, 0.250022}

Since this timing includes the time required for the construction of the triangles, this construction of the parallelization makes the computation about 2.5 times as fast.

Addendum

This addendum contains some mathematical refinements, including a proof that the probability indeed is 1/4.

The first refinement has to do with the construction of the triangles. There is no loss in generality when we assume the first pont of the triangle is at the top of the circle: A={0,1}. We then choose only two random points B={b1,b2} and C={c1,c2} on the unit circle. The condition that the origin is inside the triangle then amounts to b1 and c1 having opposite sign and the intersection of the line BC with the vertical axis is negative. The latter gives:

Assuming[b1 c1<0, Simplify[(z /. Solve[ (1-\[Lambda]){b1,b2}+\[Lambda]{c1,c2}=={0,z}, {\[Lambda], z}][[1]]) <0]]

(b1-c1) (-b2 c1+b1 c2)<0

which can be further simplified to

b1(-b2 c1+b1 c2)<0

The function origininside is now simpler:

origininside=With[{Part=Compile`GetElement},
  Compile[{{tr, _Real, 2}},
    If[tr[[1,1]]tr[[2,1]]>0, 0,
      Boole[tr[[1,1]](-tr[[1,2]] tr[[2,1]]+tr[[1,1]] tr[[2,2]])<0]],
    CompilationTarget->"C",
    Parallelization->False,
    RuntimeAttributes->{Listable}]]

The testing on 8 subkernels:

LaunchKernels[8]

(In the following command I use Henric's construction for random points on the circle; this turns out to be faster than using RandomPoint as mentioned by @Sjoerd Smit.)

nn=2500;
N@Total[ParallelTable[
  Total[origininside@With[{angles=RandomReal[2 Pi,{nn,2,1}]},
    Join[Cos[angles],Sin[angles],3]]],{10000}]]/10000/nn // AbsoluteTiming

{0.546887,0.25012}

This is about 25% faster.

The fact that the probability that the origin is inside a random triangle on the unit circle is 1/4 can be seen in the following way.

(Draw by hand a circle with a point A at the top and a point B at the right hand side. Let A' and B' be the opposite points of A and B. I would have placed this picture here in this answer, but I have no idea how to do that.)

The first point A is the top of the circle. Then we choose randomly the point B on the circle. Without loss of generality we may assume that B is at the right hand side of the circle, as in in the figure.

Now we choose the point C randomly on the unit circle. For the origin being inside the trangle, the point C must be on the arc A'B', This arc is as long as the arc AB.

Otherwise stated: we choose the point x randomly in the interval [0,π]. Then we choose y randomly in the interval[0,2π] and the criterion is that y in a subinterval of length x.

Therefore, the problem of finding the probability that the origin is inside a random triangle on the unit circle is equivalent with the following one:

What is the probability that when we randomly choose a number x in the interval [0,1], the random number yin the interval [0,1] will be less than x/2?

Obviously, the answer is

Integrate[x/2, {x, 0, 1}]

1/4

For the testing, we construct a list of 250000001 random numbers, so that we have 25000000 pairs. Then we can proceed in a lot ways, for example by partitioning the list into pairs with an increment of 1, and then count how many pairs have the property that the first element is larger than two times the second argument.

Here are two ways:

nn = 25000000; lst = RandomReal[1, {nn + 1}];

N[Total[UnitStep[ListCorrelate[{1, -2}, lst]]]/nn] // AbsoluteTiming

{0.45397, 0.250025}

lst1 = Most[lst];
Rest[lst];
N[Total[UnitStep[lst1 - 2 lst2]]/nn ] // AbsoluteTiming

{0.20445, 0.250025}

This is indeed very fast, without having used any C-compilation or parallel computing.

$\endgroup$
  • $\begingroup$ That's a very good observation! Indeed, I am quite surprised as my experiences with ParallelTable and its relatives have been rather disillusioning... $\endgroup$ – Henrik Schumacher Dec 26 '17 at 13:12
1
$\begingroup$
cf = With[{code = 
    Sign@Det@{{x1,y1},{x2,y2}}==Sign@Det@{{x2,y2},{x3,y3}}==Sign@Det@{{x3,y3},{x1,y1}}//Boole
    }, 
   Compile[{{A, _Real, 1}}, 
    Block[{x1 = A[[1]], x2 = A[[2]], x3 = A[[3]], y1 = A[[4]], y2 = A[[5]], y3 = A[[6]]},
     code
     ], RuntimeAttributes -> {Listable}
    ]
   ];

n = 10^7;

With[{angles = RandomReal[2 Pi, {n, 3}]},
  Join[Cos[angles], Sin[angles], 2]] // cf // Mean // N // AbsoluteTiming

{2.04003, 0.249968}

Module[{n, angles, x1, x2, x3, y1, y2, y3, s1, s2, s3},
    n = 10^7;
    angles = RandomReal[2 Pi, {3, n}];
    {x1, x2, x3} = Cos[angles];
    {y1, y2, y3} = Sin[angles];
    s1 = Sign[x1 y2 - x2 y1];
    s2 = Sign[x2 y3 - x3 y2];
    s3 = Sign[x3 y1 - x1 y3];
    Unitize[s1 - s2]~BitOr~Unitize[s2 - s3]
    ] // Mean // 1.0 - # & // AbsoluteTiming

{1.44645, 0.250057}

$\endgroup$
1
$\begingroup$

Since I wanted to try this out on a computer with very limited memory, generating the whole list of triangles all at once was out of the question. Instead, I use a Do[] loop to accumulate a running count of the triangles that passed the criterion, and then divide by the total number of samples to estimate the probability. Thus:

BlockRandom[SeedRandom[42, Method -> "MersenneTwister"]; (* for reproducibility *)
            n = 1*^6; (* number of samples *) s = 0;
            Do[s += Boole[Graphics`PolygonUtils`InPolygonQ[
                          Transpose[Through[{Cos, Sin}[RandomReal[2 π, 3]]]], {0, 0}]],
               {n}];
            N[s/n, 20]]
   0.25003500000000000000

where I had used the undocumented point-in-polygon test featured in this answer.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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