bugGNU Octave - Bugs: bug #54342, rand() produces different results...

 
 

bug #54342: rand() produces different results on octave 4.4.0 compared to earlier versions

Submitter:  Ian McCallion <ianmcc>
Submitted:  Fri 20 Jul 2018 03:29:54 PM UTC
   
 
Category:  Octave Function Severity:  2 - Minor
Priority:  3 - Low Item Group:  Regression
Status:  Confirmed Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * 4.4.0 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 08 Jan 2020 04:02:45 PM UTC, comment #21: 

This is probably worth a longer discussion on the Octave Maintainer's list.  It is not just the floor/round difference with Matlab.  The way Octave calculates the seed, given an input like 42, is also different.  There is another bug report about that somewhere.  And I fixed a large bug in the random number generator when creating numbers of class "single", but this has the ability to throw off the "double" random number generator so that it doesn't create the same numbers as in previous versions of Octave.

Rik <rik5>
Group administrator
Wed 08 Jan 2020 01:25:34 PM UTC, comment #20: 

Matlab and Octave seem to use a different rounding strategy with non integer positive seeds: Octave using round() (as a consequence of casting to uint32?) and Matlab using floor().


octave:1> rand('state',42)
octave:2> rand
ans = 0.6394
octave:3> rand('state',42.4)
octave:4> rand
ans = 0.6394
octave:5> rand('state',42.6)
octave:6> rand
ans = 0.038552
octave:7> rand('state',43)
octave:8> rand
ans = 0.038552



>> rand('state',42)
>> rand
ans =
    0.1591
>> rand('state',42.4)
>> rand
ans =
    0.1591
>> rand('state',42.6)
>> rand
ans =
    0.1591
>> rand('state',43)
>> rand
ans =
    0.3448


Guillaume <gyom>
Sun 22 Jul 2018 10:52:50 PM UTC, comment #19: 

Yes, muddled.  That's why I wouldn't be worried about compatibility with its v5 syntax.

The complicated part--and I'm not saying that this is a issue here, for the reason I just stated--is that I thought I had read that Matlab's recent behavior is to have just one r.n.g. for all such functions.  Conflating compatibility across the board will mean


octave:6> rand("state",0)
octave:7> randn("state",0)
octave:8> x1 = rand()
x1 =  0.84442
octave:9> x2 = rand()
x2 =  0.75795
octave:10> rand("state",0)
octave:11> x1 = rand()
x1 =  0.84442
octave:12> y1 = randn()
y1 = -1.2248
octave:13> x2 = rand()
x2 =  0.75795


isn't consistent; the calling of y1=randn() should alter the next outcome of x2=rand().

Dan Sebald <sebald>
Sun 22 Jul 2018 10:32:53 PM UTC, comment #18: 

Matlab was very muddled in their transition to new random number generators, and presumably that confusion will be forced on to Octave as well.

As I understand it, our "state" is the equivalent of their "twister" and IS the recommended random number generator to use.


Rik <rik5>
Group administrator
Sun 22 Jul 2018 07:51:49 PM UTC, comment #17: 

I'm not sure it is worthwhile trying to aim for compatibility with Matlab v5, as the random number generators for that version apparently aren't statistically in line with randomness and Mathworks suggests not using that anymore.  So maybe it is best to just make Octave 4.4.1 rand() "seed" and "state" consistent with Octave 4.2.2 rand() and consider it deprecated at that point.  (Add a note to rand() documentation that "seed"/"state" will be deprecated and things will be considerably different in the future, not only in syntax but also default seed.)  Hopefully there aren't too many users who've gotten used to the Octave 4.4.0 'state' behavior for rand().

