10
$\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.

Improve this question
$\endgroup$
8
  • $\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, 2015 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, 2015 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, 2015 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, 2015 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, 2015 at 22:15

1 Answer 1

Reset to default
6
$\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$
3
  • $\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, 2017 at 19:16
  • $\begingroup$ In this version, the randomSeed for Sobol doesn't work. ( See here for a simple fix: mathematica.stackexchange.com/questions/34847/… $\endgroup$
    – Martin Roberts
    Jan 2, 2022 at 10:26
  • $\begingroup$ But how to solve the problem that at high dimensions the quasi-random properties practically disappear, and the correlation increases? $\endgroup$
    – dtn
    Apr 4 at 7:46

Your Answer

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

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