% find bond-graph parameters whole whole cell with multiple modules
% present:
%  cAMP; LRGbinding; B1AR; PKA; PLB; Inhib1; GsProtein;

% based on Pan cardiac AP 2018 code

clear;
% clc;
close all;

%% Set directories
current_dir = cd;
Idx_backslash = find(current_dir == filesep);
parent_dir = current_dir(1:Idx_backslash(end-1));
main_dir = current_dir(1:Idx_backslash(end));
output_dir = [current_dir filesep 'output' filesep];

%% Define constants
R = 8.314; % unit J/mol/K
T = 310;
F = 96485;

cNao_st = 140; % unit mM
cNai_st = 10; % unit mM
cKo_st = 5.4;
cKi_st = 145;

N_A = 6.022e23;

A_cap = 1.534e-4; % Unit cm^2

%% Define volumes (unit pL)
V_myo = 34.4;
V_e = 5.182;     % external volume
V_SR = V_myo*0.035; % SR volume
V_di = V_myo*0.0539; % diadic volume
V = struct();
V.V_myo = V_myo;
V.V_SR = V_SR;
V.V_di = V_di;
V.V_e = V_e;

%% Load stoichiometric matrices and kinetic rate constants
array_names = {'cAMP'; 'LRGbinding'; 'beta1AR'; 'PKA'; 'PLB'; 'Inhib1'; 'GsProtein'};
num_subsystems = length(array_names);
sys_struct = cell(num_subsystems,1);
rxnIDs = {};
Knames = {};
keepKname = {};
for i_system = 1:num_subsystems
    sys_name = array_names{i_system};
    sys_dir = [main_dir sys_name '\MATLAB\'];
    cd(sys_dir)
    sys_struct{i_system}.name = sys_name;
    forward_mat_path = ['data\all_forward_matrix.txt'];
    reverse_mat_path = ['data\all_reverse_matrix.txt'];
    N_f = load_matrix(forward_mat_path);
    N_r = load_matrix(reverse_mat_path);
    sys_struct{i_system}.N_f = N_f;
    sys_struct{i_system}.N_r = N_r;
    
%     disp(array_names{i_system});
    dims = struct();
    dims.num_rows = size(N_f,1);
    dims.num_cols = size(N_f,2);
    I = eye(dims.num_cols);
    M = [I transpose(N_f); I transpose(N_r)];
    [k_kinetic, N_cT, K_C, W] = kinetic_parameters(M, true, dims, V);
    sys_struct{i_system}.kfkr = k_kinetic;
    sys_struct{i_system}.Kc = K_C;
    sys_struct{i_system}.N_cT = N_cT;
%         sys_struct{i_system}.W = W(dims.num_cols+1:end);
    fID = fopen('data\rxnID.txt', 'r');
    rxnID = textscan(fID,'%s', 'delimiter','\n');
    fclose(fID);
    rxnIDs = [rxnIDs; rxnID{1}];

    fID = fopen('data\Kname.txt', 'r');
    Kname = textscan(fID,'%s', 'delimiter','\n');
    fclose(fID);
    Knames = [Knames; Kname{1}]; 
end

cd(current_dir)
% relations between submodule to whole module

sys_struct{1}.I_vec = 1:15;           % cAMP
sys_struct{2}.I_vec = 16:21;          % LRGbinding
sys_struct{3}.I_vec = [19 21 22:29];  % B1AR
sys_struct{4}.I_vec = [2 30:42];      % PKA
sys_struct{5}.I_vec = [43 28 44:49];  % PLB
sys_struct{6}.I_vec = [28 48 50:53];  % Inhib1
sys_struct{7}.I_vec = [17 18 20 19 21 14 54:55];  % GsProtein
num_rows = sys_struct{size(sys_struct,1)}.I_vec(end);

N_f = [];
N_r = [];

for i_system = 1:num_subsystems
%     disp(array_names{i_system});
    T = calcT(sys_struct{i_system}.I_vec,num_rows);

    sys_struct{i_system}.T = T;
    
    N_f = [N_f T*sys_struct{i_system}.N_f];
    N_r = [N_r T*sys_struct{i_system}.N_r];
end

N_fT = transpose(N_f);
N_rT = transpose(N_r);

N = N_r - N_f;
N_T = N_rT - N_fT;

num_cols = size(N,2);
I = eye(num_cols);

M = [I N_fT; I N_rT];
M_rref = rref(M);

%% Set up the vectors for kinetic rate constants
%  cAMP; LRGbinding; B1AR; PKA; PLB; Inhib1; GsProtein

kf = [];
kr = [];
W = ones(size(N,2),1);
for i_system=1:num_subsystems
    nrx = length(sys_struct{i_system}.kfkr)/2;
%     disp(array_names{i_system});
%     disp(nrx)
    kf = [kf, sys_struct{i_system}.kfkr(1:nrx)'];
    kr = [kr, sys_struct{i_system}.kfkr(nrx+1:end)'];
%     W = [W; sys_struct{i_system}.W]; 
end
k_kinetic = [kf, kr]';
W = [ones(size(N,2),1);V_myo*ones(15,1);V_myo*ones(14,1);V_myo*ones(13,1);
    V_myo*ones(7,1);V_myo*ones(4,1);V_myo*ones(2,1)];

lambdaW = exp(pinv(M)*log(k_kinetic));
lambda = lambdaW./W;
kappa = lambda(1:size(N,2));
K = lambda(size(N,2)+1:end);

save([output_dir 'whole_saucerman03_BG_parameters.mat'],'kappa','K','k_kinetic');

%% Checks
N_rref = rref(N);
R_mat = null(N,'r');

k_est = exp(M*log(lambdaW));
diff = (k_kinetic - k_est)./k_kinetic;
error = sum(abs(diff));

K_eq = kf./kr;
zero_est = R_mat'*K_eq';


%% print outputs
for ik=1:length(kappa)
    fprintf('var kappa_%s: fmol_per_sec {init: %g, pub: out};\n',rxnIDs{ik},kappa(ik));
end
Kunique = {};
for ik=1:length(Knames)
    if ~any(strcmp(Kunique,Knames{ik}))
        Kunique = [Kunique; Knames{ik}];
    end
end
for ik = 1:length(Kunique)
    fprintf('var K_%s: per_fmol {init: %g, pub: out};\n',Kunique{ik},K(ik));
end
%% FUNCTION SCRIPTS
function matrix = load_matrix(stoich_path)
    stoich_file_id = fopen(stoich_path,'r');
    stoich_file_data = textscan(stoich_file_id,'%s','delimiter','\n');
    fclose(stoich_file_id);
    num_rows = length(stoich_file_data{1});
    num_cols = sum(stoich_file_data{1}{1} == ',')+1;
    matrix = zeros(num_rows,num_cols);
    for i_row = 1:num_rows
        line_str = stoich_file_data{1}{i_row};
        line_split = regexp(line_str,',','split');

        for i_col = 1:num_cols
            matrix(i_row,i_col) = str2double(line_split{i_col});
        end
    end
end


function T = calcT(I_vec,num_rows)
    num_cols = length(I_vec);
    T = zeros(num_rows,num_cols);
    for i = 1:num_cols
        T(I_vec(i),i) = 1;
    end
end