Then in future versions, look towards consistency with all the other numerical programs using the twister code (i.e., the link in Comment #15).  Maybe then add a "seedlegacy" option to whatever routine seeds the future generator that will reset the r.n.g. to the seed used in Octave 4.4.1/4.2.2.  That way, there's a good chance that users can run their legacy code in the future and be able to repeat results of old.

Moving forward will be a bit confusing because of underlying issues of how the generator is used, but this way there'd be a straightforward solution in the short term.

Dan Sebald <sebald>
Sun 22 Jul 2018 07:22:14 PM UTC, comment #16: 

Hi Rik,

Yes indeed, matlab compatibility is a more worthwhile discussion topic than earlier Octave release compatibility. Compatibility can take other forms than imitating Matlab though. E.g. To enable cooperation with a matlab-using co-worker on my project we mimic Octave 4.2.2 semantics everywhere using:

    if (isoctave())
        %% Mimic octave 4.2.2 when using Octave 4.4.0
        if state == Inf
            rand('state',0);
        else
            rand('state',floor(state));
        end
        rvec=rand(size(abin));
    else
        %% Mimic Octave 4.2.2 when using matlab
        mt=octave_twister_seed(state);
        ss=RandStream('mt19937ar');
        ss.set('State', mt);
        rvec=ss.rand(size(abin));
    end


I can probably provide the .m source for octave_twister_seed(state) for you if you think that gets you any further forwards, though as it is not my own work I would need to check out of courtesy first.

Ian

Ian McCallion <ianmcc>
Sun 22 Jul 2018 04:17:37 PM UTC, comment #15: 

The question for me is whether we should purposefully break backward compatibility, and move to Matlab compatibility so that the same random number seed value would produce the same random number streams on both Octave and Matlab.  This could be important for people who want to directly take Matlab code and run it in Octave.

See bug #42557: rand seed for twister algorithm differs from Matlab.

Rik <rik5>
Group administrator
Sun 22 Jul 2018 05:56:27 AM UTC, comment #14: 

Rik,

Thanks. FWIW, my software generates seed states that are all positive and I can confirm that

  •   using floor(seed) instead of seed
  •   using 0 when seed==Inf

restores compatibility between 4.4.0 and 4.2.2, for all the cases I've tried so far.

Regarding "how far to go in preserving continuity", for me at least, there is no need to go further since I now have code that accesses my "old" data. Other users may be unable or unwilling to change their code and will be stuck. So I believe the bug should be fixed. A quick unscientific check suggests that in fact for +ve values Octave 4.4.0 uses round(seed) rather than floor(seed) so from the outside anyway, it looks like it should be an easy and safe fix to apply.

Thanks for your help.

Ian

Ian McCallion <ianmcc>
Sun 22 Jul 2018 05:06:15 AM UTC, comment #13: 
Dmitri A. Sergatskov <dasergatskov>
Sun 22 Jul 2018 04:55:59 AM UTC, comment #12: 

This syntax apparently has been deprecated in Matlab.  It might be difficult finding any v4/v5 documentation, which would probably just state that the v value should be integer.

Dan Sebald <sebald>
Sun 22 Jul 2018 04:15:19 AM UTC, comment #11: 

Certainly looks like Octave versions 4.2 and below used the equivalent of


rand ("state", floor (v)))


But, that is only for positive values of v.  For negative values, they wrapped around and the mapping is


intmax ("uint32") - floor (v);


This is getting complex.  I'm re-opening the bug report, but there should probably be a discussion on the Octave Maintainer's list about how far we want to go in preserving continuity.



Rik <rik5>
Group administrator
Sat 21 Jul 2018 02:09:03 AM UTC, comment #10: 

Ah, yes, you're right.  Doing the conversion prior to the function call isn't quite the same thing.

Dan Sebald <sebald>
Sat 21 Jul 2018 01:03:28 AM UTC, comment #9: 

That case is covered too.  Your example gets caught by the same issue which is that the state vector returned is of uint32_t.  See below.


octave:2> format long
octave:3> v = rand ("state");
octave:4> class (v)
ans = uint32
octave:5> size (v)
ans =

   625     1

octave:6> v(:) = Inf;
octave:7> v(1:3)
ans =

  4294967295
  4294967295
  4294967295

octave:8> rand ("state", v)
octave:9> rand
ans =    3.215697975380883e-01
octave:10> v = Inf (625,1);
octave:11> v(1:3)
ans =

   Inf
   Inf
   Inf

octave:12> rand ("state", v)
octave:13> rand
ans =    1.464507627751816e-01


