%main function to generate results presented in the paper

clear
clc
warning off

% add the path that contains the structure of network models
addpath('Benchmark_Networks');
% add the path that contains all CCELL codes
addpath('../CCELL');
% add the path that contain Bayesian compressive sensing package
addpath('../CCELL/SB2_Release_200');

%parameter settings
T1 = 5;% number of time points in each perturbation
ss = 300;% number of species in the network
os = ss/2;% number of combined measurements
genNewC = 0;% if ==1,then Low-density parity-check code will generate new C
genNewXY = 0;% if ==1,then generate new x,y
noiseR = 0.0001;% assumed variance of measurement noise
noiseQ = 0.0001;% assumed variance of process noise
N_per_all = 2:1:7;% different number of perturbations
% 7 different network models to be investigated in the paper
A_name = {'A_JAK_STAT4.mat','A_ERK11.mat','A_p38_39.mat','A_p38_ERK_50.mat','A_p38_ERK_GR_53.mat','A_p38_ERK_GR_JAK_58.mat'};
Num_non_zeros = [4,11,39,50,53,58];% number of non-zero entries in transition matrices of those 7 models
N_sparsity = length(A_name);% number of network models: 7

%% network inference
% for different number of perturbations
for N_per = N_per_all
    N_per
    T = N_per*T1;% total number of time points
    Perf = cell(N_sparsity,1);% stores the results
    % for different pathway
    for i = 1:N_sparsity
        % use the same seed to generate random number
        seed = 2;
        rand('state', seed);
        randn('state', seed);
        %load the transition matrix A
        A_i = load(A_name{i});
        Sname = fieldnames(A_i);
        A_path = A_i.(Sname{1});
        %get the number of non-zero entries of current A
        num0 = Num_non_zeros(i);
        %generate the whole transition matrix with some rows having only one
        %non-zero entry corresponding to time-invariant species
        A_whole = eye(ss);
        A_whole(1:num0,1:num0) = A_path;
        A = A_whole;
        %call pathway inference program to get results
        if genNewXY==1% generate new x and y
            x_all=[];
            y_all =[];
        elseif genNewXY==0 % read in x and y
            load(strcat('data/x_all_sparsity',num2str(i),'_perturbation',num2str(N_per),'.mat'))
            load(strcat('data/y_all_sparsity',num2str(i),'_perturbation',num2str(N_per),'.mat'))
        end
        [Result_i  x_all y_all]= pathway_inference_benchmark(N_per,A,ss,os,T,T1,genNewC,noiseR,noiseQ,x_all,y_all);
        Perf{i} = Result_i;
        if genNewXY==1% save newly generated x and y
            save(strcat('data/x_all_sparsity',num2str(i),'_perturbation',num2str(N_per),'.mat'),'x_all')
            save(strcat('data/y_all_sparsity',num2str(i),'_perturbation',num2str(N_per),'.mat'),'y_all')
        end
    end
    %save results
    save(strcat('result/Purterbation',num2str(N_per),'Perf.mat'),'Perf')
end

% plot results
plot_result