function All_seqs = BRM_list_rand_algo(no_As, no_Bs, no_Cs, no_Ds, no_of_seqs)

% You make a shuffled list of As, Bs, Cs, and Ds, with no_As, no_Bs, no_Cs 
% and no_Ds number of each item.  
% Algorithm: Remove the repeated element of all doubles to obtain S.  The 
% repeats constitute a list, R. The elements of R are then inserted into 
% RANDOM locations in S.  
% Each time you insert a random elt from R into S, you then reshuffle S.
% You then randomize the numbers of the insertion locations in S.
% For each element in S, you run down this rand index list to find a place
% to insert the elt from R. 
% example:
% S = BRM_list_rand_algo(0, 6, 12, 12, 10000);
% will generate 10,000 sequences of 30 elements, each containing 0 As, 6 Bs, 
% 12 Cs and 12 Ds.
% Generalizing this code to more item types is straightforward.

global MAX_CYCLES
MAX_CYCLES = 10000;

tic
All_seqs = [];

% We create a list with the number of 1s, 2s, 3s, and 4s indicated in the 
% parameter list (Note: the labels 1,2,3,4 replace A, B, C, and D).  
Orig_W = [1*ones(1,no_As), 2*ones(1,no_Bs), 3*ones(1,no_Cs), 4*ones(1, no_Ds)];
len_W = length(Orig_W);

seq_no = 1;
while seq_no <= no_of_seqs
  % first we remove repeats.  Start with the no-repeat list W and the list
  % of repeats, R.
  W = Orig_W;
  
% We shuffle W.
  W(randperm(len_W)) = W;
  W_shift = [W(end), W(1:end-1)];
  indices_items_to_remove = find([eq(W, W_shift)] == 1);
  R = W(indices_items_to_remove);
  W(indices_items_to_remove) = [];
  S  = W;
% S now has no repeats and its first and last elts are different.
% R contains the repeated elements to be inserted

% We "shuffle" S -- actually, cut S like a deck of cards, as we will do 
% throughout the algorithm.  Since the first and last elts are different, 
% shuffling in this manner produces no new repeats
  cut_point = mod(seq_no, length(S));
  S =[S(1+cut_point:end), S(1:cut_point)];
  
% a counter allows us to get out of the while loop if we loop for too long
  ctr = 0;
  
  while ~isempty(R)
    
% if we cycle too long, break, and try again
    ctr = ctr+1;
    if ctr > MAX_CYCLES
      break;
    end;
    
% pick a random element from R
    item_index = ceil(length(R)*rand);
    cur_item = R(item_index);

% shuffle S (exactly like cutting a deck of cards)
    S_len = length(S);
    cut_point = ceil(length(S)*rand);
    S =[S(1+mod(cut_point, S_len):end), S(1:cut_point)];

% randomize the list of indices of S    
    S_indices = randperm(S_len);

% go down this list of indices until one of them is found for which the
% element from R can be inserted, i.e., the neighbors on either side of the
% elt. from R are not the same as the item from R.  Insert the elt. from R
% at this point.
% If the insertion point is the last item in S, then we add the elt. from R
% either to the beginning or to the end of S.
    for rand_index = S_indices
      prev_item = S(rand_index);
      succ_item = S(mod(rand_index, S_len) + 1);
% the mod(rand_index, S_len) is necessary if the rand_index is the index of
% the last item in S.  This wraps everything around.  Thus the successor of
% the last item in S is the first item in S.

% check that the immediate neighors of the cur_item (i.e., prev_item and 
% succ_item) taken from R are not the same as the cur_item.  If they are
% not the same and the index is not the index of the final elt in S, insert
% the cur_item between prev_item and succ_item in S.  If the insertion
% point is the final item of the list, half the time insert after the end
% of S, half the time before the beginning of S.
      if ne(cur_item, prev_item) && ne(cur_item, succ_item)
         if rand_index ~= S_len
          S = [S(1:rand_index), cur_item, S(rand_index+1:end)];
        else
          if rand > 0.5
            S = [S, cur_item];
          else
            S = [cur_item, S];
          end;
        end;
% when the insertion of the item from R is completed, remove it from R and
% leave the loop in order to pick a new item from R to insert
        R(item_index) = [];
        break;
       end;
    end;
  end;
 
% Build up a matrix of all of the sequences found
  if length(S) == len_W
    All_seqs = [All_seqs; S];
    seq_no = seq_no + 1;
  end; 

end;
toc
return;


