Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- % See also: http://oeis.org/A209805
- % This MATLAB code can be used to calculate the rows of the number triangle OEIS A209805,
- % i.e. the number of k-block noncrossing partitions of an n-element set up to rotation.
- % This is a brute force method and may cause memory problems for n > 12.
- % All objects are created and than counted,
- % so this code can also be used to find the actual partitions.
- % Written by Tilman Piesk, March 2012.
- %
- %
- % The functions "SetPartition" and "Bell" by Bruno Luong from Matlab Central are used,
- % as well as the following two functions:
- function [output] = crossing(A,B)
- % Two sets crossing?
- % Enter sets as row vectors.
- cross = 0;
- if min(B) < min(A)
- H = A ; A = B ; B = H ;
- end
- x2 = min(B) ;
- if max(A) > x2
- x3 = min( setdiff(A,1:x2) ) ;
- if max(B) > x3
- cross = 1 ;
- end
- end
- output = cross ;
- end
- function [y] = perm2mat(x)
- % Input: Permutation of elements 1 ... n as row vector
- % Output: Permutation matrix (to permute rows, e.g. elements of a column vector)
- wide = size(x,2) ;
- I = eye(wide) ;
- for k = 1:wide
- P( k , 1:wide ) = I( x(k) , : ) ;
- end
- y = P ;
- end
- % The row number is called "num" and has to be entered before calculation.
- % To explain the code some outputs for num = 6 have been included as comments.
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % num =
- % 6
- bell = Bell(num) ;
- % bell =
- % 203
- powtwo = zeros(num,1) ;
- for m=1:num
- powtwo(m) = 2^(m-1) ;
- end
- % powtwo =
- % 1
- % 2
- % 4
- % 8
- % 16
- % 32
- %%%%% All partitions:
- part1 = SetPartition(num) ;
- % part1{20}
- % ans =
- % [1x3 double] [1x3 double]
- % part1{20}{1}
- % ans =
- % 1 2 4
- % part1{20}{2}
- % ans =
- % 3 5 6
- % part1{50}
- % ans =
- % [1x2 double] [3] [1x2 double] [5]
- % part1{50}{1}
- % ans =
- % 1 2
- % part1{50}{2}
- % ans =
- % 3
- % part1{50}{3}
- % ans =
- % 4 6
- % part1{50}{4}
- % ans =
- % 5
- %%%%% How many blocks does each partition have?
- part1long = zeros(bell,1) ;
- for m=1:bell
- part1long(m) = size(part1{m},2) ;
- end
- % part1long(20)
- % ans =
- % 2
- % part1long(50)
- % ans =
- % 4
- %%%%% How many non-singleton blocks does each partition have?
- partlong = zeros(bell,1) ;
- for m=1:bell
- M=0;
- for n=1:part1long(m)
- if size(part1{m}{n},2) > 1
- M = M+1 ;
- end
- end
- partlong(m) = M ;
- end
- % partlong(20)
- % ans =
- % 2
- % partlong(50)
- % ans =
- % 2
- %%%%% All partitions, singleton blocks are not shown:
- part = cell(bell,1) ;
- for m=1:bell
- i = 1 ;
- for n=1:part1long(m)
- N = part1{m}{n} ;
- if size(N,2) > 1
- part{m}{i} = N ;
- i = i+1 ;
- end
- end
- end
- % part{50}
- % ans =
- % [1x2 double] [1x2 double]
- % part{50}{1}
- % ans =
- % 1 2
- % part{50}{2}
- % ans =
- % 4 6
- clear part1 % part1 is not needed anymore.
- %%%%% Which partitions are crossing?
- cross = zeros(bell,1) ;
- for a=1:bell
- c = 0 ;
- if partlong(a) > 1
- A = nchoosek(1:partlong(a),2); % index numbers to get all pairs of blocks
- for b=1:size(A,1)
- B = A(b,:) ;
- B1 = part{a}{B(1)};
- B2 = part{a}{B(2)};
- c = crossing(B1,B2) ; % The function crossing tells whether
- if c == 1 % two blocks are crossing (1) or not (0).
- break
- end
- end
- end
- cross(a) = c ;
- end
- % cross(20)
- % ans =
- % 1
- % cross(50)
- % ans =
- % 0
- %%%%% The blocks of the partitions shall be represented by binary vectors:
- %%%%% E.g. the block {1,2,6} by the vector 110001.
- matcell = cell(bell,1) ;
- for m=1:bell
- if cross(m) == 0 % If the partition is crossing matcell(m) is left empty.
- high = partlong(m) ;
- M = zeros( high , num ) ;
- for n=1:high
- N = part{m}{n} ;
- for p=1:size(N,2)
- M( n , N(p) ) = 1 ;
- end
- end
- matcell{m} = M ;
- end
- end
- % matcell{20}
- % ans =
- % []
- % matcell{50}
- % ans =
- % 1 1 0 0 0 0
- % 0 0 0 1 0 1
- clear part % part is not needed anymore.
- %%%%% Turning the partitions by permuting the binary vectors:
- %%%%% 110001 into 100011, 000111, 001110 etc.
- %%%%% First generate permutation matrices:
- permcell = cell(num,1) ;
- perm = perm2mat([num 1:(num-1) ]) ;
- for m=1:num
- M = perm^(m-1) ;
- permcell{m} = M ;
- end
- % permcell{2}
- % ans =
- % 0 0 0 0 0 1
- % 1 0 0 0 0 0
- % 0 1 0 0 0 0
- % 0 0 1 0 0 0
- % 0 0 0 1 0 0
- % 0 0 0 0 1 0
- % permcell{3}
- % ans =
- % 0 0 0 0 1 0
- % 0 0 0 0 0 1
- % 1 0 0 0 0 0
- % 0 1 0 0 0 0
- % 0 0 1 0 0 0
- % 0 0 0 1 0 0
- %%%%% Permute binary vectors and turn results into decimal numbers,
- %%%%% e.g. vector 110000 into number 3.
- matcellbig = cell(bell,num) ;
- for m=1:bell
- if cross(m) == 0
- for n=1:num
- matcellbig{m,n} = matcell{m} * permcell{n} * powtwo ;
- end
- end
- end
- % matcellbig{50,1}
- % ans =
- % 3
- % 40
- % matcellbig{50,2}
- % ans =
- % 33
- % 20
- % matcellbig{50,3}
- % ans =
- % 48
- % 10
- % matcellbig{50,4}
- % ans =
- % 24
- % 5
- % matcellbig{50,5}
- % ans =
- % 12
- % 34
- % matcellbig{50,6}
- % ans =
- % 6
- % 17
- clear matcell % matcell is not needed anymore.
- %%%%% These numbers have to be sorted and put into a matrix,
- %%%%% so it can be checked whether two partitions generate the same matrix:
- matcellsort=cell(bell,1);
- for m=1:bell
- if cross(m) == 0
- N = zeros( partlong(m) , num ) ;
- for n=1:num
- N(1:partlong(m),n) = sort(matcellbig{m,n}) ;
- end
- matcellsort{m} = sortrows(N') ;
- end
- end
- % matcellsort{50}
- % ans =
- % 3 40
- % 5 24
- % 6 17
- % 10 48
- % 12 34
- % 20 33
- clear matcellbig % matcellbig is not needed anymore.
- %%%%% The top row can represent the whole matrix:
- toprowscell = cell(bell,1) ;
- for m=1:bell
- if cross(m) == 0
- M = zeros( 1 , max(partlong) ) ;
- M( 1:partlong(m) ) = matcellsort{m}(1,:) ;
- toprowscell{m} = [ part1long(m) M ] ;
- end
- end
- % toprowscell{50}
- % ans =
- % 4 3 40 0
- % The 4 is the number of blocks in the 50-th partition.
- % 3 and 40 are the top row of the matrix above.
- % All rows must have the same length
- % to form a matrix in the next step.
- clear matcellsort % matcellsort is not needed anymore.
- toprows = cell2mat( toprowscell ) ;
- % whos toprows
- % Name Size Bytes Class Attributes
- %
- % toprows 132x4 4224 double
- %%%%% How many different rows in toprows?
- toprowsint = intersect(toprows,toprows,'rows') ;
- % toprowsint =
- % 1 63 0 0
- % 2 3 60 0
- % 2 7 56 0
- % 2 31 0 0
- % 3 3 12 48
- % 3 3 24 36
- % 3 3 28 0
- % 3 3 44 0
- % 3 3 52 0
- % 3 3 56 0
- % 3 5 56 0
- % 3 15 0 0
- % 3 23 0 0
- % 3 27 0 0
- % 4 3 12 0
- % 4 3 20 0
- % 4 3 24 0
- % 4 3 36 0
- % 4 3 40 0
- % 4 5 40 0
- % 4 7 0 0
- % 4 11 0 0
- % 4 13 0 0
- % 4 21 0 0
- % 5 3 0 0
- % 5 5 0 0
- % 5 9 0 0
- % 6 0 0 0
- %%%%% How often does each number (of blocks) appear in the left column of toprowsint?
- row = zeros(1,num) ;
- for m=1:num
- row(m) = sum( toprowsint(:,1)==m ) ;
- end
- row
- % row =
- % 1 3 10 10 3 1
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
