Posted By


croftmedia on 03/16/13

Tagged


Statistics


Viewed 352 times
Favorited by 0 user(s)

Simplicial Homology


/ Published in: MatLab
Save to your folder(s)

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
  1. %%
  2. %% BASE.M
  3. %%
  4.  
  5.  
  6. function ending = base(simplices)
  7. sSize = size(simplices);
  8. sSize = sSize(2);
  9. out = {};
  10. for i=1:sSize
  11. % Cycle through each simplex.
  12. simplex = char(simplices(i));
  13. % Need to also go through with each vertex removed.
  14. out = [out,total_split(simplex)];
  15. end
  16. % Now remove permutations of our simplices.
  17. output = {};
  18. out = unique(out);
  19. sSize = size(out);
  20. sSize = sSize(2);
  21. for i=1:sSize
  22. % Cycle through each simplex
  23. simplex = char(out(i));
  24. % Go through each member of output, check this isn't a
  25. % permutation of it.
  26. thisSize = {};
  27. oSize = size(output);
  28. oSize = oSize(2);
  29. for k=1:oSize
  30. if length(char(output(k)))==length(simplex)
  31. thisSize{end+1} = char(output(k));
  32. end
  33. end
  34. oSize = size(thisSize);
  35. oSize = oSize(2);
  36. perm = 0;
  37. for j=1:oSize
  38. if permutation_of(char(simplex),char(thisSize(j)))==1
  39. perm = 1;
  40. break
  41. end
  42. end
  43. if perm==0
  44. % Not a permutation, add it to output!
  45. output{end+1} = simplex;
  46. end
  47. end
  48. overall = permutation_of('ba','bc');
  49. ending = unique(output);
  50. end
  51.  
  52. function truth = permutation_of(a,b)
  53. % Tests if a is a permutation of b
  54. if length(a)~=length(b)
  55. % Not the same length, cannot be a permutation
  56. truth = 0;
  57. return
  58. else
  59. % Same length
  60. for i=1:length(a)
  61. ok = 0;
  62. for j=1:length(b)
  63. if a(i)==b(j)
  64. ok = 1;
  65. break
  66. end
  67. end
  68. if ok==0
  69. break
  70. end
  71. end
  72. if ok==1
  73. % Found all the characters
  74. truth = 1;
  75. else
  76. truth = 0;
  77. end
  78. end
  79. end
  80.  
  81. function outage = split(simplex)
  82. outage = {};
  83. for j=1:length(simplex)
  84. for k=j:length(simplex)
  85. outage{end+1} = simplex(j:k);
  86. end
  87. end
  88. end
  89.  
  90. function outputtage = total_split(simplex)
  91. outputtage = split(simplex);
  92. if length(char(simplex))>1
  93. % Run again on each split
  94. for i=1:length(simplex)
  95. outputtage = [outputtage, split(sub_str(simplex,i+1,i-1))];
  96. outputtage = [outputtage, total_split(...
  97. sub_str(simplex,i+1,i-1))];
  98. end
  99. end
  100. end
  101.  
  102. function sub = sub_str(stringy,i,j)
  103. if j==0
  104. j=length(stringy);
  105. end
  106. if i==j
  107. sub = stringy(i);
  108. else
  109. if i<j
  110. sub = stringy(i:j);
  111. else
  112. if i>j
  113. % Now things get interesting
  114. sub = [stringy(i:length(stringy)) stringy(1:j)];
  115. end
  116. end
  117. end
  118. end
  119.  
  120. %%
  121. %% HOMOLOGY.M
  122. %%
  123.  
  124. function groups = homology(basis)
  125. sSize = size(basis);
  126. sSize = sSize(2);
  127. max = 0;
  128. for i=1:sSize
  129. if length(char(basis(i)))>max
  130. max = length(char(basis(i)));
  131. end
  132. end
  133. if max==0
  134. % Given an empty basis
  135. groups = '';
  136. else
  137. groups = {};
  138. for i=0:max
  139. groups{end+1} = ['H_',[num2str(i)],['=',...
  140. [ sing_homology(basis,i)]]];
  141. end
  142. end
  143. end
  144.  
  145. function sing = sing_homology(basis,k)
  146. sSize = size(basis);
  147. sSize = sSize(2);
  148. max = 0;
  149. for i=1:sSize
  150. if length(char(basis(i)))>max
  151. max = length(char(basis(i)));
  152. end
  153. end
  154. max = max-1;
  155. if k<0
  156. sing = '0';
  157. elseif k>max
  158. % We're too high.
  159. sing = '0';
  160. else
  161. % Find nullity of d_k
  162. chain(basis,k-1);
  163. chain(basis,k);
  164. mk_0 = bound(basis,k);
  165. % Need to check it's not a 1x1 matrix
  166. sSize = size(mk_0);
  167. if all(sSize(1)==1 & sSize(2)==1)
  168. % It is 1x1.
  169. null = mk_0;
  170. else
  171. null = smith(mk_0);
  172. end
  173. null = kernel(null);
  174. if k==max
  175. % Don't find range.
  176. ran=0;
  177. else
  178. % Find range of map
  179. mk_1 = bound(basis,k+1);
  180. % Need to check it's not a 1x1 matrix
  181. sSize = size(mk_1);
  182. if all(sSize(1)==1 & sSize(2)==1)
  183. % It is 1x1.
  184. m_1 = mk_1;
  185. else
  186. m_1 = smith(mk_1);
  187. end
  188. ran = image(m_1);
  189. end
  190. betti = null-ran;
  191. sing = '';
  192. for i=1:betti
  193. if isempty(sing)
  194. sing = 'Z';
  195. else
  196. sing = [sing,'(+)Z'];
  197. end
  198. end
  199. if ran~=0
  200. % There are groups to add.
  201. sSize = size(m_1);
  202. if ~all(sSize(1)==1 & sSize(2)==1)
  203. m_1 = sym(m_1);
  204. maple('with','LinearAlgebra');
  205. m_1 = maple('Diagonal',m_1);
  206. m_1 = double(m_1);
  207. end
  208. sSize = size(m_1);
  209. sSize = sSize(1);
  210. for i=1:sSize
  211. m_1 = abs(m_1);
  212. if m_1(i)~=0
  213. % We have a group to add.
  214. if m_1(i)~=1
  215. group = ['Z_',num2str(m_1(i))];
  216. if isempty(sing)
  217. sing = group;
  218. else
  219. sing = [sing,['(+)',group]];
  220. end
  221. end
  222. end
  223. end
  224. end
  225. if k==0
  226. % Add one Z to it.
  227. if isempty(sing)
  228. sing = 'Z';
  229. else
  230. sing = [sing,'(+)Z'];
  231. end
  232. end
  233. if isempty(sing)
  234. sing = '0';
  235. end
  236. end
  237. end
  238.  
  239. function count = image(matrix)
  240. sSize = size(matrix);
  241. columns = sSize(2);
  242. rows = sSize(1);
  243. count = 0;
  244. for i=1:rows
  245. % Go through each row
  246. ok = 0;
  247. for j=1:columns
  248. if matrix(i,j)~=0
  249. ok = 1;
  250. end
  251. end
  252. if ok==1
  253. count = count+1;
  254. end
  255. end
  256. end
  257.  
  258. function counter = kernel(matrix)
  259. % Takes nxm matrix, counts columns that are all zero.
  260. sSize = size(matrix);
  261. columns = sSize(2);
  262. rows = sSize(1);
  263. counter = 0;
  264. for i=1:columns
  265. % Go through each column.
  266. ok = 1;
  267. for j=1:rows
  268. if matrix(j,i)~=0
  269. ok = 0;
  270. break
  271. end
  272. end
  273. if ok==1
  274. counter = counter + 1;
  275. end
  276. end
  277. end
  278.  
  279. function elem = chain(basis,k)
  280. elem = {};
  281. sSize = size(basis);
  282. sSize = sSize(2);
  283. for i=1:sSize
  284. if length(char(basis(i)))==(k+1)
  285. elem{end+1} = char(basis(i));
  286. end
  287. end
  288. end
  289.  
  290. function N = bound(basis,k)
  291. kelem = chain(basis,k-1);
  292. kpluselem = chain(basis,k);
  293. m = size(kelem);
  294. m = m(2);
  295. n = size(kpluselem);
  296. n = n(2);
  297. if k==0
  298. % It's the trivial map.
  299. N = ones(1,n);
  300. elseif isempty(kpluselem)
  301. % It's the highest map.
  302. N = zeros(m,1);
  303. else
  304. N = zeros(m,n);
  305. % We're just dealing with vertices and edges
  306. for i=1:m
  307. % Cycle through each k+1-simplex
  308. ksimp = kelem(i);
  309. for j=1:n
  310. % Cycle through each k-simplex
  311. kplussimp = kpluselem(j);
  312. N(i,j) = boundary_entry(ksimp,kplussimp);
  313. end
  314. end
  315. end
  316. end
  317.  
  318. function boundary_val = boundary_entry(a,b)
  319. a = char(a);
  320. b = char(b);
  321. for i=1:length(a)
  322. ok = 0;
  323. for j=1:length(b)
  324. if a(i)==b(j)
  325. ok = 1;
  326. break
  327. end
  328. end
  329. if ok==0
  330. break
  331. end
  332. end
  333. if ok==1
  334. % Found all the characters
  335. % Work out which character is NOT in the simplex.
  336. for i=1:length(b)
  337. ok = 0;
  338. for j=1:length(a)
  339. if b(i)==a(j)
  340. ok = 1;
  341. break
  342. end
  343. end
  344. if ok==0
  345. % This character not in a.
  346. notin_index = i;
  347. end
  348. end
  349. % Create subcomplex without the i'th vertex
  350. contsimp = '';
  351. for i=1:length(b)
  352. if i~=notin_index
  353. contsimp = [contsimp,b(i)];
  354. end
  355. end
  356. % Now calculate number of transpositions to get from contsimp to a!
  357. counter = 0;
  358. while contsimp~=a
  359. for i=1:length(a)
  360. if contsimp(i)~=a(i)
  361. % We've got a boundering vertex!
  362. for j=(i+1):length(a)
  363. % Find the otherone equalling a(i)
  364. counter = counter + 1;
  365. if contsimp(j)==a(i)
  366. contsimp(j) = contsimp(i);
  367. contsimp(i) = a(i);
  368. break
  369. end
  370. end
  371. end
  372. end
  373. end
  374. if mod(counter,2)==0
  375. counter = 1;
  376. else
  377. counter = -1;
  378. end
  379. boundary_val = (-1)^(notin_index-1)*counter;
  380. else
  381. boundary_val = 0;
  382. end
  383. end
  384.  
  385. function D=smith(A);
  386. A = sym(A);
  387. maple('with','LinearAlgebra');
  388. D = maple('SmithForm',A);
  389. D = double(D);
  390. end

Report this snippet


Comments

RSS Icon Subscribe to comments

You need to login to post a comment.