naught101 / sobol_seq

python implementation of Sobol' sequence generator
MIT License
67 stars 29 forks source link

i4_sobol_generate ignores skip parameter #14

Closed lbrichards closed 4 years ago

lbrichards commented 6 years ago

in i4_sobol_generate() I think you need to add + skip on line 131 as follows:

        seed = j + 1 + skip

Currently skip has no effect.

naught101 commented 6 years ago

What version of the code are you looking at?

https://github.com/naught101/sobol_seq/blob/6ac1799818a1b9a359a5fc86517173584fe34613/sobol_seq/sobol_seq.py#L127

lbrichards commented 6 years ago

Version 0.1.2 is what I installed with pip.

It contains the following:

def i4_sobol_generate(dim_num, n, skip=1): """ i4_sobol_generate generates a Sobol dataset.

Parameters:
  Input, integer dim_num, the spatial dimension.
  Input, integer N, the number of points to generate.
  Input, integer SKIP, the number of initial points to skip.

  Output, real R(M,N), the points.
"""
r = np.full((n, dim_num), np.nan)
for j in range(n):
    seed = j + 1
    r[j, 0:dim_num], next_seed = i4_sobol(dim_num, seed)

return r

On Jun 4, 2018, at 12:26 PM, naught101 notifications@github.com wrote:

What version of the code are you looking at?

https://github.com/naught101/sobol_seq/blob/6ac1799818a1b9a359a5fc86517173584fe34613/sobol_seq/sobol_seq.py#L127 https://github.com/naught101/sobol_seq/blob/6ac1799818a1b9a359a5fc86517173584fe34613/sobol_seq/sobol_seq.py#L127 — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/naught101/sobol_seq/issues/14#issuecomment-394224808, or mute the thread https://github.com/notifications/unsubscribe-auth/AKgQPIjWuA9dQjqmC6wxU-S9sBE2HRNNks5t5KjxgaJpZM4UX-rT.

naught101 commented 6 years ago

This was fixed at https://github.com/naught101/sobol_seq/issues/8

Install the git version, it includes the fix. @CnrLwlss might know something about updating the pypi version - I don't use pypi.

lbrichards commented 6 years ago

Understood, and thanks for the clarification. My default way of installing dependencies is pip, but I can of course clone in this case. BTW, this is probably not the appropriate place to discuss, and I hope you don't mind my asking, but the i4_sobol function contains a self, as though someone intended it to be a class. Indeed, it seems to be a function wanting to be a class, since unless there is a class instance the self.seed_save will not be remembered from one function call to the next. Someone, it seems, intended for this function to remember its most recent seed upon each call. I presume this was to save computation time when a very large seed is used. I would be happy to propose some refactoring, if you're open to such an idea. Essentially I would propose a rewrite as

class i4_sobol(object):
    def __init__(self, dim_num, seed ):

and add a __call__ method allowing it to retain its function-like behavior.

naught101 commented 6 years ago

I basically uploaded this because I wanted to make some minor changes and make it easier for others to install it. Now, looking through the code, I'm a bit surprised I didn't already do that. I unfortunately don't have time to work on this at the moment, but I'm more than happy to accept pull requests, ideally with backwards compatibility.

It would be especially sensible to move all those globals into class attributes. Not sure where you're seeing a self, but yes, I think that would be a great idea.

The i4_* syntax is a hang-over from Fortran, it could also be removed, although it would be good to leave those functions there as wrappers with deprecation warnings.

lbrichards commented 6 years ago

Kudos to you for taking the first initiative. I would love to contribute to such a worthy cause. I have found several variants of this code, but none very Pythonic, and none -- as far as I can tell -- takes advantage of the Gray code efficiency boost by Antonov and Saleev. For my simulations I need to find long pairs of Sobol sequences that have approximately zero correlation, and this is very time consuming. Not sure how much the Gray code implementation would speed things up, but it's worth a try. And making things less Fortran-like and more Pythonic is also nice. I will aim to send you some pull requests after I escape some deadlines...