Using 'time' in observable expressions
显示 更早的评论
Hello,
Is it possible to select only certain timepoints for observable calculations, for example:
mdl.addobservable('maxA','max(A(time>(3*7*24)))','Units','dimensionless');
In my example, A is a parameter that is set by a repeated assignment rule.
When I run the model in R2021b, I get the following error, whih suggests incorrect dimensions (but my observable expression should return a scalar, no?):
------------------------------------------------------
When running sbiosteadystate:
Error using sbiosteadystate (line 128)
sbiosteadystate encountered a model compilation error.
Caused by:
Error using SimBiology.internal.verifyHelper
--> Error reported from Expression Validation:
The result of evaluating observable 'maxA' has the wrong size. It must either be a scalar or a column vector of the same length as the time vector.
-------------------------------------------------------
When running sbiosimulate:
Error using SimBiology.internal.simulate
--> Error reported from Expression Validation:
The result of evaluating observable 'maxA' has the wrong size. It must either be a scalar or a column vector of the same length as the time vector.
----------------------------------------------------
Please advise.
Thank you,
Abed
采纳的回答
Florian Augustin
2022-2-3
0 个投票
Hi Abed,
when evaluating expressions of observables, SimBiology first tries to determine whether the result is a scalar value or a vector. I therefore recommend to write observables' expressions in a way that they either always return a scalar value or always a vector of values - for any possible input values of referenced model components and time.
In your case, the expression max(A(time > 3*7*24))) can be a vector whose length differs from the length of the time vector, e.g. if none of the values in time exceeds 3*7*24, or if the maximum is attained at multiple time points. You can write a helper function that returns the max. value and protects against invalid edge cases:
function maxValue = getMaxValueAfterSpecifiedTime(values, time, threshold)
% Get max value (call unique to ensure scalar or empty output).
maxValue = unique(max(values(time > threshold)));
% Protect against empty maxValue by returning nan.
if isempty(maxValue)
maxValue = nan;
end
end
You can then call this helper as follows:
% Load/build your model
mdl = sbmlimport("radiodecay");
% Observable referencing time point
observable = addobservable(mdl, "maxX", "getMaxValueAfterSpecifiedTime(x, time, 2)");
set(observable, "Active", true);
% Set simulation options
configset = getconfigset(mdl);
configset.RuntimeOptions.StatesToLog = "x";
% Simulate and plot model state and active observable
simData = sbiosimulate(mdl);
sbioplot(simData);
When you run this particular example, you may see that maxX differs from the actual max value as shown in the plots. You can add interpolation to your helper function to increase the accuracy of the reported max value; an example how this could look like is attached to this post.
I hope this helps.
Best,
-Florian
4 个评论
Abed Alnaif
2022-2-3
Florian Augustin
2022-2-3
Hi Abed,
I'm glad it helped!
The line where I set the Active property on the observable is actually not necessary, the default value for this property is already true. I only included it to make it obvious at first sight what to expect to see in the plot when running the code (active observables and states (values of species, parameters, compartments) listed in the configuration set's StatesToLog property).
Best,
Florian
Arthur Goldsipe
2022-4-26
I just had an idea for how to avoid writing a separate function getMaxValueAfterSpecifiedTime. In many cases, I think you could instead change
max(expression)
to
max(vertcat(nan,expression))
to make sure that max is never called on an empty vector.
Abed Alnaif
2022-4-27
更多回答（0 个）
社区
更多回答在 SimBiology Community
类别
在 帮助中心 和 File Exchange 中查找有关 Scan Parameter Ranges 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
