% rna_loops generates an RNA loop of length about n with about L additional loops inserted at random places

function [Pi] = rna_loops(n,L)

Pair = ['AU'; 'CG'; 'AA'; 'CC'; 'GG'; 'UU'; 'AC'; 'AG'; 'CU'; 'GU'; 'AX'; 'CX'; 'GX'; 'UX']; 
                                     % all possible base pairs
                                     % X will be deleted at the end
Freq = [40 40 2 2 2 2 2 2 2 2 1 1 1 1];    % base pair relative frequencies
Freq = Freq/sum(Freq);               % make sure Freq sums to 1

Loop = ['AAAG'; 'GGAG'; 'CGGU'; 'CUAC'];   % possible loop ends
L    = length(Loop(:,1));            % number of possible loops

Pi   = Loop(rando(ones(L,1)/L),:);   % start with a random loop end

for i=1:round(n/2),
  b = rando(Freq);
  s  = rando([1 1]/2);               % code to insert pair as AU or UA
  Pi = [Pair(b,s) Pi Pair(b,3-s)];   % add this pair to front and back of Pi
  if (rand < 2*L/n)                  % on average, add L loops
    if (rand < 0.5)
      Pi = [canonical_loop(rando(ones(n,1)/n)) Pi]; % prepend another loop
    else
      Pi = [Pi canonical_loop(rando(ones(n,1)/n))]; % append another loop
    end
  end

end

Pi = Pi(find(Pi ~= 'X'));            % remove the X's