• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Noncrossing partitions up to rotation

mate2code Feb 13th, 2013 212 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
2. % This MATLAB code can be used to calculate the rows of the number triangle OEIS A209805,
3. % i.e. the number of k-block noncrossing partitions of an n-element set up to rotation.
4. % This is a brute force method and may cause memory problems for n > 12.
5. % All objects are created and than counted,
6. % so this code can also be used to find the actual partitions.
7. % Written by Tilman Piesk, March 2012.
8. %
9. %
10. % The functions "SetPartition" and "Bell" by Bruno Luong from Matlab Central are used,
11. % as well as the following two functions:
12.
13.
14. function [output] = crossing(A,B)
15. % Two sets crossing?
16. % Enter sets as row vectors.
17. cross = 0;
18. if min(B) < min(A)
19.     H = A ; A = B ; B = H ;
20. end
21. x2 = min(B) ;
22. if max(A) > x2
23.     x3 = min(  setdiff(A,1:x2)  ) ;
24.     if max(B) > x3
25.         cross = 1 ;
26.     end
27. end
28. output = cross ;
29. end
30.
31.
32. function [y] = perm2mat(x)
33. % Input: Permutation of elements 1 ... n as row vector
34. % Output: Permutation matrix (to permute rows, e.g. elements of a column vector)
35. wide = size(x,2) ;
36. I = eye(wide) ;
37. for k = 1:wide
38.     P( k , 1:wide ) = I( x(k) , : ) ;
39. end
40. y = P ;
41. end
42.
43.
44. % The row number is called "num" and has to be entered before calculation.
45. % To explain the code some outputs for num = 6 have been included as comments.
46. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47.
48.
49. %                            num =
50. %                                 6
51.
52.
53.
54. bell = Bell(num) ;
55.
56. %                            bell =
57. %                               203
58.
59.
60.
61. powtwo = zeros(num,1) ;
62. for m=1:num
63.     powtwo(m) = 2^(m-1) ;
64. end
65.
66. %                            powtwo =
67. %                                 1
68. %                                 2
69. %                                 4
70. %                                 8
71. %                                16
72. %                                32
73.
74.
75.
76. %%%%%  All partitions:
77.
78. part1 = SetPartition(num) ;
79.
80. %                            part1{20}
81. %                            ans =
82. %                                [1x3 double]    [1x3 double]
83. %                            part1{20}{1}
84. %                            ans =
85. %                                 1     2     4
86. %                            part1{20}{2}
87. %                            ans =
88. %                                 3     5     6
89.
90. %                            part1{50}
91. %                            ans =
92. %                                [1x2 double]    [3]    [1x2 double]    [5]
93. %                            part1{50}{1}
94. %                            ans =
95. %                                 1     2
96. %                            part1{50}{2}
97. %                            ans =
98. %                                 3
99. %                            part1{50}{3}
100. %                            ans =
101. %                                 4     6
102. %                            part1{50}{4}
103. %                            ans =
104. %                                 5
105.
106.
107.
108. %%%%%  How many blocks does each partition have?
109.
110. part1long = zeros(bell,1) ;
111. for m=1:bell
112.     part1long(m) = size(part1{m},2) ;
113. end
114.
115. %                            part1long(20)
116. %                            ans =
117. %                                 2
118.
119. %                            part1long(50)
120. %                            ans =
121. %                                 4
122.
123.
124.
125. %%%%%  How many non-singleton blocks does each partition have?
126.
127. partlong = zeros(bell,1) ;
128. for m=1:bell
129.     M=0;
130.     for n=1:part1long(m)
131.         if size(part1{m}{n},2) > 1
132.             M = M+1 ;
133.         end
134.     end
135.     partlong(m) = M ;
136. end
137.
138. %                            partlong(20)
139. %                            ans =
140. %                                 2
141.
142. %                            partlong(50)
143. %                            ans =
144. %                                 2
145.
146.
147.
148. %%%%%  All partitions, singleton blocks are not shown:
149.
150. part = cell(bell,1) ;
151. for m=1:bell
152.     i = 1 ;
153.     for n=1:part1long(m)
154.         N = part1{m}{n} ;
155.         if size(N,2) > 1
156.             part{m}{i} = N ;
157.             i = i+1 ;
158.         end
159.     end
160. end
161.
162. %                            part{50}
163. %                            ans =
164. %                                [1x2 double]    [1x2 double]
165. %                            part{50}{1}
166. %                            ans =
167. %                                 1     2
168. %                            part{50}{2}
169. %                            ans =
170. %                                 4     6
171.
172. clear part1   %  part1 is not needed anymore.
173.
174.
175.
176. %%%%%  Which partitions are crossing?
177.
178. cross = zeros(bell,1) ;
179. for a=1:bell
180.     c = 0 ;
181.     if partlong(a) > 1
182.        A = nchoosek(1:partlong(a),2);   %  index numbers to get all pairs of blocks
183.              for b=1:size(A,1)
184.                  B = A(b,:) ;
185.                  B1 = part{a}{B(1)};
186.                  B2 = part{a}{B(2)};
187.                  c = crossing(B1,B2) ;   %  The function crossing tells whether
188.                  if c == 1               %  two blocks are crossing (1) or not (0).
189.                      break
190.                  end
191.              end
192.     end
193.     cross(a) = c ;
194. end
195.
196. %                            cross(20)
197. %                            ans =
198. %                                 1
199.
200. %                            cross(50)
201. %                            ans =
202. %                                 0
203.
204.
205.
206. %%%%%  The blocks of the partitions shall be represented by binary vectors:
207. %%%%%  E.g. the block {1,2,6} by the vector 110001.
208.
209. matcell = cell(bell,1) ;
210. for m=1:bell
211.     if cross(m) == 0            %  If the partition is crossing matcell(m) is left empty.
212.        high = partlong(m) ;
213.        M = zeros( high , num ) ;
214.        for n=1:high
215.            N = part{m}{n} ;
216.            for p=1:size(N,2)
217.                M( n , N(p) ) = 1 ;
218.            end
219.        end
220.        matcell{m} = M ;
221.     end
222. end
223.
224. %                            matcell{20}
225. %                            ans =
226. %                                 []
227.
228. %                            matcell{50}
229. %                            ans =
230. %                                 1     1     0     0     0     0
231. %                                 0     0     0     1     0     1
232.
233. clear part   %  part is not needed anymore.
234.
235.
236.
237. %%%%%  Turning the partitions by permuting the binary vectors:
238. %%%%%  110001 into 100011, 000111, 001110 etc.
239.
240. %%%%%  First generate permutation matrices:
241.
242. permcell = cell(num,1) ;
243. perm = perm2mat([num 1:(num-1) ]) ;
244. for m=1:num
245.     M = perm^(m-1) ;
246.     permcell{m} = M ;
247. end
248.
249. %                            permcell{2}
250. %                            ans =
251. %                                 0     0     0     0     0     1
252. %                                 1     0     0     0     0     0
253. %                                 0     1     0     0     0     0
254. %                                 0     0     1     0     0     0
255. %                                 0     0     0     1     0     0
256. %                                 0     0     0     0     1     0
257.
258. %                            permcell{3}
259. %                            ans =
260. %                                 0     0     0     0     1     0
261. %                                 0     0     0     0     0     1
262. %                                 1     0     0     0     0     0
263. %                                 0     1     0     0     0     0
264. %                                 0     0     1     0     0     0
265. %                                 0     0     0     1     0     0
266.
267.
268.
269. %%%%%  Permute binary vectors and turn results into decimal numbers,
270. %%%%%  e.g. vector 110000 into number 3.
271.
272. matcellbig = cell(bell,num) ;
273. for m=1:bell
274.     if cross(m) == 0
275.        for n=1:num
276.            matcellbig{m,n} = matcell{m} * permcell{n} * powtwo ;
277.        end
278.     end
279. end
280.
281. %                            matcellbig{50,1}
282. %                            ans =
283. %                                 3
284. %                                40
285. %                            matcellbig{50,2}
286. %                            ans =
287. %                                33
288. %                                20
289. %                            matcellbig{50,3}
290. %                            ans =
291. %                                48
292. %                                10
293. %                            matcellbig{50,4}
294. %                            ans =
295. %                                24
296. %                                 5
297. %                            matcellbig{50,5}
298. %                            ans =
299. %                                12
300. %                                34
301. %                            matcellbig{50,6}
302. %                            ans =
303. %                                 6
304. %                                17
305.
306. clear matcell   %  matcell is not needed anymore.
307.
308.
309.
310. %%%%%  These numbers have to be sorted and put into a matrix,
311. %%%%%  so it can be checked whether two partitions generate the same matrix:
312.
313. matcellsort=cell(bell,1);
314. for m=1:bell
315.     if cross(m) == 0
316.        N = zeros( partlong(m) , num ) ;
317.        for n=1:num
318.            N(1:partlong(m),n) = sort(matcellbig{m,n}) ;
319.        end
320.        matcellsort{m} = sortrows(N') ;
321.     end
322. end
323.
324. %                            matcellsort{50}
325. %                            ans =
326. %                                 3    40
327. %                                 5    24
328. %                                 6    17
329. %                                10    48
330. %                                12    34
331. %                                20    33
332.
333. clear matcellbig   %  matcellbig is not needed anymore.
334.
335.
336.
337. %%%%%  The top row can represent the whole matrix:
338.
339. toprowscell = cell(bell,1) ;
340. for m=1:bell
341.     if cross(m) == 0
342.        M = zeros( 1 , max(partlong) ) ;
343.        M( 1:partlong(m) ) = matcellsort{m}(1,:) ;
344.        toprowscell{m} = [ part1long(m) M ] ;
345.     end
346. end
347.
348. %                            toprowscell{50}
349. %                            ans =
350. %                                 4     3    40     0
351.
352.                                   %  The 4 is the number of blocks in the 50-th partition.
353.                                         %  3 and 40 are the top row of the matrix above.
354.                                                     %  All rows must have the same length
355.                                                     %  to form a matrix in the next step.
356.
357.
358. clear matcellsort   %  matcellsort is not needed anymore.
359.
360.
361.
362. toprows = cell2mat( toprowscell ) ;
363.
364. %                            whos toprows
365. %                              Name           Size            Bytes  Class     Attributes
366. %
367. %                              toprows      132x4              4224  double
368.
369.
370.
371. %%%%%  How many different rows in toprows?
372.
373. toprowsint = intersect(toprows,toprows,'rows') ;
374.
375. %                            toprowsint =
376. %                                 1    63     0     0
377. %                                 2     3    60     0
378. %                                 2     7    56     0
379. %                                 2    31     0     0
380. %                                 3     3    12    48
381. %                                 3     3    24    36
382. %                                 3     3    28     0
383. %                                 3     3    44     0
384. %                                 3     3    52     0
385. %                                 3     3    56     0
386. %                                 3     5    56     0
387. %                                 3    15     0     0
388. %                                 3    23     0     0
389. %                                 3    27     0     0
390. %                                 4     3    12     0
391. %                                 4     3    20     0
392. %                                 4     3    24     0
393. %                                 4     3    36     0
394. %                                 4     3    40     0
395. %                                 4     5    40     0
396. %                                 4     7     0     0
397. %                                 4    11     0     0
398. %                                 4    13     0     0
399. %                                 4    21     0     0
400. %                                 5     3     0     0
401. %                                 5     5     0     0
402. %                                 5     9     0     0
403. %                                 6     0     0     0
404.
405.
406.
407. %%%%%  How often does each number (of blocks) appear in the left column of toprowsint?
408.
409. row = zeros(1,num) ;
410. for m=1:num
411. row(m) = sum( toprowsint(:,1)==m ) ;
412. end
413.
414. row
415.
416. %                            row =
417. %                                 1     3    10    10     3     1
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top