% function P = get_rand_seq(n_seq,n_trials,reward_range,corr_tol,fun,include_fun)
%
% This function can be used to create multiple sequences that look like a
% random walk, but that
% (a) pass through the whole range rather than doing so by chance
% and 
% (b) all have the same frequency spectrum (i.e. the same rate of change). 
% The correlation range can be set to a specific value denoting how much
% these sequences are allowed to correlate.
% Please check the sequences, because just by chance specific unwanted
% patterns can emerge (such as random fluctuation with a specific shape).
%
% Please note: This function only checks for correlations across trials for
% all frequencies together, not for correlations across trials for a given
% frequency. This means it is possible that specific frequencies alone are
% highly correlated. For example, an observer may not perceive specific
% changes, but only very slow changes (equivalent to a low-pass filter).
% In this case, the sequences may be perceived as being correlated although
% they are not.
%
% Input:
%   n_seq: Number of sequences
%   n_trials: Number of trials for each sequence (I'd recommend creating
%       shorter blocks that all have the same sequence length, e.g. when
%       you have 200 trials, maybe create four sequence pairs of 50 trials)
%   reward_range: This is a 1x2 vector denoting the range of rewards we
%       will get later for our sequence (e.g. [0.2 0.8]). Please note that
%       we use Fourier transformation, i.e. that these ranges might be
%       exceeded. We control for that by slightly reducing the input range.
%       If you care about exactly reaching the extreme values, then
%       normalize your results (which will however influence the frequency
%       spectrum). (default: [0.05 0.95])
%   corr_tol: This is a 1x2 vector denoting the range in which the
%       sequences are allowed to be correlated (default: [-0.05 0.05]).
%   fun (optional): Per default, all values within reward_range are equally
%       likely (i.e. a linear function). If you want to change that, you
%       can manually enter the function yourself. In this case,
%       reward_range and n_trials will not play a role.
%       [default: linspace(reward_range(1),reward_range(2),n_trials) ]
%   include_fun (optional): 0 or 1. Whether to include the sequence in fun (default: 0)
%
% Example:
% P = get_rand_seq(2,64,[0.2 0.8])
% plot(P), axis([1 64 0 1])
% % check power spectrum difference
% pw1 = abs(fft(P(:,1))); pw2 = abs(fft(P(:,2)));
% % Should be close to zero
% range(pw1-pw2)
%
% 2014/03 Martin Hebart

function P = get_rand_seq(n_seq,n_trials,reward_range,corr_tol,fun,include_fun)

max_iter = 100000; % number of iterations for search (per sequence)

% set defaults
if ~exist('reward_range','var') || isempty(reward_range)
    reward_range = [0.05 0.95];
end

if ~exist('corr_tol','var') || isempty(corr_tol)
    corr_tol = [-0.05 0.05];
end

if ~exist('include_fun','var') || isempty(include_fun)
    include_fun = 0;
end

run_checks; % see function at end

% adjust reward range
m = mean(reward_range);
reward_range = 0.95*(reward_range-m)+m;
if ~exist('fun','var') || isempty(fun)
    x = linspace(reward_range(1),reward_range(2),n_trials)'; % default value
else
    x = fun;
end

% run fourier transformation, get power and phase
fftx = fft(x);
powerx = abs(fftx);
phasex = angle(fftx);

if include_fun % should we include the original function?
    P(:,1) = fun;
else
    % carry out phase scrambling and include scrambled into output
    randphase = angle(fft(rand(size(x)))); 
    % this below is the same as above, but shifting all phases by +/- pi/2 providing us with uncorrelated sequences within each frequency band
%     randphase(1) = phasex(1); randphase(2:end) = pi*round(rand(size(randphase(2:end))))-pi/2; 
    scrphase = phasex + randphase;
    P(:,1) = real(ifft(powerx.*exp(sqrt(-1)*(scrphase))));
    
end


for i_seq = 2:n_seq
    
    ct = 0; % counter for bailing out
    
    ind = tril(true(i_seq),-1); % index for comparison of correlations
    
    while 1
        
        ct = ct+1;
        
        % do same as above, keep power
        randphase = angle(fft(rand(size(x))));
%         randphase(1) = phasex(1); randphase(2:end) = pi*round(rand(size(randphase(2:end))))-pi/2;
        scrphase = phasex + randphase;
        p = real(ifft(powerx.*exp(sqrt(-1)*(scrphase))));
        
        % correlate all sequences
        r = corr([P p]);
        
        % check if all is in right range, otherwise repeat
        if all(r(ind)>corr_tol(1)) && all(r(ind)<corr_tol(2))
            P = [P p]; %#ok<AGROW> % update P
            break
        end
        
        % check if number of maximum iterations is exceeded
        if ct > max_iter
            fprintf('With given parameters couldn''t find %i sequences that match.\nPlease adjust parameters or edit function and change variable max_iter!\n',n_seq)
            P = [];
            return
        end
        
    end
    
end

if any(P(:)>1)
    warning('GET_RAND_SEQ:LARGER1','Some entries were larger than 1 and will be truncated to 0. This may slightly affect the frequency, mean, and std')
    P(P(:)>1) = 1;
end

if any(P(:)<0)
    warning('GET_RAND_SEQ:SMALLER0','Some entries were smaller than 0 and will be truncated to 0. This may slightly affect the frequency, mean, and std')
    P(P(:)<0) = 0;
end



function run_checks
    
    if n_seq < 2
        error('Minimum number of sequences is 2')
    end
    
    if any(reward_range<0|reward_range>1)
        error('Variable reward_range must be in range 0 to 1')
    end
    
    if reward_range(1) == reward_range(2)
        error('Range of variable reward_range must be larger than 0')
    end
    
    if reward_range(1)>reward_range(2)
        reward_range = reward_range([2 1]);
    end
    
    
    if corr_tol(1)>corr_tol(2)
        corr_tol = corr_tol([2 1]);
    end
    
    if exist('fun','var') && ~isempty(fun)
        if size(fun,1)==1
            fun = fun';
        end
        if ~isempty(n_trials) && n_trials ~= length(fun)
            warning('RUN_CHECKS:N_TRIALS_IGNORED','Input variable ''n_trials'' ignored. Using length of ''fun'' as number of trials.')
        end
        if reward_range(1) ~= min(fun) || reward_range(2) ~= max(fun)
            warning('RUN_CHECKS:REWARD_RANGE_IGNORED','Input variable ''reward_range'' ignored. Using range of ''fun'' as range.')
        end
    end
    
end

end