Restart random numbers with parallel loop

I want to be able to restart my simulation from where I left off and also for it to be reproducible

if isfile("State.mat")
    % If we previously saved the state, start from there
    load("State");
    globalStream.State = rngState;
    ...
else
    % Otherwise start from scratch
    rng(42,'twister');
    globalStream = RandStream.getGlobalStream;
    rngState = globalStream.State;
    ...

    save("State", "current_time", "B", "rngState");
    dlmwrite("data/tvol_sim_mat.csv", tvol_mat(1, :), 'precision', 15);
end

...

for t = 1:2

    for i = 1:n

        % Simulates the system forward

    end

    globalStream = RandStream.getGlobalStream;
    rngState = globalStream.State;

    save("State", "current_time", "B", "rngState");
end

But I want to replace for by parfor. I've read the docs and tried

if isfile("State.mat")
    % If we previously saved the state, start from there
    load("State");
    s = RandStream('mlfg6331_64', 'Seed', parameters_for_inference.filter.seed);
    options = statset('UseParallel', true,  'Streams', s, 'UseSubstreams', true);
    s.State = rngState;
    ...
else
    % Otherwise start from scratch
    s = RandStream('mlfg6331_64', 'Seed', parameters_for_inference.filter.seed);
    options = statset('UseParallel', true,  'Streams', s, 'UseSubstreams', true);
    rngState = s.State;

    save("State", "current_time", "B", "rngState");
end

...

for t = 1:2

    parfor i = 1:n

        % Simulates the system forward

    end

    rngState = s.State;

    save("State", "current_time", "B", "rngState");
end

But I get this error when trying to restart

Error using matlab.internal.math.RandStream_getset_mex
State array class is invalid for a mlfg6331_64 generator.

Error in RandStream/subsasgn (line 646)
                    matlab.internal.math.RandStream_getset_mex('state',a.StreamID,b);

Error in emacsrunregion (line 23)
    evalin('base',evalTxt);

Is it possible to do what I want in Matlab? Apologies if this is obvious but I am very new to Matlab. Many thanks :)

Answer 1

[Edric Ellis]

The error that you're seeing there I think is because you've got an old version of your "State.mat" file - the error is the one that you receive when you try to set the State property of a 'mlfg6331_64' generator using the State extracted from a 'twister' generator. So maybe you simply need to delete or rename your old "State.mat" file.

Also, can I suggest you review this page which describes how to use Parallel Computing together with Statistics and Machine Learning Toolbox. Basically, if you're using statset with option 'UseParallel' -> true, then you do not need to write the parfor loop yourself to get parallel execution of the stats functions. (Although perhaps you're not calling stats functions inside your parfor loop?)

[Dominic Steinitz]

Thanks very much - I think that must have been the problem - just knowing what I wanted to do was indeed possible helped.

Thanks for the pointer and I have already read that page. I am using someone else's model so I can't be fully sure but afaics, the model just uses random numbers and no stats functions. I am running a particle filter (not the Matlab one as it makes assumptions about the (Markov) probability kernel) on the model. I'd like to run it on a GPU but that might be a bit ambitious given how new I am to Matlab.

Thanks once again :)

[Edric Ellis]

Ah I see. In which case this page will tell you more about reproducing random numbers inside parfor. That statset call is only relevant for stats functions.

[Dominic Steinitz]

Thanks very much for this but I have e.g.

r_vec_pd = makedist('Normal', 'mu', 0.0, 'sigma', r_sigma);
...
random(r_vec_pd, 1, 1))

And there doesn't seem to be a way to say which stream to use for the call to random?

Also I think this means having to go through someone else's model and change the code. In Haskell, I would create as many generators as I need and then run each thread with its own generator with no changes to the model. I hope I have misunderstood :)

[Edric Ellis]

For methods that don't accept the RandStream as an input, you can set the "global" stream using RandStream.setGlobalStream. You'd want to wrap your code something like this (untested...)

oldStream = RandStream.getGlobalStream();
restoreStream = onCleanup(@() RandStream.setGlobalStream(oldStream));
RandStream.setGlobalStream(someStream);
% Now, "random" will use "someStream"; when "restoreStream" goes out of
% scope, MATLAB's global stream will be restored.