So if you actually initialize with a state vector of Inf it does return the same as a state vector of 0 or -Inf.

Rik <rik5>
Group administrator
Fri 20 Jul 2018 11:26:51 PM UTC, comment #8: 

Sorry, it took a long time to build, and I was then distracted...  Works nicely.

Should there be any consideration for the other path, i.e., where state V is an actual vector?  I recall that being a different condition in the code.  To illustrate:


octave:13> v = rand ('state');
octave:14> v(:) = 0; rand ('state', v); rand
ans =  0.14645
octave:15> v(:) = -Inf; rand ('state', v); rand
ans =  0.14645
octave:16> v(:) = Inf; rand ('state', v); rand
ans =  0.32157


Dan Sebald <sebald>
Fri 20 Jul 2018 08:10:20 PM UTC, comment #7: 

I pushed a change to the stable branch here (https://hg.savannah.gnu.org/hgweb/octave/rev/c759aa39c23a).

Marking as fixed and closing report.

Rik <rik5>
Group administrator
Fri 20 Jul 2018 06:47:51 PM UTC, comment #6: 

Attached is a small diff that works for me.

(file #44603)

Rik <rik5>
Group administrator
Fri 20 Jul 2018 06:31:49 PM UTC, comment #5: 

@Dan: "state" and "seed" are totally different.  "seed" sets the seed for the older generation of random number generators AND switches Octave to use them when rand() is called.  It should be avoided.

I haven't looked at this, but I think this could be quite easy to resolve.  The incoming array to set internal state is now uint32 type.  But check out the behavior of uint32 for various inputs


octave:8> uint32 (-Inf)
ans = 0
octave:9> uint32 (NaN)
ans = 0
octave:10> uint32 (Inf)
ans = 4294967295


As expected, the constructor for an int maps values which are outside the range to either the min or the max representable value.  My guess is that previously Inf was getting mapped to 0.

Rik <rik5>
Group administrator
Fri 20 Jul 2018 06:23:27 PM UTC, comment #4: 

@Ian: You can change your encoding to use "-Inf" which has the same behavior across versions.  For example,


perl -i.bak -pe 's/rand ("state", Inf)/rand ("state", -Inf)/ig' NAME_OF_SCRIPT



Rik <rik5>
Group administrator
Fri 20 Jul 2018 06:16:20 PM UTC, comment #3: 

The following seems to return the value that Inf used to:


octave:13> rand('state', NaN)
octave:14> rand
ans =  0.84442


I'm not sure what this means though, as the documentation indicates that one should do:


octave:27> v = rand ('state');
octave:28> rand ('state', v);


Perhaps the word 'state' and 'seed' are synonymous when v above is a scalar.  So, let's check 'seed':


octave:33> rand ('seed', Inf)
octave:34> rand
ans =  0.86444
octave:35> rand ('state', Inf)
octave:36> rand
ans =  0.63536


OK, so 'seed' is sort of like what 'state' used to do, but:


octave:39> rand ('seed', NaN)
octave:40> rand
ans =  0.92987
octave:41> rand ('state', NaN)
octave:42> rand
ans =  0.84442


So, 'seed' and 'state' are different entities.

Is there some consistency between versions for you with regard to "rand ('seed', s)"?

Just to get pointed in the right direction in the code, it seems that the 'state' option works it's way to (oct-rand.cc):


  void rand::set_internal_state (const uint32NDArray& s)
  {
    octave_idx_type len = s.numel ();

    const uint32_t *sdata = reinterpret_cast <const uint32_t *> (s.data ());

    if (len == MT_N + 1 && sdata[MT_N] <= MT_N && sdata[MT_N] > 0)
      set_mersenne_twister_state (sdata);
    else
      init_mersenne_twister (sdata, len);
  }


I've verified that


octave:2> rand ('state', NaN)
C++ std::cerr: Case 2
octave:3> rand ('state', Inf)
C++ std::cerr: Case 2
octave:4> rand ('state', [1 2])
C++ std::cerr: Case 2


all call the second function init_mersenne_twister(), whereas the full length vector V case calls the first function set_mersenne_twister_state();


  /* initializes state[MT_N] with a seed */
  void init_mersenne_twister (const uint32_t s)
  {
    int j;
    state[0] = s & 0xffffffffUL;
std::cerr << "init_mersenne_twister uint32 value: " << s << "\n";


This is an integer value, so these floating point numbers need to be converted.  I suspect this has something to do with that conversion from floating point to integer.  Can anyone see what the issue is from the following results?


octave:33> rand ('state', 1); rand
ans =  0.13436
octave:34> rand ('state', 0.5); rand
ans =  0.13436
octave:35> rand ('state', .499999); rand
ans =  0.84442
octave:36> rand ('state', -Inf); rand
ans =  0.84442
octave:37> rand ('state', -1); rand
ans =  0.84442
octave:38> rand ('state', NaN); rand
ans =  0.84442
octave:39> rand ('state', Inf); rand
ans =  0.63536


Could it be that in older versions the -Inf, -1, 0, NaN, and Inf all get converted to unsigned 0, whereas now the Inf gets converted to unsigned L where "L" means the largest unsigned 32-bit value (or something similar in concept by which the mantissa resolution comes into play)?  Might this be a C++11 thing?

Dan Sebald <sebald>
Fri 20 Jul 2018 05:42:40 PM UTC, comment #2: 

Thanks for responding. My problem is that I have lots of data stored after encoding using rand (some by chance seeded using rand('state',inf) and which needs to be decoded before use. I currently can't decode it on Octave 4.4.0. Your 'workaraound' does not work for me at least.

If you can tell me how to create the same rand state as rand('state', Inf) I could cope. But this bug currently prevents me moving on to Octave 4.4.0.

Ian

Ian McCallion <ianmcc>
Fri 20 Jul 2018 05:05:51 PM UTC, comment #1: 

Confirmed, that there is a difference when using 'Inf' as a seed.  On the other hand, 'Inf' is not generally a seed value I would choose.  For finite state values like 1, -5, pi, the random number generator produces equivalent results.  I am marking the bug as a regression, but lowering the priority since there is the very easy workaround of not using Inf as a seed state.

For reference with respect to other non-finite values:

Octave 4.2.2


octave:2> format long
octave:3> rand ("state", -Inf)
octave:4> rand
ans =  0.844421851525048
octave:5> rand ("state", NaN)
octave:6> rand
ans =  0.844421851525048
octave:7> rand ("state", Inf)
octave:8> rand
ans =  0.844421851525048


Octave 4.4.0


octave:2> format long
octave:3> rand ("state", -Inf)
octave:4> rand
ans =    8.444218515250481e-01
octave:5> rand ("state", NaN)
octave:6> rand
ans =    8.444218515250481e-01
octave:7> rand ("state", Inf)
octave:8> rand
ans =    6.353574441341173e-01



Rik <rik5>
Group administrator
Fri 20 Jul 2018 03:29:54 PM UTC, original submission:  

These two statements:

   rand('state', Inf)
   rand

produce 0.63536 on Octave 4.4.0. But they produce 0.84442 on Octave 4.2.0 and 3.8.2.

Maybe rand state is not specified for Inf, but given it works it is a shame it is not compatible between versions.

CIan

Ian McCallion <ianmcc>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #44603:  randstate.diff added by rik5 (870B - text/x-patch)

 

Digest:
   bug dependencies.

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by gyom (Posted a comment)
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by sebald (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by ianmcc (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 12 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-07-22 rik5 Dependencies- Depends on bugs #42557
    2018-07-22 rik5 StatusFixed Confirmed
        Open/ClosedClosed Open
    2018-07-20 rik5 StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
    2018-07-20 rik5 Attached File- Added randstate.diff, #44603
        StatusConfirmed Patch Submitted
    2018-07-20 rik5 Severity3 - Normal 2 - Minor
        Priority5 - Normal 3 - Low
        Item GroupIncorrect Result Regression
        StatusNone Confirmed
        Operating SystemMicrosoft Windows Any

    Back to the top

    Powered by Savane 3.13-cf05.
    Corresponding source code