RSC_OpticFlow/code/PosMapOrdered.m

function [ idx, pos_map_ordered, pos_map, pos_map_trial, pos_map_trial_norm, pos_map_norm_min, pos_map_norm_max]...
    = PosMap_ordered( tcs, deconv, behavior, varargin )

% position firing maps are calculated for the entire session
% then reshape to position bins x trials
% check PlotPosMapTcs.m

% default
nBins = 100;
nSmooth = 3;
n_lowpass = 0.3;
n_medf = 10;
use_deconv = false;
use_lowpass = true;
unit_list = 1:size(deconv,2);

% parse input parameters
% Parse parameter list
for i = 1:2:length(varargin),
    if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property']);
    end
    switch(lower(varargin{i}))
        case 'nsmooth'
            nSmooth = varargin{i+1};
            if ~isscalar(nSmooth) || length(nSmooth) > 2,
                error('Incorrect value for property ''nSmooth''');
            end
		case 'unit_list',
            unit_list = varargin{i+1};
            if ~isvector(unit_list),
                error('Incorrect value for property ''unit_list''');
            end
        case 'use_deconv',
			use_deconv = varargin{i+1};
            if ~islogical(use_deconv) || length(use_deconv) > 2,
				error('Incorrect value for property ''use_deconv''');
            end
        case 'use_lowpass'
            use_lowpass = varargin{i+1};
            if ~islogical(use_lowpass)
                error('Incorrect value for property ''use_lowpass''');
            end

		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) '''']);
    end
end

nUnit = length(unit_list);
trial = behavior.trial;
nTrial = max(trial);
nBins_alltrial = [nBins*nTrial 1];
ts_beh = behavior.ts;
pos_cum = behavior.pos_cum;
ok = trial > 0;
ts_beh = ts_beh(ok);
pos_cum = pos_cum(ok);
pos_cum = (pos_cum-pos_cum(1))/(pos_cum(end)-pos_cum(1));
if isrow(ts_beh)
    ts_beh = ts_beh';
end
if isrow(pos_cum)
    pos_cum = pos_cum';
end

ts_pos_beh = [ts_beh, pos_cum./max(pos_cum)];

ts_tcs = (tcs.tt)';
ok = ts_tcs >= ts_beh(1) & ts_tcs <= ts_beh(end);
ts_tcs = ts_tcs(ok);
if isrow(ts_tcs)
    ts_tcs = ts_tcs';
end

if use_lowpass
    ratio = lowpass(tcs.ratio, n_lowpass);
else
    ratio = medfilt1(double(tcs.ratio), n_medf);
end

pos_map = zeros(nUnit, nBins(1));
pos_map_norm_min = zeros(nUnit, 1);
pos_map_norm_max = zeros(nUnit, 1);
pos_map_ordered = zeros(nUnit, nBins(1));

pos_map_trial = zeros(nTrial, nBins, nUnit);
pos_map_trial_norm = zeros(nTrial, nBins, nUnit);

unit_order = [];
if ~use_deconv
    ratio = ratio(ok, :);
    y = ratio;
else
    deconv = deconv(ok, :);
    y = deconv;
end

for iunit = 1:nUnit

    id_unit = unit_list(iunit);
    y_this = y(:, id_unit);
        
    map = Map(ts_pos_beh, [ts_tcs, y_this], 'smooth', nSmooth, 'nBins', nBins_alltrial); 
    pos_map_tmp = reshape(map.z, nBins, nTrial);

    pos_map_tmp = pos_map_tmp';

    pos_map_trial(:,:,iunit) = pos_map_tmp;

    [max_val, ~] = max(pos_map_tmp, [], 2);
    [min_val, ~] = min(pos_map_tmp, [], 2);

    norm_val = max_val - min_val;
    ok = norm_val > 0.1;   % silent trials, normalized by 1 (i.e. 0/1 = 0)
    pos_map_tmp(~ok, :) = 0;

    norm_val(~ok) = 1;
    pos_map_tmp = bsxfun(@minus, pos_map_tmp, min_val);
    pos_map_trial_norm(:,:,iunit) = bsxfun(@rdivide, pos_map_tmp, norm_val);

    a = mean(pos_map_tmp,1);
    pos_map(iunit, :) = a;
    pos_map(iunit, :) = (a-min(a))/(max(a)-min(a));
    pos_map_norm_min(iunit, :) = min(a);
    pos_map_norm_max(iunit, :) = max(a);

    [~, idx] = max(pos_map(iunit, :));
    unit_order = cat(1, unit_order, [id_unit, idx]);
    
end
    
[~, idx] = sort(unit_order(:, 2));

for iunit = 1:nUnit
    i = idx(iunit);
    pos_map_ordered(iunit, :) = pos_map(i, :);
end

end