Commit 9fcbbe36 authored by Mattias Fält's avatar Mattias Fält
Browse files

Avoid extra allocation in sobol generation

parent e6c4b901
......@@ -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
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment