Posted By

croftmedia on 03/16/13


Tagged

matlab complex simplicial homology


Versions (?)

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.

  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  

You need to login to post a comment.