/ 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.
Expand |
Embed | Plain Text
Copy this code and paste it in your HTML
%% %% BASE.M %% function ending = base(simplices) sSize = sSize(2); out = {}; % Cycle through each simplex. % Need to also go through with each vertex removed. out = [out,total_split(simplex)]; end % Now remove permutations of our simplices. output = {}; sSize = sSize(2); % Cycle through each simplex % Go through each member of output, check this isn't a % permutation of it. thisSize = {}; oSize = oSize(2); for k=1:oSize end end oSize = oSize(2); perm = 0; 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'); end function truth = permutation_of(a,b) % Tests if a is a permutation of b % Not the same length, cannot be a permutation truth = 0; return else % Same length ok = 0; 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 = {}; end end end function outputtage = total_split(simplex) outputtage = split(simplex); % Run again on each split outputtage = [outputtage, total_split(... end end end end else else % Now things get interesting end end end end %% %% HOMOLOGY.M %% function groups = homology(basis) sSize = sSize(2); max = 0; end end % Given an empty basis groups = ''; else groups = {}; end end end function sing = sing_homology(basis,k) sSize = sSize(2); max = 0; 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 % It is 1x1. null = mk_0; else end 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 % It is 1x1. m_1 = mk_1; else m_1 = smith(mk_1); end end betti = null-ran; sing = ''; if isempty(sing) sing = 'Z'; else sing = [sing,'(+)Z']; end end if ran~=0 % There are groups to add. m_1 = sym(m_1); maple('with','LinearAlgebra'); m_1 = maple('Diagonal',m_1); end sSize = sSize(1); % We have a group to add. 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 columns = sSize(2); rows = sSize(1); count = 0; % Go through each row ok = 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. columns = sSize(2); rows = sSize(1); counter = 0; % Go through each column. ok = 1; ok = 0; break end end if ok==1 counter = counter + 1; end end end function elem = chain(basis,k) elem = {}; sSize = sSize(2); end end end function N = bound(basis,k) kelem = chain(basis,k-1); kpluselem = chain(basis,k); m = m(2); n = n(2); if k==0 % It's the trivial map. elseif isempty(kpluselem) % It's the highest map. else % We're just dealing with vertices and edges % Cycle through each k+1-simplex % Cycle through each k-simplex end end end end function boundary_val = boundary_entry(a,b) ok = 0; 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. ok = 0; 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 = ''; if i~=notin_index end end % Now calculate number of transpositions to get from contsimp to a! counter = 0; while contsimp~=a % We've got a boundering vertex! % Find the otherone equalling a(i) counter = counter + 1; break end end end end end 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); end