Setting a seed to value of getSeed() doesn't produce expected draws

[deruncie]

This is similar to Issue #7 , but I don't think quite the same:

My goal is to make a MCMC sampler reproducible even if it is stopped in the middle and then re-started.

As an example, I tried drawing 4 normal values, but stopping after two, getting the current seed, re-applying it to zsetseed(), and then drawing two more values. But I get a different output than drawing all 4 at once. Here is an example:

I can set the Ziggurat seed like this:

> zsetseed(bitXor(123456789,3))
> zgetseed()
[1] 3

Draw 4 normal values:

> zrnorm(4)
[1] -0.07785444  1.20543968  0.53777089  1.44539251
> zgetseed()
[1] 1909268553

Try doing this in two steps:
1 draw 2 values
2 get current seed
3 reset zsetseed() to saved seed
4 draw 2 more values

> zsetseed(bitXor(123456789,3))
> zrnorm(2)
[1] -0.07785444  1.20543968
> zgetseed()
[1] 201886211
> zsetseed(bitXor(123456789,201886211))
> zrnorm(2)
[1] -0.7090703  2.3388461
> zgetseed()
[1] 1909268553

As you can see, even though the seed as returned by zgetseed() is the same at the end, the two draws after re-starting are different than when I drew all 4 at once. In baseR, I can save .Random.seed and then reset it and restart where I left off.

Is this possible to do? Or is there a work-around?
I can't find in the code why this doesn't work, but I'm not very experienced at reading C++.

Thanks.

[eddelbuettel]

I can reproduce this. It may have to do with the bitshift you tried to undo with the bitXor. So I am not sure this is entirely achievable (to reset "in stream"), and I had pretty much taken this as a feature of Ziggurat.

For your MCMC use maybe some of the newer stream RNGs may be better. Ziggurat is good, and simple and maybe too simple for your use here.

[eddelbuettel]

Maybe you can get it to work by undoing the bitmap with unsigned long v = 123456789;. We may be loosing something when converting from R as R only has (signed) int32 and double.

Where is this bitXor() coming from?

[deruncie]

bitXor() is from the bitops package. I started with the base bitwXor but that fails with big integers

[eddelbuettel]

Ok. I just played with that the C++ level and then checked with bitops::bitXor() and they do the same :-/

So mid-stream reset may just not be a thing. Unless ... you hack it locally and add another accessor that literally just passes jsr down (without bit-mod'ing it). Should take you no more than 5 minutes.

[deruncie]

OK. Thanks for the help. I'll look up the stream RNG and try the jsr thing. The zrnormR is still ~4x faster than rnorm, so I may just stick to that.

[eddelbuettel]

Maybe 10 minutes. I guess you need a new member function fir ZiggurartLZLLV, maybe

  void resetSeedDirectly(const uint32_t jsrseed) {
         jsr = jsrseed;
  }

and then hook into that similar to the seed getter/setter:

// [[Rcpp::export]]
void zresetseed(unsigned long int s) {
    zigg.resetSeedDirectly(s);
    return;
}

That is clear 'there be dragons' territory as we no long guarantee the other values are good, but the usual caveats it may work.

[deruncie]

I think this (almost) works for me:

I think the trick is that z, w and jcong also have to be set:

I add this to Ziggurart.h

        Rcpp::IntegerVector getpars() {
            Rcpp::IntegerVector pars(4);
            pars[0] = jsr;
            pars[1] = z;
            pars[2] = w;
            pars[3] = jcong;
            return pars; 
        }
        void setPars(Rcpp::IntegerVector pars) {
            jsr = pars[0];
            z = pars[1];
            w = pars[2];
            jcong = pars[3];
        }

and the this to ziggurat.cpp:

// [[Rcpp::export]]
Rcpp::IntegerVector zgetpars() {
    return zigg.getpars();
}
// [[Rcpp::export]]
void zsetpars(Rcpp::IntegerVector pars) {
    zigg.setPars(pars);
    return;
}

Though I think the IntegerVector isn't right because the int32 isn't long enough?

[eddelbuettel]

But if you are instantiated then they should already have a value? The Ziggurat object persists for the session, not the call. Sure we need them?

(Types are tough. R has bit64 which offers integer64 but cleverly resuing the 64 bits of a double. We have no unsigned. You could just static_cast<>() in/from double in a NumericVector. But then again I do not know we really need this.)

[eddelbuettel]

I take that back.

I think you may be right -- I keep looking at RNOR but in the 'fixup' case for the tail the others matter.

Did you test it yet?

[deruncie]

It's that these other values also change as you move through the chain:

values are: jsv, z, w, jcong:

> pars = zgetpars()
> zrnorm(1);zgetpars()
[1] 2.101104
[1] 1925554971  713322877  618681873 1549783358
> zrnorm(1);zgetpars()
[1] 1.827481
[1]  309504328 1074071241  396603440 2813037357
> zrnorm(1);zgetpars()
[1] -1.388457
[1] 2888488517   64231542  811878051 2242876048
> zrnorm(1);zgetpars()
[1] -1.213473
[1]  715369396  231500858  325506388 2326561751
> zsetpars(pars)
> zrnorm(4)
[1]  2.101104  1.827481 -1.388457 -1.213473

Is this code right (modified from above):

Ziggurart.h:

std::vector<uint32_t> getPars() {
            std::vector<uint32_t> pars;
            pars.push_back(jsr);
            pars.push_back(z);
            pars.push_back(w);
            pars.push_back(jcong);
            return pars; 
        }
        void setPars(std::vector<uint32_t> pars) {
            jsr = pars[0];
            z = pars[1];
            w = pars[2];
            jcong = pars[3];
        }

ziggurat.cpp:

// [[Rcpp::export]]
Rcpp::NumericVector zgetpars() {
    return Rcpp::wrap(zigg.getPars());
}
// [[Rcpp::export]]
void zsetpars(Rcpp::NumericVector pars_) {
    std::vector<uint32_t> pars;
    for(int i = 0; i < 4; i++){
        uint32_t par_i = static_cast<uint32_t>(pars_[i]);
        pars.push_back(par_i);
    }
    zigg.setPars(pars);
    return;
}

[eddelbuettel]

Yes, that works:

R> library(RcppZiggurat)
R> zsetseed(123)
R> zrnorm(4)     # four reference values
[1] -0.460512  2.456883  0.117593  0.560411
R> zrnorm(4)     # four new values
[1] 1.1066074 0.0959253 0.5994839 1.7240729
R> zsetseed(123)
R> zrnorm(4)     # back to reference
[1] -0.460512  2.456883  0.117593  0.560411
R> zsetseed(123)
R> zrnorm(2)     # first two
[1] -0.460512  2.456883
R> P <- zgetpars()   # store params
R> zrnorm(2)   # next six as expected
[1] 0.117593 0.560411
R> zrnorm(2)
[1] 1.1066074 0.0959253
R> zrnorm(2)
[1] 0.599484 1.724073
R> zsetpars(P)
R> zrnorm(2)    # as expected too
[1] 0.117593 0.560411
R> 

I'll push a new branch in a moment.