function [perm_matrix,perm_orders] = fmri_gen_perms(initial_order,num_perms);

%
% USAGE: [perm_matrix,perm_orders] = fmri_gen_perms(initial_order,num_perms);
%
%   Generate the order matrix for permutation test  
%
%   INPUT:
%	initial_order - Nx1 vector, where N is the number of event onsets,
%			specifies the event type for each onset.  The
%			elements in the vector are the indices of the
%                       event types.  The values are in the range of 
%			[ 1..T ].  1 for Evt#1, ..., n for Evt#T.
%       num_perms - number of permutations will be performed 
%
%   OUTPUT:
%	perm_matrix - NxM matrix, where N is the number of event onsets, 
%                     and M is equal to num_perms.  Each column contains 
%		      the indices of event types for the onsets. 
%	perm_orders - Matrix with the same size as perm_matrix.  It stores
%		      the values of permuation order with respect to the
%                     initial_order. That is, for the ith permutation, 
%                        perm_matrix(:,i) = initial_order(perm_orders(:,i))
%
%   Example:
%
%      initial_order = [ 1 1 1 2 2 2 3 3 3 4 4 4 ];
%      [pm,po] = fmri_gen_perms(initial_order,50); 
%
%
%   -- Created Nov 2000 by Wilkin Chau, Rotman Research Institute
%		

  initial_order = initial_order(:);    % make sure it is column vector
  normalized_order = initial_order / sqrt(initial_order'  * initial_order);
  
  num_onsets = length(initial_order);
  perm_orders = zeros(num_onsets,num_perms);
  perm_matrix = zeros(num_onsets,num_perms);
  
  for i=1:num_perms,
     new_order_found = 0;
     while (~new_order_found)
        new_order = randperm(num_onsets);

        %  make sure the new permutation order is not a repeat one
        %  
	d = abs(perm_matrix(:,1:i-1) - initial_order(new_order)*ones(1,i-1)); 
        matched_idx = find (sum(d) == 0);
        if isempty(matched_idx) 
    	   new_order_found = 1;
        end

     end;
  
     perm_orders(:,i) = new_order';
     perm_matrix(:,i) = initial_order(new_order);
  end;
  
  return;