Revision: 62844
Initial Code
Initial URL
Initial Description
Initial Title
Initial Tags
Initial Language
at March 16, 2013 02:22 by croftmedia
Initial Code
%% %% 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
Initial URL
Initial Description
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.
Initial Title
Simplicial Homology
Initial Tags
Initial Language
MatLab