# Posted By

croftmedia on 03/16/13

# Statistics

Viewed 59 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.

Copy this code and paste it in your HTML
`%%%% 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