Upload files to 'code'

code/Map.m

@@ -0,0 +1,209 @@
+function [map,stats] = Map(v,z,varargin)
+
+%Map - Map z on (x,y) where x, y and z are time-varying variables (samples).
+%
+%  Compute a continuous map, where one time-varying variable z is represented
+%  as a function of one or two time-varying variables x and y. The variable z
+%  can either be a point process (typically, a list of spike timestamps) or a
+%  continuous measure (e.g. the instantaneous velocity of the animal, the
+%  spectral power of an LFP channel in a given frequency band, the coherence
+%  between two oscillating LFP channels, etc.) Typical examples of x and y
+%  include spatial coordinates and angular directions.
+%
+%  An occupancy map is also computed.
+%
+%  USAGE
+%
+%    map = Map([t x y],z,<options>)
+%
+%    t              timestamps for x and y
+%    x              x values in [0,1]
+%    y              optional y values in [0,1]
+%    z              list of timestamps or (timestamp,value) pairs
+%    <options>      optional list of property-value pairs (see table below)
+%
+%    =========================================================================
+%     Properties    Values
+%    -------------------------------------------------------------------------
+%     'smooth'      smoothing size in bins (0 = no smoothing, default = 2)
+%     'nBins'       number of horizontal and vertical bins (default = [50 50])
+%     'minTime'     minimum time spent in each bin (in s, default = 0)
+%     'maxGap'      z values recorded during time gaps between successive (x,y)
+%                   samples exceeding this threshold (e.g. undetects) will not
+%                   be interpolated; also, such long gaps in (x,y) sampling
+%                   will be clipped to 'maxGap' to compute the occupancy map
+%                   (default = 0.100 s)
+%     'type'        'linear' if z is linear (default), 'circular' otherwise
+%    =========================================================================
+%
+%  OUTPUT
+%
+%    map.x          x bins
+%    map.y          y bins
+%    map.z          average map (z continuous)
+%                   or rate map (z point process)
+%    map.count      count map (z point process)
+%    map.time       occupancy map (in s)
+%
+%  NOTES
+%
+%    x values are arranged in columns and y values in rows in all output matrices
+%    (e.g. 'map.z').
+%
+%  SEE
+%
+%    See also MapStats, FiringMap, PlotColorMap, Accumulate.
+
+% Copyright (C) 2002-2011 by Michael Zugaro
+%
+% This program is free software; you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation; either version 3 of the License, or
+% (at your option) any later version.
+
+% Check number of parameters
+if nargin < 2 | mod(length(varargin),2) ~= 0,
+  error('Incorrect number of parameters');
+end
+
+% Check parameter sizes
+if size(v,2) < 2,
+	error('Parameter ''[t x y]'' should have at least 2 columns');
+end
+if (size(z,2) < 1 || size(z,2) > 2) && ~isempty(z),
+	error('Parameter ''z'' should have 1 or 2 columns');
+end
+
+% Default values
+maxGap = 0.1;
+map.x = [];
+map.y = [];
+map.count = [];
+map.time = [];
+map.z = [];
+smooth = 1;
+nBins = [100 1];
+minTime = 0;
+type = 'linear';
+
+% 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 'smooth',
+			smooth = varargin{i+1};
+			if ~isdvector(smooth) | length(smooth) > 2,
+				error('Incorrect value for property ''smooth''');
+			end
+
+		case 'nbins',
+			nBins = varargin{i+1};
+			if ~isivector(nBins) | length(nBins) > 2,
+				error('Incorrect value for property ''nBins''');
+			end
+
+		case 'mintime',
+			minTime = varargin{i+1};
+			if ~isdscalar(minTime,'>=0'),
+				error('Incorrect value for property ''minTime''');
+			end
+
+		case 'maxgap',
+			maxGap = varargin{i+1};
+			if ~isdscalar(maxGap,'>=0'),
+			error('Incorrect value for property ''maxGap''');
+			end
+
+		case 'type',
+			type = lower(varargin{i+1});
+			if ~isstring(type,'circular','linear'),
+				error('Incorrect value for property ''type''');
+			end
+
+		otherwise,
+			error(['Unknown property ''' num2str(varargin{i}) '''']);
+
+  end
+end
+
+if isempty(v), return; end
+
+% Some info about x, y and z
+pointProcess = (isempty(z) | size(z,2) == 1);
+t = v(:,1);
+x = v(:,2);
+if size(v,2) >= 3,
+	y = v(:,3);
+else
+	y = [];
+end
+
+% Number of bins for x and y
+nBinsX = nBins(1);
+if length(nBins) == 1,
+	nBinsY = nBinsX;
+	nBins(2) = nBins;
+else
+	nBinsY = nBins(2);
+end
+
+% Bin x and y
+x = Bin(x,[0 1],nBinsX);
+if ~isempty(y),
+	y = Bin(y,[0 1],nBinsY);
+end
+
+% Duration for each (X,Y) sample (clipped to maxGap)
+dt = diff(t);dt(end+1)=dt(end);dt(dt>maxGap) = maxGap;
+
+if pointProcess,
+	% Count occurrences for each (x,y) timestamp
+	n = CountInIntervals(z,[t t+dt]);
+%  	n = CountInIntervals(z,[t(1:end-1) t(2:end)]);
+%  	n(end+1) = 0;
+else
+	% Interpolate z at (x,y) timestamps
+	[z,discarded] = Interpolate(z,t,'maxGap',maxGap);
+	if isempty(z), return; end
+	if strcmp(type,'circular'),
+		range = isradians(z(:,2));
+		z(:,2) = exp(j*z(:,2));
+	end
+	n = 1;
+end
+
+% Computations
+if isempty(y),
+	% 1D (only x)
+	map.x = linspace(0,1,nBinsX);
+	map.count = Smooth(Accumulate(x,n,nBinsX),smooth)';
+	map.time = Smooth(Accumulate(x,dt,nBinsX),smooth)';
+	if pointProcess,
+		map.z = map.count./(map.time+eps);
+	else
+		map.z = Smooth(Accumulate(x(~discarded),z(:,2),nBinsX),smooth)';
+		map.z = map.z./(map.count+eps);
+	end
+else
+	% 2D (x and y)
+	map.x = linspace(0,1,nBinsX);
+	map.y = linspace(0,1,nBinsY);
+	map.count = Smooth(Accumulate([x y],n,nBins),smooth)';
+	map.time = Smooth(Accumulate([x y],dt,nBins),smooth)';
+	if pointProcess,
+		map.z = map.count./(map.time+eps);
+	else
+		map.z = Smooth(Accumulate([x(~discarded) y(~discarded)],z(:,2),nBins),smooth).';
+		map.z = map.z./(map.count+eps);
+	end
+end
+
+% Circular z
+if strcmp(type,'circular'), map.z = wrap(angle(map.z),range); end
+
+% Discard regions with insufficient sampling
+map.z(map.time<=minTime) = 0;
+

code/PosMapOrdered.m

@@ -0,0 +1,140 @@
+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

code/PositionDecoding.m

@@ -0,0 +1,148 @@
+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);
+
+