%PHdemo -- persistent homology algorithm (ELZ/ZC): visual demo
%
% [2011-may-02] Vin de Silva, Pomona College.

%% data

% abstract complex:
S = {[1], [2], [3], [4], [1 2], [2 3], [3 4], [1 4], [1 3], [1 3 4]};

% vertex positions:
X = [0 0; 0 1; 1 1; 1 0]';

%% runtime parameters

step_pause = inf;     % pause this many seconds; Inf = wait for mouseclick
substep_pause = inf;  % pause this many seconds; Inf = wait for mouseclick

explore = 1;    % after calculation, enter explore mode? 1=yes, 0=no.

%% graphical parameters

% colors
figure_col = [0.8 0.8 0.8];

% binary colors
bitoff_col = [0.7 0.7 0.7];
bit0_col = [0 0 0];
bit1_col = [1 1 1];

% leading term colors (chain/cycle/boundary)
C_col = [1 0 0];
Z_col = [0 1 0];
B_col = [1 1 0];

C_col2 = [1, 0.4, 0.4];
Z_col2 = [0.4, 1, 0.4];
B_col2 = [1 1 0.4];

% grid colors (Current/Buffer/plain)
current_col = C_col;
buffer_col = B_col;
grid_col = [0.7 0.7 0.7];

% shape colors
caret_col = [0 0 0];
trail_col = [0.6 0.6 0.6];

% simplicial complex colours
d_col_off = [0.6 0.6 0.6; 0.65 0.65 0.65; 0.7 0.7 0.7];
d_col_on = [0.2 0.2 0.2; 0.3 0.3 0.3; 0.4 0.4 0.4];


% shapes
gridX = [0 0 1 1 0];
gridY = [0 1 1 0 0];

bitX = [0.1 0.1 0.9 0.9 0.1];
bitY = [0.1 0.9 0.9 0.1 0.1];

caretX = [0.1 0.5 0.9 0.5 0.1];
caretY = [0.5 0.1 0.5 0.9 0.5];



%% organise the data, compute the boundary matrix D

n = length(S);

dS = zeros(n,1); % dimension of each cell
for j = 1:n
    S{j} = sort(S{j});
    dS(j) = length(S{j})-1;  
end

Smat = zeros(n,max(dS)+1);
for j = 1:n
    Smat(j, 1:(dS(j)+1)) = S{j};
end

% calculate boundary matrix
D = zeros(n);

for j = find(dS(:) > 0)'    
    for k = 0:dS(j)
        Sjk = S{j}([1:k, k+2:end]); % the kth face of simplex j
        found = find(ismember([dS, Smat(:, 1:dS(j))], [dS(j)-1, Sjk], 'rows'));
        if length(found)==1
            D(j, found) = (-1)^k;
        elseif length(found) > 1
            error('S has repeated cells.')
        elseif (length(found) < 1)
            % no problem here
        end
    end
end

if any(any(D^2))
    error('Must have D^2 = 0.')
end

%% create matrix plot, store handles
figure(47), clf
set(47, 'position', [10 100 800 700])

HaMain = axes('position', [0.05 0.05 0.95 0.95]);
axis equal off
axis([-2 n+3 0 n+2])

% coefficient matrix
HMbg = zeros(n);
HMfg = zeros(n);
HMtxt = zeros(1,n);

for j = 1:n
    for k = 1:j
        HMbg(j,k) = patch(j+gridX, k+gridY, 'w');
        set(HMbg(j,k), 'facecolor', grid_col);

        HMfg(j,k) = patch(j+bitX, k+bitY, 'k');
        set(HMfg(j,k), 'facecolor', bitoff_col);
    end

    jlabel = sprintf(' %d', S{j});
    jlabel = jlabel(2:end);
    HMtxt(j) = text(j+0.5,0.7, jlabel);
    set(HMtxt(j), 'horizontalalignment', 'right')
    set(HMtxt(j), 'rotation', 90)

end

% buffer (boundary of current chain)
HBbg = zeros(1,n);
HBfg = zeros(1,n);
HLtxt = zeros(1,n); % left labels
HRtxt = zeros(1,n); % right labels

