How to retrieve simulation data

[Haotian Wanght]

Hi Tjalling,
How to retrieve the simulation data points? It appears that there is only figure provided.

[Tjalling]

The trick is to call calc_and_plot with an output:

Xall = calc_and_plot(...);

The variable Xall will then contain the same values as shown in the plot. It will be a cell array, so you'll have to figure out what is what. I think that Xall{1} will contain a matrix for the first state variable, with time in the first column, and all treatments in the next columns. But please check.

[Haotian Wanght]

Hi Tjalling,

Thanks for your reply. I tried with Xall = calc_and_plot(...);. but matlab showed an error.

Then I tried 'xall = sim_and_plot(par,opt_sim);', but xall is a figure instead of a matrix.

[Tjalling]

What is the error you got? If your call to calc_and_plot runs without error (and produces a plot), I don’t see how adding an output variable could produce an error ...

[Haotian Wanght]

Hi Tjialling,
Here are the codes.

#1 I tried Xall = calc_and_plot(...);

Code: Select all

%% Initial values for the state variables
% Initial states, scenarios in columns, states in rows. First row are the
% 'names' of all scenarios.

% these scenarios specify different starting conditions: N(0)
X0mat(1,:) = [0.00004 0.000045 0.00005]; % scenario identifiers
X0mat(2,:) = [0 0 0]; % initial value for state 1
X0mat(3,:) = [0 0 0]; % initial value for state 1
X0mat(4,:) = [0 0 0]; % initial value for state 1

   
%% Initial values for the model parameters
% Model parameters are part of a 'structure' for easy reference. 

% syntax: par.name = value
par.ku   = 4339;    
par.ke   = 0.2693;   
par.ksr  = 1.223;    
par.ksl  = 0.09231;  
par.krs  = 0.1181;   

%% What to calculate and what to plot
% Specify what to calculate and what to plot.

glo.t = linspace(0,10,50); % time vector for the model curves in days

% specify the y-axis labels for each state variable
glo.ylab{1} = 'leaf concentration';
glo.ylab{2} = 'stem concentration';
glo.ylab{3} = 'root concentration';

% specify the x-axis label (same for all states)
glo.xlab    = 'time (t)';
glo.leglab1 = 'C_w = '; % legend label before the 'scenario' number
glo.leglab2 = ' ppm'; % legend label after the 'scenario' number

% set options for call_deri here (we don't do prelim_checks for simulations)
glo.useode   = 1; % calculate model using ODE solver (1)
glo.eventson = 0; % events function on (1) or off (0)
glo.stiff    = 0; % ODE solver 0) ode45 (standard), 1) ode113 (moderately stiff), 2) ode15s (stiff)

%% Calculations and plotting
% Also try plottype 5 for this example! See the SIMbyom package for more
% detailed information on the simulator options.
 
opt_sim.plottype = 3; % select a plotting option
% 1) 3d plot, 
% 2) 2d plot, 
% 3) states in subplots, 
% 4) scenarios in subplots,
% 5) dx/dt versus x
% 6) 2d plot for each scenario separate

sim_and_plot(par,opt_sim); % call the script which calculates and plots (simulation only)
Xall = calc_and_plot(...);

Than I got an error
>> byom_plant_sim
Error: File: byom_plant_sim.m Line: 109 Column: 1
Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for
mismatched delimiters.

Afterward, I tried Xall = calc_and_plot(par,opt_sim);
The figure can be plotted, but there was still an error

Reference to non-existent field 'zvd'.

Error in calc_and_plot (line 27)
zvd_plot = opt_plot.zvd; % plotting of zero-variate data (if defined)

Error in byom_plant_sim (line 108)
Xall = calc_and_plot(par,opt_sim);

[Tjalling]

instead of the ..., put the code that was there before . Or better, just put Xall= before the call to calc_and_plot.

[Haotian Wanght]

Thanks for the reply.

But it appears that there is no 'calc_and_plot(par_out,opt_plot)' in the simulation. In the simulation file, it shows the 'sim_and_plot(par,opt_sim)'

I tried several times with , but there are still errors.

Code: Select all

%% Calculations and plotting
% Also try plottype 5 for this example! See the SIMbyom package for more
% detailed information on the simulator options.
 
opt_sim.plottype = 3; % select a plotting option
% 1) 3d plot, 
% 2) 2d plot, 
% 3) states in subplots, 
% 4) scenarios in subplots,
% 5) dx/dt versus x
% 6) 2d plot for each scenario separate

Xall = calc_and_plot(par,opt_sim); 
sim_and_plot(par,opt_sim); % call the script which calculates and plots (simulation only)
%Xall = sim_and_plot(par,opt_sim);

Code: Select all

>> byom_plant_sim
Warning: Multiple folders with BYOM in the name have been found above this location; lowest is used! 
Reference to non-existent field 'zvd'.

Error in calc_and_plot (line 27)
zvd_plot = opt_plot.zvd;     % plotting of zero-variate data (if defined)

Error in byom_plant_sim (line 107)
Xall = calc_and_plot(par,opt_sim);

[Tjalling]

Ah, you are using sim_and_plot ... then this does not work. There is no easy way to get the data out, as plotting is done ‘on the fly’ so that it’s kind of animated.

You could use a regular script file (the ones that fit), and simply don’t fit anything. There’s an option for that ... I will have to check later, and add it to this post. That way, you would be able to use calc_and_plot, with extra output.

Edit: take a look at the example file in BYOM names byom_bioconc_start.m. There is the following option:

opt_optim.fit = 1; % fit the parameters (1), or don't (0)

If you set it to zero, there will be no fitting, just simulation (but only states-vs-time, and no animation).

[Haotian Wanght]

Yes and Thanks so much for the help.

I tried Xall = calc_and_plot(par_out,opt_out) in the fit file, and it did work out with a matrix .

The reason I asked for acquiring simulation data is because I am wondering if the BYOM could do some reverse dosimetry work?
For example, I could measure the concentration of a chemical in one-compartment (i.e. blood), then I could predict the external concentration (i.e. the chemical concentration in the water phase).

Now, I just gave some initial value of water concentration, and then compared the simulation values with measured values, by which judge if the given value is approximate to the expected one.

[Tjalling]

The reason I asked for acquiring simulation data is because I am wondering if the BYOM could do some reverse dosimetry work?
For example, I could measure the concentration of a chemical in one-compartment (i.e. blood), then I could predict the external concentration (i.e. the chemical concentration in the water phase).

Something like that would be possible. If your external concentration is constant, it would be quite simple: make the external concentration a model parameter, so you can fit it (with a CI if needed). If the external concentration is time-varying, things are more awkward, though you could try to come up with a function for the external concentration, and make the parameters of that function model parameters to be fitted. I hope this gives you some ideas of what to try.