Issues with doses/variants in parfor loop in SimBiology

I have a question about runing the sbproj models in parfor loop

1. I load project and set up the solver:

sbioloadproject('modelName.sbproj','m1')
csObj = getconfigset(m1);
csObj.SolverType = 'sundials';
csObj.SolverOptions.AbsoluteToleranceScaling = 1;
csObj.SolverOptions.AbsoluteTolerance = 1E-16;
csObj.SolverOptions.RelativeTolerance = 1E-12;
csObj.MaximumWallClock = 60;

2. I define dose and variant:

dose1 = getdose(m1,'myDose1');
v1 = getvariant(m1,'myVariant1');

3. accelerate the model

sbioaccelerate(m1,csObj,v1,dose1);

4. simulate in parfor loop

m1ParConst = parallel.pool.Constant(m1);
parfor i = 1:r
    [...]
    % dose1 = getdose(m1ParConst.Value,'rela_3mgpkg_80kg_q2wx31');
    % v1 = getvariant(m1,'myVariant1');
    [t,y] = sbiosimulate(m1ParConst.Value,[],v1,dose1);
    [...]
end

Question:

Q1. Is this correct? If I include dose1 and v1 definition (i.e. if I uncomment dose1 and v1 in parfor loop) I get an error:

Undefined function 'getdose' for input arguments of type 'double'.
Error in parallel_function>make_general_channel/channel_general (line 923)
                    O = F(C{:});
Error in remoteParallelFunction (line 46)
        out = parallel.internal.pool.serialize(feval(channel, channelArgs{:}));

Q2: are the solver settings respected? If I try to include csOBj in 'sbiosimulate' call I get an error as well.

[Sietse Braakman]

Which version of MATLAB are you using?

[emjey]

R2018b, I have to check if I can switch to later version but for now I have this setup and it would be great if I can get it to work.

Update: I am installing R2019b now and keep you posted what happens.

Answer 1

[Arthur Goldsipe]

First, what version of SimBiology are you using? Prior to R2019b, there was a problem that prevented objects from getting sent to workers properly. That is usually the cause of errors like "Undefined function 'getdose' for input arguments of type 'double'."

Second, it can be difficult to get details right with acceleration and parallelization when using sbioaccelerate and sbiosimulate. For example, sbioaccelerate is a "local" operation that is not preserved when transferring a model across workers. So if you want an accelerated model on a worker, you have to run sbioaccelerate on that worker.

As an alternative, I recommend using a SimFunction instead. It takes care of all of these details and is generally faster and more efficient. Let me show you some sample code. For now, I will assume you want to simulate varying the values of parameters p1 and p2, and look at the simulation results for y1 and y2. (You can also vary the dosing, or specify your alternate paramaeters using variants.)

1. Load project and set up the solver:

sbioloadproject('modelName.sbproj','m1')
csObj = getconfigset(m1);
csObj.SolverType = 'sundials';
csObj.SolverOptions.AbsoluteToleranceScaling = 1;
csObj.SolverOptions.AbsoluteTolerance = 1E-16;
csObj.SolverOptions.RelativeTolerance = 1E-12;
csObj.MaximumWallClock = 60;
tStop = csObj.StopTime;

2. Define dose and parameter values:

dose1 = getdose(m1,'myDose1');
doseTable = getTable(dose1);
p1 = [1;2;3];
p2 = [4;5;6];

3. Create the SimFunction (ignoring any warning about the dose):

func = createSimFunction(m1, ["p1", "p2"], ["y1", "y2"], dose1, "UseParallel", true);

4. Simulate in parallel:

[timeCell, yCell] = func([p1 p2], tStop, doseTable);

[emjey]

Thanks Arthur, but I am not sure how I can translate the code below (Morris global SA based on calculation of elementary effects). I can live without acceleration but I have to be able to pass csObj, variants and doses. Any suggestions would be very apreciated.

sbioloadproject(['SBPROJ/',modelName],'m1')
csObj = getconfigset(m1);
% solver setup
csObj.SolverType = 'sundials';
csObj.SolverOptions.AbsoluteToleranceScaling = 1;
csObj.SolverOptions.AbsoluteTolerance = 1E-16;
csObj.SolverOptions.RelativeTolerance = 1E-12;
csObj.MaximumWallClock = 60;

% dose and variant
dose1 = getdose(m1,'myDose');
v1 = getvariant(m1,'myVariant');
            
% elementary effects array
EE = zeros(r,p);
m1ParConst = parallel.pool.Constant(m1);
% run r trajectories
parfor i = 1:r
   
    f = zeros(1,(p+1));
    % calculate a parameter array for i-th trajecory 
    C = morrisCmatrix(p,l,paramsDF,delta);

    for j = 1:(p+1)
        try
            % assign model parameters 
            % columns in C matrix = parameter sets
            for k = 1:size(C,2) 
                blub = sbioselect(m1ParConst.Value, 'Name', deblank(paramsDF.parameter(k)));
                blub.Value = C(j,k);
            end
            [t,y] = sbiosimulate(m1ParConst.Value,csObj,v1,dose1);
                
            f(j) = y(end,varNo);
        catch me
            disp( getReport( me, 'extended', 'hyperlinks', 'on' ) )
            f(j) = NaN;
        end
    end
    
    v = zeros(1,p);
    for j = 1:p
        indexOfCm = find(diff(C(:,j)) ~= 0);
        v(j) = (f(indexOfCm) - f(indexOfCm+1))/delta
    end
    EE(i,:) = v;
