AMICI-dev / AMICI

Advanced Multilanguage Interface to CVODES and IDAS
https://amici.readthedocs.io/
Other
108 stars 30 forks source link

Steady-state sensitivities in Matlab: sllh issue #2037

Open aidinbii opened 1 year ago

aidinbii commented 1 year ago

Hi, could you please help me out In Matlab when I run steady-state sensitivities in the case of adjoint sensitivities

options_ss.sensi_meth = 'adjoint';

sllh vector has all zeroes.

while in forward sensitivities

options_ss.sensi_meth = 'forward';

it calculates non-zero values

I provide random data

D.Y = [rand(1,30)*50];
D.Sigma_Y = rand(1, 30);
FFroehlich commented 1 year ago

Thanks for bringing this up, this is probably some issue that options are not passed correctly in the matlab interface.@plakrisenko could you please have a look at https://github.com/AMICI-dev/AMICI/blob/master/src/interface_matlab.cpp and see whether we missed adding some of the options to the matlab interface that were added to the python interface?

FFroehlich commented 1 year ago

Hi aidinbii, are you sure that the integration was actually successful, i.e. that no integration error occured?

aidinbii commented 1 year ago

@FFroehlich, yes, the posteq_status is -4; 1; 0 (for both cases: forward & adjoint)

FFroehlich commented 1 year ago

okay, could you please try with the latest version of amici (0.17.0) that was just released and, if the error persists, upload an example script that reproduces the error?

aidinbii commented 1 year ago

@FFroehlich , it did not help

Here is the script:

clc
clear
close all

%% get tuned parameters
load parameter_Set_healthy_ini_cond_SS_SCRNA_new

p = exp(parameters.MS.par(:,1)); % exp bc they are in log scale
k = initial_condition;

%% Simulations

% indices of inactive cells
inact_ind = [1;7;16;20];

% indices of active cells
act_ind = [2;3;8;9;11;14;17;21;22;23];

% indices of proteins
prot_ind = [4;5;6;10;12;13;15;18;19;24];

%%
t = linspace(0,1000,3000);
options_simul = amioption('sensi',1,'sensi_meth','adjoint', ...
    'maxsteps',1e6);
modelName = 'simulate_IgG4_ver1_SCRNAseq_MM1';

%% Random data
D.Y = [rand(1,10)*50];
D.Sigma_Y = rand(1, 10);

%% test adjoint sensi.
out = searchSteadyState(modelName, p, k, options_simul, D);

where the function searchSteadyState:

function sol = searchSteadyState(modelName, param, k_ss, options_ss, varargin)
t = Inf;

[~] = which(modelName); % fix for inaccessability problems

model = str2func(modelName);

if nargin > 4
    D_ss = varargin{1};
    sol = model(t,param,k_ss,D_ss,options_ss);
else
    sol = model(t,param,k_ss,[],options_ss);
end

%% Diagnosis
disp(['The status is ', num2str(sol.diagnosis.posteq_status')]);
disp(['Timepoint reached ', num2str(sol.diagnosis.posteq_t')]);
end