Posted By

croftmedia on 03/16/13

Tagged

Versions (?)

Last Edited at 03/16/13 02:22am

Statistics

Viewed 360 times

Favorited by 0 user(s)

Simplicial Homology

/ Published in: MatLab

Simple program to compute simplicial homology groups of a simplicial complex. Includes "base.m", makes a collection of simplices into a simplicial complex. Written by Paul Horrocks as part of final year undergraduate mathematics project on Computational Homology. Requires Maple Kernel installed to run smith normal form operation.

Highlighted TextExpand

%%
%% BASE.M
%%

function ending = base(simplices)
    sSize = size(simplices);
    sSize = sSize(2);
    out = {};
    for i=1:sSize
        % Cycle through each simplex.
        simplex = char(simplices(i));
        % Need to also go through with each vertex removed.
        out = [out,total_split(simplex)];
    end
    % Now remove permutations of our simplices.
    output = {};
    out = unique(out);
    sSize = size(out);
    sSize = sSize(2);
    for i=1:sSize
       % Cycle through each simplex
       simplex = char(out(i));
       % Go through each member of output, check this isn't a
       % permutation of it.
       thisSize = {};
       oSize = size(output);
       oSize = oSize(2);
       for k=1:oSize
           if length(char(output(k)))==length(simplex)
               thisSize{end+1} = char(output(k));
           end
       end
       oSize = size(thisSize);
       oSize = oSize(2);
       perm = 0;
       for j=1:oSize
           if permutation_of(char(simplex),char(thisSize(j)))==1
               perm = 1;
               break
           end
       end
       if perm==0
           % Not a permutation, add it to output!
           output{end+1} = simplex;
       end
    end
    overall = permutation_of('ba','bc');
    ending = unique(output);
end

function truth = permutation_of(a,b)
    % Tests if a is a permutation of b
    if length(a)~=length(b)
        % Not the same length, cannot be a permutation
        truth = 0;
        return
    else
        % Same length
        for i=1:length(a)
            ok = 0;
            for j=1:length(b)
                if a(i)==b(j)
                    ok = 1;
                    break
                end
            end
            if ok==0
                break
            end
        end
        if ok==1
            % Found all the characters
            truth = 1;
        else
            truth = 0;
        end
    end
end

function outage = split(simplex)
    outage = {};
    for j=1:length(simplex)
        for k=j:length(simplex)
            outage{end+1} = simplex(j:k);
        end
    end
end

function outputtage = total_split(simplex)
    outputtage = split(simplex);
    if length(char(simplex))>1
        % Run again on each split
        for i=1:length(simplex)
            outputtage = [outputtage, split(sub_str(simplex,i+1,i-1))];
            outputtage = [outputtage, total_split(... 
            sub_str(simplex,i+1,i-1))];
        end
    end
end

function sub = sub_str(stringy,i,j)
    if j==0
        j=length(stringy);
    end
    if i==j
        sub = stringy(i);
    else
        if i<j
            sub = stringy(i:j);
        else
            if i>j
                % Now things get interesting
                sub = [stringy(i:length(stringy)) stringy(1:j)];
            end
        end
    end
end

%%
%% HOMOLOGY.M
%%

function groups = homology(basis)
    sSize = size(basis);
    sSize = sSize(2);
    max = 0;
    for i=1:sSize
        if length(char(basis(i)))>max
            max = length(char(basis(i)));
        end
    end
    if max==0
        % Given an empty basis
        groups = '';
    else
        groups = {};
        for i=0:max
            groups{end+1} = ['H_',[num2str(i)],['=',...
            [ sing_homology(basis,i)]]];
        end
    end
end

