Posted By

Tagged

Versions (?)

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

Posted By: junaidasmat on July 26, 2017: Hy, this matlab code is not working. How we can run this code. If we have to run functions one by one or this is the full code.Kindly help me with it ASAP

Comment:

You need to login to post a comment.

Posted By

Tagged

Versions (?)

Simplicial Homology

Related snippets

Comments

Choose a language for easy browsing: