RSC_OpticFlow/code/PositionDecoding.m

function [pos_err, pos_real] = PositionDecoding(pos_map_trial, tcs, deconv, behavior, varargin);

% input::
% pos_map_trial:
% occupancy-normalized position activity map in the format of nTrial x nBin x nUnit

% output::
% pos_err:
% absolute difference between decoded and actual position in cm
% pos_real:
% actual position in cm

order_idx = 1:size(tcs.ratio, 2);
unit_list = 1:size(tcs.ratio, 2);
nTrial = max(behavior.trial);
spd_thrsd = 3;     % speed threshold: 3 cm/s
tau = 0.5;         % time window
tracklen = 90;     % in cm
ifplot = false;

% 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 'tracklen'
            tracklen = varargin{i+1};
        case 'tau'
            tau = varargin{i+1};
            if ~isvector(tau),
                error('Incorrect value for property ''unit_list''');
            end
        case 'unit_list'
            unit_list = varargin{i+1};
            if ~isvector(unit_list),
                error('Incorrect value for property ''unit_list''');
            end
        case 'order_idx'
            order_idx = varargin{i+1};
            if ~isvector(order_idx),
                error('Incorrect value for property ''order_idx''');
            end
        case 'ifplot'
            ifplot = varargin{i+1};
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) '''']);
    end
end

ts_tcs = tcs.tt;
ts_beh = behavior.ts;
pos = behavior.pos_norm * tracklen;  % scale to track length
trial = behavior.trial;
spd = behavior.speed;

nX = round(tau ./ tcs.tt(2));
nNeuron = length(unit_list);
deconv = deconv(:, unit_list);

deconv_sm_size = 10;
deconv_sm = Smooth(double(deconv), [deconv_sm_size, 0]);

deconv_ordered = deconv_sm(:, order_idx);
% deconv_ordered(deconv_ordered < 0.3) = 0;

pos_tuning_curve = squeeze(mean(pos_map_trial(1:2:end,:,order_idx),1)); % use only odd trials for training
mltp_factor = exp(-tau * sum(pos_tuning_curve, 2));
nBin = length(mltp_factor);

pos_decoded = [];
ts_decoded = [];
trial_decoded = [];

for iTrial = 1:1:nTrial
    
    ok = trial == iTrial;
    ts_this = ts_beh(ok);
    ok2 = ts_tcs > ts_this(1) & ts_tcs < ts_this(end);
    len_this = sum(ok2);
    if len_this > 0
        dcv_this = deconv_ordered(ok2, :);
        ts_this = ts_tcs(ok2);
        for iX = 1:ceil(len_this/nX)
            st = (iX - 1) * nX + 1;
            ed = iX * nX;
            if ed > len_this
                ed = len_this;
            end
            ts_decoded = cat(1, ts_decoded, mean(ts_this(st:ed)));
            dcv_tmp = mean(dcv_this(st:ed, :), 1);
            dcv_mat_tmp = repmat(dcv_tmp, nBin, 1);
            prob_pos = prod(pos_tuning_curve.^dcv_mat_tmp, 2) .* mltp_factor;
            [~, idx] = max(prob_pos);
            pos_decoded = cat(1, pos_decoded, idx);  % need to be normalized to track length
            trial_decoded = cat(1, trial_decoded, iTrial * ones(length(idx),1));
        end
    end
    
end

spd_decoded = interp1(ts_beh, spd, ts_decoded, 'linear');
ok3 = spd_decoded > spd_thrsd;
st_run = find(diff([0; ok3]) == 1);  % start of a running period
ed_run = find(diff([0; ok3]) == -1);

if length(st_run) > length(ed_run)
    st_run = st_run(1:end-1);
end

if ifplot

    hold on
    for iTrial = 1:nTrial
        ok = trial==iTrial;
        plot(ts_beh(ok), pos(ok)/tracklen*nNeuron, 'b');
    end
    for iRun = 1:length(st_run)
        idx_this = st_run(iRun):ed_run(iRun);
        ts_decoded_this = ts_decoded(idx_this);
        pos_decoded_this = pos_decoded(idx_this)/nBin*nNeuron;
        pos_break = find(abs(diff(pos_decoded_this))>round(nNeuron*0.8));
        if isempty(pos_break)
            plot(ts_decoded_this, pos_decoded_this, 'r')
        else
            pos_break = [0; pos_break; length(pos_decoded_this)];
            for i = 1:length(pos_break)-1
                idx_now = pos_break(i)+1:pos_break(i+1);
                plot(ts_decoded_this(idx_now), pos_decoded_this(idx_now), 'r')
            end          
        end
    end
    hold off
end

ts_decoded = ts_decoded(ok3);
pos_decoded = pos_decoded(ok3) / nBin * tracklen;
trial_decoded = trial_decoded(ok3);

pos_true = interp1(ts_beh, pos, ts_decoded, 'linear');

ok4      = rem(trial_decoded, 2) == 0;   % calculate errors for only the even trials
pos_err  = abs(pos_decoded(ok4) - pos_true(ok4));
pos_real = pos_true(ok4);
trial_decoded = trial_decoded(ok4);