end

Update: I am installing R2019b now and keep you posted what happens and if I can run my original code.

[Arthur Goldsipe]

I don't see any reason why you can't use SimFunction. If you were modifying configuration settings between runs, then it might not make sense to use a SimFunction. But the only things you're changing between runs are parameter values, so I would definitely recommend that approach. I can show you how to rewrite this code to use a SimFunction instead, and I think you'll find it simpler and definitely much faster to run. But before I take the time to do that, I wanted to see whether it was worth the effort. If you just want to use your original code, then I expect that upgrading to R2019b will eliminate the errors you were seeing.

[emjey]

First I have to find out the reason for the error:

Warning: Reported from ODE Compilation:
Following error occurred while testing a MEX file created for SimBiology acceleration: Accelerated results for initialAssignment
rules did not match the unaccelerated results. Check that initialAssignment rule expressions and any user-defined functions are
valid and supported for code generation as described in the Language, Function, and Object Support in MATLAB documentation. If
this problem persists, contact Technical Support. Click here for more information.
(Type "warning off SimBiology:CodeGeneration:MexFileEvaluationFailed" to suppress this warning.)
> In sbioaccelerate (line 71)

I didn't get it with the same model using R2018b.

[Arthur Goldsipe]

You should be able to ignore that warning. It indicates a potential loss of performance. But since you're only seeing the warning for initial assignment rules, the performance impact should be negligible to non-existant.

The reason you didn't see this message before is because this is a new check that I introduced in R2019b to catch potential problems where the behavior of the accelerated code differs from the unaccelerated code. After customers started using R2019b, I discovered a few issues with that check. One of the issues was fixed in R2019b update 4, but there are still a few more issues I need to fix.

If you can answer a few questions, it would help me figure out what else I need to fix for this issue. Are you using R2019b update 4? Do some of your initial assignment rules involve raising a number to a fractional power?

[emjey]

'9.7.0.1296695 (R2019b) Update 4', and yes, we have fractional power in initial assignments in the model

[Arthur Goldsipe]

Ok, it sounds like you're running into the issue I know about. If you want any more help (for example, in updating your sample code to use a SimFunction), just let me know.

[emjey]

Apologies Arthur for the delay. For start I would like to be able to use efficienlty the R2019b version but I ran into problems.

In the meantime I've experimented with 'parfor' simulation of the 'AntibacterialPKPD.sbproj' model which I extended a bit, added one long dosing, 'myDose', and increased the simulation time accordingly. I did that to test/learn acceleration and parallel setup.

As you can see from the figures, there is no speedup in the accelerated version. For 25 loops, each with 11 simulations, takes

• 5.72 min with an apparently accelerated model
• 6.0 min with a non-accelerated model

It must be something I get wrong in the script, AntibacterialPKPD_PARFOR.m. Attached 1) sbproj file, 2) parameter file, 3) M-file.

Any comments would be very appreciated.

[Sietse Braakman]

This problem does not exist in MATLAB/SimBiology version R2019b and later.

Prior to 19b, SimBiology objects are not automatically transferred to the workers when a parfor loop is used. You can circumvent this storing the necessary objects in the parallel.pool.Constant. Once you have 'unpacked the parallel.pool.Constant values for each worker, you can accelerate the model inside the parfor loop. Alternatively, you can create a simfunction outside the parfor loop, attach that to the parallel.pool.Constant and then use the simfunction inside the parfor loop. I have attached an updated version of your script. Note that I expected the simfunction to give me the same values of Drug.Central at t = 3360 as the sbiosimulate approach but this doesn't seem to be the case. You will need to check what's going on there, but at least you'll be able to see how to use the simfunction in this case. Perhaps I switched a parameter order or dimension of C around?

I am posting the relevant additions/changes below. Let me know if that works.

loopTimeTot = zeros(r,p+1); 
params = sbioselect(m1,'Where','Name','==',paramsDF.parameter);
simfun = createSimFunction(m1,paramsDF.parameter,varName,d1.TargetName);
c = parallel.pool.Constant({m1,simfun});
accTimes = nan(r,1);
parfor i = 1:r
    %%%%%%%%%%% SETUP %%%%%%%%%%%
    model = c.Value{1}; % grab model from parallel.pool.Constant
    dose = getdose(model,'myDose'); % extract dose from model
    configset = getconfigset(model); % extract configset from model
    % start acceleration, note that the acceleration is only performed once
    % on each worker.
    startAcc=tic;
    sbioaccelerate(model,configset,dose);
    accTimes(i)=toc(startAcc);
    % inputs for simfunction
    simfunction = c.Value{2}; % grab simfunction from parallel.pool.Constant
    doseTable = getTable(dose); % create a table from the dose object
    stopTime = configset.StopTime; % grab the stoptime from the configset
    %%%%%%%%%%%%%%%%%%

...

    [t,y] = simfunction(C,stopTime,doseTable);
end

[emjey]

Thank you, it works indeed for the example model. Acceleration brings the simulation time down to 20% of the default one!