daily pastebin goal
64%
SHARE
TWEET

Noncrossing partitions up to rotation

mate2code Feb 13th, 2013 163 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. % See also: http://oeis.org/A209805
  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. OK, I Understand
 
Top