8
$\begingroup$

Introduction


I am now trying to program an algorithm to achieve an uniform distribution of points within a space of dimension n. I started working in 2D.

The idea is to achieve an uniform distribution of points with a heuristic method.The important aim is to get a fast algorithm but not necessarily an extremely precise one.

I would like to have a better uniform distribution of points than a statistical method could provided. In fact, the idea to use a heuristic to solve this problem was suggested to me by a sphere packing problem.

I first tried a heuristic method based on the electrostatic repulsion principle. Without success.

The new method I am trying now is the coupling of a Delaunay triangulation with the calculation of the centroid associated polygons.


About the method


After more work on the method, I can say that here:

  • The speed of execution depends largely on the algorithm responsible for the Delaunay triangulation.

  • The precision and homogeneity in the distribution depends largely on the number of available point on the border of the study area.

I leave here the code that I developed, if it may be useful to anyone interested in the project. The code has not been optimized.



Informations


  • The dimension of my 2D space are :
    x ∈ [0,1]
    y ∈ [0,1]
  • The list of point ends with the convex hull of the search space. For my example, the end of the list are {....,{0,0},{0,1},{1,0},{1,1}}.



Code


Here's the code:

(* Fonction n1 *)
Cut[x_, y_] :=
  Block[
        {step1, step2, Res},
         step1 = (Abs[Subtract @@ x])/y;
         step2 = Range[x[[1]], x[[2]], step1];
         Res = Partition[step2, 2, 1]
       ];

