From 9fcbbe3656ea031993ca19c983b2e84b02d8cef4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mattias=20F=C3=A4lt?= <mattiasf@control.lth.se>
Date: Tue, 1 Dec 2015 18:00:25 +0100
Subject: [PATCH] Avoid extra allocation in sobol generation
---
src/sobol.jl | 85 +++++++++++++++++++++++++---------------------------
1 file changed, 40 insertions(+), 45 deletions(-)
diff --git a/src/sobol.jl b/src/sobol.jl
index 5ab9290..db027c0 100644
--- a/src/sobol.jl
+++ b/src/sobol.jl
@@ -28,69 +28,64 @@ function sobol(dim::Int64, N::Int64)
X = Array(Float64,(N,dim))
MeM = InitSobol(dim)
- for j = 1:N
- X[j,:], nextseed, MeM = getNewSobolVector(dim, nextseed, MeM)
- end
+ nextseed = getNewSobolArray!(dim, nextseed, MeM, X, N)
return X, nextseed, MeM
end
-function sobol(X, nextseed, MeM)
+function sobol!(X, nextseed, MeM)
(N,dim) = size(X)
- for j = 1:N
- X[j,:], nextseed, MeM = getNewSobolVector(dim, nextseed, MeM)
- end
+ nextseed = getNewSobolArray!(dim, nextseed, MeM, X, N)
- return X, nextseed, MeM
+ return nextseed
end
-function getNewSobolVector(dim, seed, MeM)
-
- seed::Int64 = max(floor(seed),0)
-
- if seed == 0
- l = 1
- MeM.lastq = zeros(1,dim)
- elseif seed == MeM.seed_save + 1
- l = smallest0bit(seed)
- elseif seed <= MeM.seed_save
- MeM.seed_save = 0
- l = 1
- MeM.lastq[1:dim] = 0
+function getNewSobolArray!(dim, seed, MeM, X, N)
+ for j = 1:N
+ seed::Int64 = max(floor(seed),0)
+ if seed == 0
+ l = 1
+ MeM.lastq = zeros(1,dim)
+ elseif seed == MeM.seed_save + 1
+ l = smallest0bit(seed)
+ elseif seed <= MeM.seed_save
+ MeM.seed_save = 0
+ l = 1
+ MeM.lastq[1:dim] = 0
- for seed_temp = MeM.seed_save:seed-1
- l = smallest0bit(seed_temp)
- for i = 1:dim
- MeM.lastq[i] = MeM.lastq[i] $ MeM.v[i,l] # bitxor
+ for seed_temp = MeM.seed_save:seed-1
+ l = smallest0bit(seed_temp)
+ for i = 1:dim
+ MeM.lastq[i] = MeM.lastq[i] $ MeM.v[i,l] # bitxor
+ end
end
- end
- l = smallest0bit(seed)
+ l = smallest0bit(seed)
- elseif (MeM.seed_save + 1) < seed
+ elseif (MeM.seed_save + 1) < seed
- for seed_temp = (MeM.seed_save + 1):(seed - 1)
- l = smallest0bit(seed_temp)
- for i = 1:dim
- MeM.lastq[i] = MeM.lastq[i] $ MeM.v[i,l] # bitxor
+ for seed_temp = (MeM.seed_save + 1):(seed - 1)
+ l = smallest0bit(seed_temp)
+ for i = 1:dim
+ MeM.lastq[i] = MeM.lastq[i] $ MeM.v[i,l] # bitxor
+ end
end
- end
- l = smallest0bit(seed)
-
- end
+ l = smallest0bit(seed)
- qrv = Array(Float64,dim)
+ end
- for i = 1 : dim
- qrv[i] = MeM.lastq[i] * MeM.recipd;
- MeM.lastq[i] = MeM.lastq[i] $ MeM.v[i,l] # bitxor
- end
+ #qrv = Array(Float64,dim)
- MeM.seed_save = seed;
- nextseed = seed + 1;
+ for i = 1 : dim
+ X[j,i] = MeM.lastq[i] * MeM.recipd; #qrv[i] = MeM.lastq[i] * MeM.recipd;
+ MeM.lastq[i] = MeM.lastq[i] $ MeM.v[i,l] # bitxor
+ end
- return qrv, nextseed, MeM
+ MeM.seed_save = seed;
+ seed = seed + 1;
+ end
+ return seed
end
function smallest0bit(b)
@@ -211,4 +206,4 @@ function test_sobol()
plot(X[:,1],X[:,2],".m")
display(pp)
-end
\ No newline at end of file
+end
--
GitLab