AMICI icon indicating copy to clipboard operation
AMICI copied to clipboard

Steady-state sensitivities in Matlab: sllh issue

Open aidinbii opened this issue 2 years ago • 5 comments

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);

aidinbii avatar Apr 13 '23 06:04 aidinbii

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 avatar Apr 15 '23 06:04 FFroehlich

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

FFroehlich avatar May 09 '23 08:05 FFroehlich

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

aidinbii avatar May 10 '23 06:05 aidinbii

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?

FFroehlich avatar May 10 '23 07:05 FFroehlich

@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

aidinbii avatar May 17 '23 06:05 aidinbii