close all
clc

% the path of the toolbox
addpath('/MATLAB/BNPBN/');
addpath('/MATLAB/BNPBN/BNPBN');
% add path to the directory that contains TruthTable.m, 
% TruthTable_del_a_temp.m, and TruthTable_del_n_temp.m

% Boolean model of p53 from paper:
% Attractor Landscape Analysis Reveals Feedback Loops in the p53 Network 
% That Control the Cellular Response to DNA Damage, Choi et al.

Names={'ATM 1' , 'p53 2' , 'Mdm2 3' , 'MdmX 4' , ... 
       'Wip1 5' , 'cyclinG 6', 'PTEN 7 ', 'p21 8', ...
       'AKT 9', 'cyclinE 10', 'Rb 11', 'E2F1 12', ...
       'p14ARf 13' , 'Bcl2 14', 'Bax 15' , 'caspase 16' };
  
%Number of inputs
n = 16;
p = 2;
% nf = ones(1,16);
c = ones(2,16);

nv =  [4 4 10 6 2 2 2 4 2 3 5 5 4 4 3 6 ];
% normal cells
varF = [ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
         5  1  9  1  2  2  2  9  7 12  1 13 12  9 14  9 
         6  3 13  9 -1 -1 -1  3 -1  8 10  1  2 16  2 15
        12  4  1 13 -1 -1 -1  2 -1 -1 16  3  5  2 -1 14 
        -1 -1 10  3 -1 -1 -1 -1 -1 -1  3 11 -1 -1 -1 12 
        -1 -1  6  5 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  8
        -1 -1  4 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
        -1 -1  2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
        -1 -1 11 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
        -1 -1  5 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1];
        
for i = 1 : 16
   I = find(varF(:,i)==-1);
   varF(I,i) = 100;
end
    
% Base line Normal system DNA 0
BaseDNA0 = [-1 1 1 1 -1 -1 -1 -1 1 1 0 1 -1 0 -1 -1];

% Base value with damage DN 1
BaseDNA1 = [2 1 1 1 -1 -1 -1 -1 1 1 0 1 0 0 -1 -1];
     
% Inputs Interactions weight
Weight=[ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
        -2  2  1 -2  2  2  2 -1 -2  1  2 -3  1  2 -1 -1
        -2 -3 -3  2  0  0  0 -1  0 -3 -1  1 -2 -1  2  3
         1 -2 -2 -3  0  0  0  3  0  0 -3  1 -1 -2  0 -1 
         0  0 -2 -4  0  0  0  0  0  0 -1 -1  0  0  0  1 
         0  0  1  1  0  0  0  0  0  0  0  0  0  0  0 -1
         0  0  2  0  0  0  0  0  0  0  0  0  0  0  0  0
         0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0
         0  0 -2  0  0  0  0  0  0  0  0  0  0  0  0  0
         0  0  2  0  0  0  0  0  0  0  0  0  0  0  0  0];
     
W2 =zeros(10,16);     
for i = 1 : 16
   [Y,I] = sort(varF(:,i));
   varF(:,i) = Y;
   W2(:,i) = Weight(I,i);
end

for i = 1 : 16
   I = find(varF(:,i)==100);
   varF(I,i) = -1;
end
%%
% Compute table for
% F = TruthTable(n,nv,varF,BaseDNA0, W2);
F = TruthTable(n, nv,varF, BaseDNA1, W2);

% delete arrows and nodes temporarily
Cnodes = (6);% constantly express 'cyclinG 6'
F1 = TruthTable_del_n_temp(F,nv,varF,p,Cnodes,1);
Cnodes1 = (7);% 'PTEN 7 '
F2 = TruthTable_del_n_temp(F1,nv,varF,p,Cnodes1,0);
Cnodes2 = (13);% 'p14ARf 13'
F3 = TruthTable_del_n_temp(F2,nv,varF,p,Cnodes2,0);
% for fixed pt and blocking
cedges = [2 3 3  8  1  3 4 11 14
          5 2 8 16 11 11 2 12 15
          1 1 1  1  1  1 1  1   1];
% cedges = load('/Users/davidmur/papers/control_alg/mathematica/SolForMatlab.txt');      
[cand, edges] = size(cedges);
results = zeros(cand-2,1);
for j = 3:cand;
    FF1 = F3;
%     disp(j);
    for i = 1:edges
        if(cedges(j,i)==1)
            FF1 = TruthTable_del_a_temp(FF1,nv,varF,p,cedges(1,i),cedges(2,i),0);
        end
    end
    % for basin size
    [AA,Avec] = bnAsparse(FF1,varF,nv);
    % show te atractos and their respective basins
    [ab,d] = bnAttractor(Avec);
    % ab : specifies the atractors (negative numbers), and their respective
    %      basins elements( positive numbers)
    % d : istance to the atractor
    attrs =  unique(ab(ab<0));
%     disp('num of attrs');
%     disp(length(attrs))
    for k = 1:length(attrs)
        ii = find(ab==-k);
        if (length(ii)==1)
%           disp(length(ii));
            if (ii == 58772)
                basin = length(find(ab==k))+1;
                results(j-2) = 100*basin/2^16; 
%                 disp(100*basin/2^16);
%                 disp('cell death');
            end
        end
    end
%     disp('new round')
end
disp(results)