for k = 1:n
    j = -1;
    
    HBbg(k) = patch(j+gridX, k+gridY, 'k');
    set(HBbg(k), 'facecolor', buffer_col);
    HBfg(k) = patch(j+bitX, k+bitY, 'k');
    set(HBfg(k), 'facecolor', bit0_col);
    
    klabel = sprintf(' %d', S{k});
    klabel = klabel(2:end);
    HLtxt(k) = text(-1.3, k+0.5, klabel);
    set(HLtxt(k), 'horizontalalignment', 'right')
    HRtxt(k) = text(n+1.3, k+0.5, klabel);
    set(HRtxt(k), 'horizontalalignment', 'left')

end

%HBtxt0 = text(-0.5, 0.5, '\partial', 'fontsize', 16);

%% create geometric plot with figure handles

figure(48), clf
set(48, 'position', [820 400 500 400])
HaPlot = axes('position', [0.1 0.1 0.8 0.8]);
axis equal off
hold on

HS = zeros(1,n);    % handles for the cells of X
pHS = cell(1,n);    % property type for the color of the cell

% 2-cells
cells = find(dS == 2);
for a = 1: length(cells)
   ca = cells(a); 
   HS(ca) = patch(X(1,Smat(ca,[1 2 3])), X(2,Smat(ca,[1 2 3])), 'k');
   set(HS(ca), 'linestyle', 'none', 'facecolor', d_col_off(3,:));
   pHS{ca} = 'facecolor';
end

% 1-cells
cells = find(dS == 1);
for a = 1: length(cells)
   ca = cells(a);
   HS(ca) = plot(X(1,Smat(ca,[1 2])), X(2,Smat(ca,[1 2])), 'k-');
   set(HS(ca), 'linewidth', 4, 'color', d_col_off(2,:));
   pHS{ca} = 'color';
end

% 0-cells
cells = find(dS == 0);
for a = 1: length(cells)
    ca = cells(a);

    HS(ca) = plot(X(1,Smat(ca,[1])), X(2,Smat(ca,[1])), 'k.');
    set(HS(ca), 'markersize', 48, 'color', d_col_off(1,:));
    pHS{ca} = 'color';

    HV(a) = text(X(1,Smat(ca,[1])), X(2,Smat(ca,[1])), sprintf('%d', a));
    set(HV(a), 'horizontalalignment', 'center', 'color', [0.5 0.5 0.5])
end


%rescale axes
ax48 = axis;
ax48mX = mean(ax48([1 2]));
ax48mY = mean(ax48([3 4]));
newax48 = (1.2)*ax48 - (0.2)*[ax48mX ax48mX ax48mY ax48mY];
axis(newax48)

clear ax48 ax48mX ax48mY newax48


%% Main Loop

% work modulo 2
D2 = (D~=0);

% output matrices
E = D2;

status = zeros(1,n);
 % 1: cycle
 % 2: boundary
 % 3: coboundary (cycle-killer)

partner = zeros(1,n);

Hpd = zeros(1,n);
Htrail = zeros(1,n);
% handles for the persistent diagram carets and lines