(* Fonction n2 *)
CutN[x_, y_] :=
  Block[
        {Res},
         Res = Cut[#, y] & /@ x
       ];

(* Fonction n3 *)
CutB[x_, y_] :=
  Block[
        {Res},
         Res = Partition[MapThread[## &, CutN[x, y]], Length@x]
       ];

(* Fonction n4 *)
SpaceDistrib[x_, y_, z_] :=
  Block[
        {step0, step1, step2, step3, step4, Res},
         step0 = Tuples[#] & /@ CutB[x, y];
         step1 = RandomVariate[UniformDistribution[#], Abs[z]] & /@ # & /@ 
         CutB[x, y];
         step2 = MapThread[## &, #] & /@ step1;
         Res = Partition[#, Length@x] & /@ step2
       ];

(* Fonction n5 *)
SpaceDistribLimit[x_, y_, z_] :=
  Block[
        {step0, step1, step2, step3, step4, Res},
         step0 = Tuples[#] & /@ CutB[x, y];
         step1 = RandomVariate[UniformDistribution[#], 
         Abs[z - Length@step0[[1]]]] & /@ # & /@ CutB[x, y];
         step2 = MapThread[## &, #] & /@ step1;
         step3 = Partition[#, Length@x] & /@ step2;
         Res = MapThread[Flatten[{#1, #2}, 1] &, {step3, step0}]
      ];

(* Fonction n6 *)
Barycentre[x_] :=
  Block[
        {step1, step2, Res},
         step1 = Transpose@x;
         step2 = Total /@ step1;
         Res = step2/(Length@x)
       ];

(* Fonction n7 *)
PaireofD[x_] := 
   Block[
         {Paire, step1},
          Paire[y_] := {y[[1]], #1} & /@ y[[2]];
          step1 = Paire[#] & /@ x
        ];

(* Fonction n8 *)
TriangleD[x_] :=
  Block[
        {step0, step1, step2, step3},
         step0 = Flatten[PaireofD[x], 1];
         step1 = (Flatten[#] & /@ Subsets[step0, {3}]);
         step2 = Table[Gather[step1[[i]], #1 == #2 &], {i, 1, Length@step1, 1}];
         step3 = Union@(Union /@ 
       Flatten /@ Cases[step2, {{k1_, k1_}, {k2_, k2_}, {k3_, k3_}}])
       ];

(* Fonction n9 *)
BarycentreT[x_, y_] :=
  Block[
        {step0, step1, step2},
         step0 = TriangleD[x];
         step1 = y[[#]] & /@ step0;
         step2 = Barycentre[#] & /@ step1
       ];

(* Fonction n10 *)
PointToMove[x_, y_] :=
  Block[
        {step1, step2, step3, step4},
         step1 = x;
         step2 = Min[Position[step1[[All, 1]], #] & /@ (ConvexHull@y)] - 1;
         step3 = {y[[1 ;; step2]], Range[1, step2, 1]}
       ];

(* Fonction n11 *)
PolyBary[x_, y_] :=
  Block[
        {step1, step2, step3, step4},
         step1 = x;
         step2 = Min[Position[step1[[All, 1]], #] & /@ (ConvexHull@y)] - 1;
         step3 = step1[[1 ;; step2, 2]];
         step4 = y[[#]] & /@ step3
       ];

(* Fonction n12 *)
Method1[x_, y_] :=
  Block[
        {step1, step2, step3, step4, Res1, step5, step6, step7,
     step8, Res},
         step1 = PolyBary[x, y];
         step2 = Barycentre[#] & /@ step1;
         step3 = PointToMove[x, y][[1]];
         step4 = MapThread[{#1, #2} &, {step3, step2}];
         Res = {step2, step4, EuclideanDistance @@@ step4}
       ];

(* Fonction n15 *)
Stopcrit[x_, y_] :=
  Block[
        {step1, step2, step3},
         step1 = Max@x;
         step2 = If[step1 < y, True, False]
            ];


GIT1[x_] :=
  Block[
        {stop, ToIt, Stock1, Stock2, n, step1, step2, step3, 
    step4, step5, Res},
         stop = False;
         ToIt = x;
         Stock1 = {ToIt};
         Stock2 = {};
         n = 1;

         While[
               stop == False,
               step1 = DelaunayTriangulation[ToIt];
               step2 = Method1[step1, ToIt];
               step3 = Min[
       Position[step1[[All, 1]], #] & /@ (ConvexHull@ToIt)] - 1;
               ToIt = Join[step2[[1]], ToIt[[step3 + 1 ;; -1]]];
               stop = Stopcrit[step2[[3]], 0.01];
               Stock1 = Insert[Stock1, ToIt, -1];
               Stock2 = Insert[Stock2, step2[[2]], -1];
               n++
              ];
          Res = {Stock1, Stock2, n}
       ];


(* Variables de test *)
esp1 = {0, 1};
esp2 = {{0, 1}, {0, 1}};
esp3 = {{0, 1}, {0, 1}, {0, 1}};
l1 = {{0, 0}, {0, 1}, {1, 0}, {1, 1}};
P1 = Flatten[SpaceDistribLimit[esp2, 1, 30], 1];
P2 = Join[P1, {{0, 0.5}, {0.5, 1}, {1, 0.5}, {0.5, 0}}];
D1 = DelaunayTriangulation[P1];

(*********************************************)
(* Fonction n *)
Cut[esp1, 2];
Cut[esp1, 3];
(* Fonction n2 *)
CutN[esp2, 2];
CutN[esp2, 3];
(* Fonction n3 *)
CutB[esp2, 2];
(* Fonction n4 *)
SpaceDistrib[esp2, 1, 2];
SpaceDistrib[esp3, 3, 2];
(* Fonction n5 *)
SpaceDistribLimit[esp3, 1, 10];
SpaceDistribLimit[esp3, 3, 10];
(* Fonction n6 *)
Barycentre[l1];
(* Fonction n7 *)
PaireofD[D1];
(* Fonction n8 *)
TriangleD[D1];
(* Fonction n9 *)
BarycentreT[D1, P1];
(* Fonction n10 *)
PointToMove[D1, P1];
(* Fonction n11*)
PolyBary[D1, P1];
(* Fonction n12 *)
Method1[D1, P1][[1]];(* point de sortie *)
Method1[D1, P1][[2]];(* flèche représentative du mouvement *)
Method1[D1, P1][[3]];(* Norme des déplacements *)
(* Fonction n14 *)
Stopcrit[Method1[D1, P1][[3]], 5];
(* ************* *)
GIT1[P1][[1]];
GIT1[P1][[2]];
GIT1[P1][[3]];



Result


I obtained the following results:
Original distribution 1 :

enter image description here

Result 1 :

enter image description here

Original distribution 2 :

enter image description here

Result 2 :


enter image description here


Towards an alternative resolution.


The stated goal of this subject is to obtain a series of random points to cover as much as possible an n-dimensional space. The number of points in the distribution may vary. As suggested by Mr. Daniel Lichtblau and Mr. Oleksandr R. a low discrepancy sequence appears to be a viable solution. After some research, several names have emerged :
- Sobol's sequence.
- Van Der Corput's sequence.
- Halton's sequence.
- Hammersley's sequence....


Questions


What is the best sequence ?

I am interested in all pseudo-code or documentation to generate low discrepancy sequence in n dimensions.

$\endgroup$
  • $\begingroup$ Despite the fact documentation tells of 1D, 2D, 3D and nD Delaunay tetrahedralizations under DelaunayMesh, only dimensions 1 through 3 would seem to work. $\endgroup$ – kirma May 18 '15 at 16:03
  • 1
    $\begingroup$ Quick and dirty but seems to give an okay distribution over the unit hypercube: quasiUniformPoints[n_, dim_] := Mod[Outer[Times, N@Range[n], E^Range[dim]], 1]. Try ListPlot in 2D or ListPointPlot3D in 3D to see what it gives. $\endgroup$ – Daniel Lichtblau May 18 '15 at 18:20
  • 1
    $\begingroup$ A grid will be uniform but not random. A random arrangement will not be uniform. There are special distributions used for Monte Carlo sampling that are uniform but not grid-like, and perhaps this is what you want. Could you please clarify your requirements of the resulting distribution? $\endgroup$ – Oleksandr R. May 19 '15 at 13:08
  • 1
    $\begingroup$ I should mention that the method I showed is a (not very high quality) version of a low discrepancy sequence. That's is the sort of thing one uses in quasi-Monte Carlo quadrature (per @Oleksandr R's remark). $\endgroup$ – Daniel Lichtblau May 19 '15 at 14:56
  • 1
    $\begingroup$ What I coded is a not very high quality rendition of what's called a "low discrepancy sequence". So that would be the thing to look up if you want to get more detail or perhaps better quality variants. $\endgroup$ – Daniel Lichtblau May 23 '15 at 22:15
5
$\begingroup$

I post here my own solution.
I post here all the code I have written to generate low discrepancy sequences.
The algorithms presented here are not necessarily optimized.
Feel free to send me any corrections or improvements.

Do not hesitate to post here also useful code to generate other low discrepancy sequence.


Generation of a low discrepancy sequence



"Translation irrationnelle du tore"

RandSuiteTore[x_, y_] :=
  Block[
        {n, dim, Stock, Stock1, step1, step2, step3, step4, Res},

        (* Initialisation  *)
         n = 1;
         dim = y;
         Stock = {};
         Stock1 = {};
        (* Boucle *)
         While[
                n <= dim,
                step1 = RandomReal[{0, 10^2}];
                step2 = NextPrime[step1, 1];
                step3 = Table[                                           
                              N@FractionalPart[i*Sqrt@step2],
                              {i, 1, x, 1}
                             ];
                Stock = Insert[Stock, step3, -1];
                Stock1 = Insert[Stock1, step2, -1];
                n++
              ];
          step4 = Partition[MapThread[## &, Stock], Length@Stock];
          Res = {step4, Stock1}
       ];

SuiteTore[x_, y_, z_] :=
  Block[
        {n, dim, Stock, Stock1, step1, step2, step3, step4, Res},
         (* Initialisation  *)
         n = 1;
         dim = y;
         Stock = {};
         Stock1 = {};
         (* Boucle *)
         While[
               n <= dim,
               step2 = z[[n]];
               step3 = Table[                                       
                             N@FractionalPart[i*Sqrt@step2],
                             {i, 1, x, 1}
                            ];
               Stock = Insert[Stock, step3, -1];
               Stock1 = Insert[Stock1, step2, -1];
               n++
              ];
        Res = Partition[MapThread[## &, Stock], Length@Stock]
       ];


"Sobol's sequence"

Sobol[x_, y_] :=
  Block[
        {step1, step2},
         step1 = ToString@RandomReal[];
         step2 =
                BlockRandom[
                            SeedRandom[
                                       step1, 
                                       Method -> {"MKL", Method -> {"Sobol", "Dimension" -> y}}
                                      ];
                            RandomReal[{0, 1}, {x, y}]
                           ]
               ];


"Niederreiter's sequence"

Nied[x_, y_] :=
  Block[
        {step1, step2},
         step1 = ToString@RandomReal[];
         step2 =
                 BlockRandom[
                             SeedRandom[
                                        step1,  
                                        Method -> {"MKL", Method -> {"Niederreiter", "Dimension" -> y}}
                                       ];

                             RandomReal[{0, 1}, {x, y}]
                            ]
               ];


"Van der Corput's sequence"

VanDerCorput[x_, y_] := Table[
                              FromDigits[{Reverse[IntegerDigits[i, base]], 0}, y],
                              {i, x}
                             ];


"Halton's sequence"

RandHalton[x_, y_] :=
  Block[
        {step1, step2, step3, Res},
         step1 = Table[RandomReal[{0, 10^2}], {y}];
         step2 = NextPrime[step1, 1];
         step3 = Table[
                       VanDerCorput[x, i],
                       {i, step2}
                      ];
         Res = {Partition[MapThread[## &, step3], Length@step3], step2}
       ];

Halton[x_, y_, z_] :=
  Block[
        {step2, step3, Res},
         step2 = z;
         step3 = Table[
                       VanDerCorput[x, i],
                       {i, step2}
                      ];
         Res = Partition[MapThread[## &, step3], Length@step3]
       ];


To test the code


Parameters : [ numberofpoints, dimension ]  

"Translation irrationnelle du tore"

v1 = RandSuiteTore[1000, 2];
ListPlot[v1[[1]]];
v1[[2]];

With :
- v1[[1]]: The low discrepancy sequence
- v1[[2]]: The parameter used (prime number)

To repeat :

SuiteTore[1000, 2, v1[[2]]]  

"Sobol's sequence"

Sobol[1000, 2];

"Niederreiter's sequence"

Nied[1000, 2];

"Van der Corput's sequence"

VanDerCorput[1000, 2]

"Halton's sequence"

RandHalton[1000, 2]


Possibles results


2D

In 2D

3D

enter image description here


Commentary


As we could observe, the "Translation irrationnelle du tore" algorithm generate sequence with variable discrepancy whose values is not necessarily guaranteed unlike conventional sequence such as Sobol's sequence.

If somebody has an algorithm to quickly calculate the discrepancy to compare sequences each other, I'm strongly interested.




| improve this answer | |
$\endgroup$
  • $\begingroup$ I'd be interested in code for 3D low-discrepancy sequences--particularly, for "open-ended" sequences, for which one does not need to decide in advance how many points to generate. $\endgroup$ – Paul B. Slater Apr 12 '17 at 19:16

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.