How to fill up a matrix with ones from their row and column distributions

Good evening,

First of all, may you have a happy new year 2017!!

I've started the year dealing with a couple of problems from a project I need to hand in. I may end up probably asking about both of them, but by now I would like to ask you just about the first one.

I will try to explain myself as well as I can:

I need to fill up a matrix with '1' and '0'. Data provided are matrix' dimensions and '1's rows and columns distributions (through two different vectors, as they can be different for each case). So we have two vectors, v and h, whose elements are defined as:

v_i = i-th element of vector v. Indicates the fraction of columns of weight i (fraction of columns over total which contain i '1's).fh_i = i-th element of vector h. Indicates the fraction of rows of weight i (fraction of rows over total which contain i '1's).

The procedure taken into account is meant as MacKay Neal algorithm. Consists of filling up the matrix in a column basis, from left to right. The weight of each column is chosen to obtain a valid distribution of '1's. The location of these '1' in each file is chosen randomly between those which are not yet full.

Please find attached my code right below. You may also find in this webpage the PDF document I am consulting right now to perform this part of my project. Everything I mentioned beforehand can be found in pages 10 and 11. MacKay algorithm is implemented in pseudocode in page 14.

    function H = MacKayNeal(N,r,v,h)
    H = zeros(N*(1-r),N);
    a = [];
    for i=1:length(v)
        for j=1:(v(i)*N)
            a = [a,i];
        end
    end
    b = [];
    for i=1:length(h)
        for j=1:(h(i)*(N*(1-r)))
            b = [b,i];
        end
    end
    for i=1:N
       c = datasample(b,length(a(i)),'Replace',false);
       for j=1:a(i)
           H(c(j),i)=1;
       end
    a = a - c;
    end
%     while 1 % cycles removed
%         for i=1:(N-1)
%             for j=(i+1):N
%                 if 1 %
%                 end
%             end
%         end
%     end              
    end

Note that values 'a' and 'b' are respectively referring to 'alpha' and 'beta' variables used in attached text's pseudocode. Their values indicate the number of '1' per column and file in each case.

In my case everything went alright until I performed substraction ' a = a - c' in line 20. I spent a whole afternoon racking my brains with it, and two questions ended up coming out:

- Why is this substraction 'a = a - c' and not 'b = b - [one unity in each corresponding position from c]? It would make much more sense for me with respect to what does appear in theory explanation prior to pseudocode in the PDF attached in the link.

- What would happen if you choose a subset of b which contains at least two equal values to perform the following assignation of '1's in H matrix? To me it seems as a perfectly feasible situation in which a_i '1's would not be added, but less of that.

If anyone may also help me out with solving the last part of the task, the one alluding in pseudocode the erasure of length 4' cycles (that simply means that there cannot be in two separate columns two pairs of '1's occupying the very same positions) I would also be really really thankful.

In my case I am working with the following data:

    N = 100
    k = 0.5
    v = [0 0.17 0.03 0.2 0 0 0.17 0.3 0.03 0.1]
    h = [0 0.1 0 0 0 0.7 0.2]

For both vectors v and h, their elements can take any value from 0 to 1 as long as their sum equals 1 in both cases, and first element equals 0 also in both cases. Also, k tan take any value between 0 and 1.

Many thanks in advance for your attention!