function sing = sing_homology(basis,k)
    sSize = size(basis);
    sSize = sSize(2);
    max = 0;
    for i=1:sSize
        if length(char(basis(i)))>max
            max = length(char(basis(i)));
        end
    end
    max = max-1;
    if k<0
        sing = '0';
    elseif k>max
        % We're too high.
        sing = '0';
    else
        % Find nullity of d_k
        chain(basis,k-1);
        chain(basis,k);
        mk_0 = bound(basis,k);
        % Need to check it's not a 1x1 matrix
        sSize = size(mk_0);
        if all(sSize(1)==1 & sSize(2)==1)
            % It is 1x1.
            null = mk_0;
        else
            null = smith(mk_0);
        end
        null = kernel(null);
        if k==max
            % Don't find range.
            ran=0;
        else
            % Find range of map
            mk_1 = bound(basis,k+1);
            % Need to check it's not a 1x1 matrix
            sSize = size(mk_1);
            if all(sSize(1)==1 & sSize(2)==1)
                % It is 1x1.
                m_1 = mk_1;
            else
                m_1 = smith(mk_1);
            end
            ran = image(m_1);
        end
        betti = null-ran;
        sing = '';
        for i=1:betti
            if isempty(sing)
                sing = 'Z';
            else
                sing = [sing,'(+)Z'];
            end
        end
        if ran~=0
            % There are groups to add.
            sSize = size(m_1);
            if ~all(sSize(1)==1 & sSize(2)==1)
                m_1 = sym(m_1);
                maple('with','LinearAlgebra');
                m_1 = maple('Diagonal',m_1);
                m_1 = double(m_1);
            end
            sSize = size(m_1);
            sSize = sSize(1);
            for i=1:sSize
                m_1 = abs(m_1);
                if m_1(i)~=0
                    % We have a group to add.
                    if m_1(i)~=1
                        group = ['Z_',num2str(m_1(i))];
                        if isempty(sing)
                            sing = group;
                        else
                            sing = [sing,['(+)',group]];
                        end
                    end
                end
            end
        end
        if k==0
            % Add one Z to it.
            if isempty(sing)
                sing = 'Z';
            else
                sing = [sing,'(+)Z'];
            end
        end
        if isempty(sing)
            sing = '0';
        end
    end
end

function count = image(matrix)
    sSize = size(matrix);
    columns = sSize(2);
    rows = sSize(1);
    count = 0;
    for i=1:rows
        % Go through each row
        ok = 0;
        for j=1:columns
            if matrix(i,j)~=0
                ok = 1;
            end
        end
        if ok==1
            count = count+1;
        end
    end
end

function counter = kernel(matrix)
    % Takes nxm matrix, counts columns that are all zero.
    sSize = size(matrix);
    columns = sSize(2);
    rows = sSize(1);
    counter = 0;
    for i=1:columns
        % Go through each column.
        ok = 1;
        for j=1:rows
            if matrix(j,i)~=0
                ok = 0;
                break
            end
        end
        if ok==1
            counter = counter + 1;
        end
    end
end

function elem = chain(basis,k)
    elem = {};
    sSize = size(basis);
    sSize = sSize(2);
    for i=1:sSize
        if length(char(basis(i)))==(k+1)
             elem{end+1} = char(basis(i));
        end
    end
end

function N = bound(basis,k)
    kelem = chain(basis,k-1);
    kpluselem = chain(basis,k);
    m = size(kelem);
    m = m(2);
    n = size(kpluselem);
    n = n(2);
    if k==0
        % It's the trivial map.
        N = ones(1,n);
    elseif isempty(kpluselem)
        % It's the highest map.
        N = zeros(m,1);
    else
        N = zeros(m,n);
        % We're just dealing with vertices and edges
        for i=1:m
            % Cycle through each k+1-simplex
            ksimp = kelem(i);
            for j=1:n
                % Cycle through each k-simplex
                kplussimp = kpluselem(j);
                N(i,j) = boundary_entry(ksimp,kplussimp);
            end
        end
    end
end

function boundary_val = boundary_entry(a,b)
    a = char(a);
    b = char(b);
    for i=1:length(a)
        ok = 0;
        for j=1:length(b)
            if a(i)==b(j)
                ok = 1;
                break
            end
        end
        if ok==0
            break
        end
    end
    if ok==1
        % Found all the characters
        % Work out which character is NOT in the simplex.
        for i=1:length(b)
            ok = 0;
            for j=1:length(a)
                if b(i)==a(j)
                    ok = 1;
                    break
                end
            end
            if ok==0
                % This character not in a.
                notin_index = i;
            end
        end
        % Create subcomplex without the i'th vertex
        contsimp = '';
        for i=1:length(b)
            if i~=notin_index
                contsimp = [contsimp,b(i)];
            end
        end
        % Now calculate number of transpositions to get from contsimp to a!
        counter = 0;
        while contsimp~=a
            for i=1:length(a)
                if contsimp(i)~=a(i)
                    % We've got a boundering vertex!
                    for j=(i+1):length(a)
                        % Find the otherone equalling a(i)
                        counter = counter + 1;
                        if contsimp(j)==a(i)
                            contsimp(j) = contsimp(i);
                            contsimp(i) = a(i);
                            break
                        end
                    end
                end
            end
        end
        if mod(counter,2)==0
            counter = 1;
        else
            counter = -1;
        end
        boundary_val = (-1)^(notin_index-1)*counter;
    else
        boundary_val = 0;
    end
end

function D=smith(A);
    A = sym(A);
    maple('with','LinearAlgebra');
    D = maple('SmithForm',A);
    D = double(D);
end

Report this snippet Tweet

Comments

Subscribe to comments

Comment:

You need to login to post a comment.