% iterate over cells
for j = 1: n
    % highlight the current column
    for k = 1: j
        set(HMbg(j,k), 'facecolor', current_col);
        set(HMfg(j,k), 'facecolor', bit0_col);
    end
    
    % place current chain in buffer C
    C = zeros(1,n);
    C(j) = 1;
    
    % place boundary of C in buffer B
    B = mod((C * D2),2);
    
    %-----------------------------------------------
    %% recursive calculation
    reduced = 0;
    while ~reduced
        % update algebraic display of C and B
        set(HMfg(j,j), 'facecolor', [1 1 1]);
        for k = 1:j-1
            set(HMfg(j,k), 'facecolor', bit1_col*C(k));        
            set(HBfg(k), 'facecolor', bit1_col*B(k));
        end
        
        % update geometric display of C and B
        for k = 1:j
           if C(k)
               set(HS(k), pHS{k}, C_col);
           elseif B(k)
               set(HS(k), pHS{k}, B_col);
           else
               set(HS(k), pHS{k}, d_col_on(dS(k)+1, :));
           end
        end
        
        drawnow
        if isinf(substep_pause)
            waitforbuttonpress
        else
            pause(substep_pause)
        end

        % three cases
        if all(B == 0)              % C is a cycle
            
            E(j,:) = C;
            status(j) = 1;
            partner(j) = Inf;
            set(HMfg(j,j), 'facecolor', Z_col);
            
            reduced = 1;
            continue
        end
        
        lead = max(find(B));
        if (status(lead) == 1)      % C spans a nonbounding cycle

            % store and display C
            E(j,:) = C;
            status(j) = 3;
            partner(j) = lead;
            set(HMfg(j,j), 'facecolor', C_col);

            % store and display B
            E(lead,:) = B;
            status(lead) = 2;
            partner(lead) = j;
            set(HMfg(lead,lead), 'facecolor', B_col);
            for k = 1: lead-1
                set(HMfg(lead,k), 'facecolor', B(k)*bit1_col);
            end
               
            % update persistence diagram 
            axes(HaMain)
            Htrail(lead) = line([lead+0.5, lead+0.5], [lead+1.05, j+0.5]);
            set(Htrail(lead), 'linewidth', 2, 'color', trail_col);            

            Htrail(j) = line([lead+0.5, j-0.05], [j+0.5, j+0.5]);
            set(Htrail(j), 'linewidth', 2, 'color', trail_col);            

            Hpd(lead) = patch(lead+caretX, j+caretY, caret_col);
              
            reduced = 1;
                                
        elseif (status(lead) == 2)  % leading term of dC is a boundary
            colead = partner(lead);
            C = C + E(colead,:);           
            C = mod(C,2);
             
            B = B + E(lead,:);
            B = mod(B,2);
        end

    end
    %-----------------------------------------------
          
    % clean up to finish iteration
    drawnow
    if isinf(step_pause)
        waitforbuttonpress
    else
        pause(step_pause)
    end
    
    
    % clear buffer B
    for k = 1:n
        set(HBfg(k), 'facecolor', bit0_col);
    end
    
    % unhighlight the current column
    for k = 1: j
        set(HMbg(j,k), 'facecolor', grid_col);
    end
    
    % uncolor geometric complex
    for k = (1: j)
        set(HS(k), pHS{k}, d_col_on(dS(k)+1, :));
    end
    
%     drawnow
%     if isinf(step_pause)
%         waitforbuttonpress
%     else
%         pause(step_pause)
%     end
end

% persistent diagram: points at infinity
for l = find(status==1)
    Htrail(l) = line([l+0.5, l+0.5], [l+1.05, n+1+0.5]);
    set(Htrail(l), 'linewidth', 2, 'color', trail_col);

    Hpd(l) = patch(l+caretX, n+1+caretY, caret_col); 
end
drawnow



%% interactive mode
if (explore)
    set(47, 'name', 'mouse to explore');
    set(48, 'name', 'keyboard to quit')
    
    % enter interactive mode
    keydown = 0;
    while ~keydown
       keydown = waitforbuttonpress;
       Hco = gco;
       
       if ismember(Hco, Hpd)
           % selected persistence diagram caret
           k = find(Hpd == Hco);
           j = partner(k);
           for a = (1: n)
               if (a > j)
                   set(HS(a), pHS{a}, d_col_off(dS(a)+1, :));
               elseif (E(k,a) & ~isinf(j))
                   set(HS(a), pHS{a}, B_col);
               elseif (E(k,a) & isinf(j))
                   set(HS(a), pHS{a}, Z_col);                   
               elseif (~isinf(j) & E(j,a))
                   set(HS(a), pHS{a}, C_col);
               else
                   set(HS(a), pHS{a}, d_col_on(dS(a)+1, :));
               end    
           end           
           
       elseif ismember(Hco, HMfg(1:n+1:n^2))
           % selected leading term of chain/cycle/boundary
           k = find(HMfg(1:n+1:n^2) == Hco);
           switch status(k)
               case 1
                   Kcol  = Z_col;
                   Kcol2 = Z_col2;
               case 2
                   Kcol  = B_col;
                   Kcol2 = B_col2;
               case 3
                   Kcol  = C_col;
                   Kcol2 = C_col2;
           end
           
           for a = (1: n)
               if (a > k)
                   set(HS(a), pHS{a}, d_col_off(dS(a)+1,:));               
               elseif (a == k)
%                   set(HS(a), pHS{a}, get(Hco, 'facecolor'));
                   set(HS(a), pHS{a}, Kcol);
               elseif E(k,a)
%                   set(HS(a), pHS{a}, bit1_col);
                   set(HS(a), pHS{a}, Kcol2);
               else
                   set(HS(a), pHS{a}, d_col_on(dS(a)+1,:));
              end
               
           end
       
       else
           % clear display
           for a = (1: n)
                set(HS(a), pHS{a}, d_col_on(dS(a)+1, :));
           end
           drawnow
           
       end        
    end    
end




