## Posted By

croftmedia on 03/16/13

# 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.

`%%%% 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    endend function outage = split(simplex)    outage = {};    for j=1:length(simplex)        for k=j:length(simplex)            outage{end+1} = simplex(j:k);        end    endend 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    endend 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    endend %%%% 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    endend 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    endend 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    endend 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    endend 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    endend 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    endend 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;    endend function D=smith(A);    A = sym(A);    maple('with','LinearAlgebra');    D = maple('SmithForm',A);    D = double(D);end` Subscribe to comments