halton sequence added

Author: mrtkp9993

src/NumericalAlgorithms.jl

@@ -23,7 +23,7 @@ export primes_1000
 export FindRoot1D, FindRootND
 export Dual
 export CalcSingleIntegral, CalcDoubleIntegral
-export LFG, RANMAR, vanderCorputSeq
+export LFG, RANMAR, vanderCorputSeq, haltonSeq
 export RunsTest
 
 end # module

src/Random.jl

@@ -79,13 +79,19 @@ function Base.rand(rng::RANMAR)
 end
 
 "van der Corput sequence, simplest one-dimensional low-discrepancy sequence over the unit interval."
-function vanderCorputSeq(n, b)
-    bn = 0
+function vanderCorputSeq(i, b)
+    bi = 0
     j = 0
-    while n != 0
-        bn += mod(n, b) / (b^(j + 1))
-        n = floor(n / b)
+    while i != 0
+        bi += mod(i, b) / (b^(j + 1))
+        i = floor(i / b)
         j += 1
     end
-    bn
+    bi
 end
+
+"Halton sequence"
+function haltonSeq(i, m)
+    hseq = [vanderCorputSeq(i, primes_1000[j]) for j = 1:m]
+    hseq
+end

test/runtests.jl

@@ -121,9 +121,13 @@ end
     seq1 = [vanderCorputSeq(i, 2) for i = 1:9]
     expected_seq1 = [0.5, 0.25, 0.75, 0.125, 0.625, 0.375, 0.875, 0.0625, 0.5625]
     seq2 = [vanderCorputSeq(i, 3) for i = 1:9]
-    expected_seq2 = [1 / 3, 2 / 3, 1 / 9, 4 / 9, 7 / 9, 2 / 9, 5 / 9, 8 / 9, 1 / 27]
+    expected_seq2 = [0.33, 0.66, 0.11, 0.44, 0.77, 0.22, 0.55, 0.88, 0.03]
     @test seq1 == expected_seq1
-    @test seq2 ≈  expected_seq2 atol = 1e-3
+    @test seq2 ≈  expected_seq2 atol = 1e-1
+
+    hseq1 = haltonSeq(3, 4)
+    expected_hseq1 = [0.75, 0.11, 0.6, 0.42]
+    @test hseq1 ≈ expected_hseq1 atol = 1e-2
 end
 
 @testset "Runs test for randomness" begin