Map

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([t1 x y],[t2 z],<options>)

    t1             timestamps for x and y
    x              x values in [0,1]
    y              optional y values in [0,1]
    t2             timestamps for z
    z              optional z values
    <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)
     'mode'        'interpolate' to interpolate missing points (< minTime)
                   (default for a continuous z) or 'discard' to discard them
                   (default for a point process)
     'maxDistance' maximal distance for interpolation (default = 5)
     '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'        three letters (one for X, one for Y and one for Z) indi-
                   cating which coordinates are linear ('l') and which are
                   circular ('c') - for 1D data, only two letters are used
                   (default 'lll')
    =========================================================================

  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.

0001 function [map,stats] = Map(v,z,varargin)
0002 
0003 %Map - Map z on (x,y) where x, y and z are time-varying variables (samples).
0004 %
0005 %  Compute a continuous map, where one time-varying variable z is represented
0006 %  as a function of one or two time-varying variables x and y. The variable z
0007 %  can either be a point process (typically, a list of spike timestamps) or a
0008 %  continuous measure (e.g. the instantaneous velocity of the animal, the
0009 %  spectral power of an LFP channel in a given frequency band, the coherence
0010 %  between two oscillating LFP channels, etc.) Typical examples of x and y
0011 %  include spatial coordinates and angular directions.
0012 %
0013 %  An occupancy map is also computed.
0014 %
0015 %  USAGE
0016 %
0017 %    map = Map([t1 x y],[t2 z],<options>)
0018 %
0019 %    t1             timestamps for x and y
0020 %    x              x values in [0,1]
0021 %    y              optional y values in [0,1]
0022 %    t2             timestamps for z
0023 %    z              optional z values
0024 %    <options>      optional list of property-value pairs (see table below)
0025 %
0026 %    =========================================================================
0027 %     Properties    Values
0028 %    -------------------------------------------------------------------------
0029 %     'smooth'      smoothing size in bins (0 = no smoothing, default = 2)
0030 %     'nBins'       number of horizontal and vertical bins (default = [50 50])
0031 %     'minTime'     minimum time spent in each bin (in s, default = 0)
0032 %     'mode'        'interpolate' to interpolate missing points (< minTime)
0033 %                   (default for a continuous z) or 'discard' to discard them
0034 %                   (default for a point process)
0035 %     'maxDistance' maximal distance for interpolation (default = 5)
0036 %     'maxGap'      z values recorded during time gaps between successive (x,y)
0037 %                   samples exceeding this threshold (e.g. undetects) will not
0038 %                   be interpolated; also, such long gaps in (x,y) sampling
0039 %                   will be clipped to 'maxGap' to compute the occupancy map
0040 %                   (default = 0.100 s)
0041 %     'type'        three letters (one for X, one for Y and one for Z) indi-
0042 %                   cating which coordinates are linear ('l') and which are
0043 %                   circular ('c') - for 1D data, only two letters are used
0044 %                   (default 'lll')
0045 %    =========================================================================
0046 %
0047 %  OUTPUT
0048 %
0049 %    map.x          x bins
0050 %    map.y          y bins
0051 %    map.z          average map (z continuous)
0052 %                   or rate map (z point process)
0053 %    map.count      count map (z point process)
0054 %    map.time       occupancy map (in s)
0055 %
0056 %  NOTES
0057 %
0058 %    x values are arranged in columns and y values in rows in all output matrices
0059 %    (e.g. 'map.z').
0060 %
0061 %  SEE
0062 %
0063 %    See also MapStats, FiringMap, PlotColorMap, Accumulate.
0064 
0065 % Copyright (C) 2002-2018 by Michaël Zugaro, 2018 by Ralitsa Todorova
0066 %
0067 % This program is free software; you can redistribute it and/or modify
0068 % it under the terms of the GNU General Public License as published by
0069 % the Free Software Foundation; either version 3 of the License, or
0070 % (at your option) any later version.
0071 
0072 % Check number of parameters
0073 if nargin < 2 | mod(length(varargin),2) ~= 0,
0074   error('Incorrect number of parameters (type ''help <a href="matlab:help Map">Map</a>'' for details).');
0075 end
0076 
0077 % Check parameter sizes
0078 if size(v,2) < 2,
0079     error('Parameter ''[t x y]'' should have at least 2 columns (type ''help <a href="matlab:help Map">Map</a>'' for details).');
0080 end
0081 if (size(z,2) < 1 || size(z,2) > 2) && ~isempty(z),
0082     error('Parameter ''z'' should have 1 or 2 columns (type ''help <a href="matlab:help Map">Map</a>'' for details).');
0083 end
0084 
0085 % Default values
0086 maxGap = 0.1;
0087 map.x = [];
0088 map.y = [];
0089 map.count = [];
0090 map.time = [];
0091 map.z = [];
0092 smooth = 2;
0093 nBins = 50;
0094 minTime = 0;
0095 type = 'lll';
0096 modePointProcess = 'discard';
0097 modeContinuous = 'interpolate';
0098 maxDistance = 5;
0099 
0100 if isempty(v) || size(v,1) < 2, return; end
0101 
0102 % Some info about x, y and z
0103 pointProcess = (isempty(z) | size(z,2) == 1);
0104 if pointProcess,
0105     mode = modePointProcess;
0106 else
0107     mode = modeContinuous;
0108 end
0109 t = v(:,1);
0110 x = v(:,2);
0111 if size(v,2) >= 3,
0112     y = v(:,3);
0113 else
0114     y = [];
0115 end
0116 
0117 % Parse parameter list
0118 for i = 1:2:length(varargin),
0119     if ~ischar(varargin{i}),
0120         error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help Map">Map</a>'' for details).']);
0121     end
0122     switch(lower(varargin{i})),
0123 
0124         case 'smooth',
0125             smooth = varargin{i+1};
0126             if ~isdvector(smooth,'>=0') || length(smooth) > 2,
0127                 error('Incorrect value for property ''smooth'' (type ''help <a href="matlab:help Map">Map</a>'' for details).');
0128             end
0129 
0130         case 'nbins',
0131             nBins = varargin{i+1};
0132             if ~isivector(nBins,'>0') || length(nBins) > 2,
0133                 error('Incorrect value for property ''nBins'' (type ''help <a href="matlab:help Map">Map</a>'' for details).');
0134             end
0135 
0136         case 'mintime',
0137             minTime = varargin{i+1};
0138             if ~isdscalar(minTime,'>=0'),
0139                 error('Incorrect value for property ''minTime'' (type ''help <a href="matlab:help Map">Map</a>'' for details).');
0140             end
0141 
0142         case 'maxgap',
0143             maxGap = varargin{i+1};
0144             if ~isdscalar(maxGap,'>=0'),
0145             error('Incorrect value for property ''maxGap'' (type ''help <a href="matlab:help Map">Map</a>'' for details).');
0146             end
0147 
0148         case 'maxdistance',
0149             maxDistance = varargin{i+1};
0150             if ~isdscalar(maxDistance,'>=0'),
0151             error('Incorrect value for property ''maxDistance'' (type ''help <a href="matlab:help Map">Map</a>'' for details).');
0152             end
0153 
0154         case 'mode',
0155             mode = lower(varargin{i+1});
0156             if ~isastring(mode,'interpolate','discard'),
0157                 error('Incorrect value for property ''mode'' (type ''help <a href="matlab:help Map">Map</a>'' for details).');
0158             end
0159 
0160         case 'type',
0161             type = lower(varargin{i+1});
0162             if (isempty(y) && ~isastring(type,'cc','cl','lc','ll')) || (~isempty(y) && ~isastring(type,'ccl','cll','lcl','lll','ccc','clc','lcc','llc')),
0163                 error('Incorrect value for property ''type'' (type ''help <a href="matlab:help Map">Map</a>'' for details).');
0164             end
0165 
0166         otherwise,
0167             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help Map">Map</a>'' for details).']);
0168 
0169   end
0170 end
0171 
0172 % Make sure x and y are normalized
0173 if max(x) > 1 || min(x) < 0,
0174     x = ZeroToOne(x);
0175     warning('Parameter ''x'' should contain values in [0 1]. The data will now be transformed accordingly.');
0176 end
0177 if ~isempty(y),
0178     if max(y) > 1 || min(y) < 0,
0179         y = ZeroToOne(y);
0180         warning('Parameter ''y'' should contain values in [0 1]. The data will now be transformed accordingly.');
0181     end
0182 end
0183 
0184 % Number of bins for x and y
0185 nBinsX = nBins(1);
0186 if length(nBins) == 1,
0187     nBinsY = nBinsX;
0188     nBins(2) = nBins;
0189 else
0190     nBinsY = nBins(2);
0191 end
0192 
0193 % Bin x and y
0194 x = Bin(x,[0 1],nBinsX);
0195 if ~isempty(y),
0196     y = Bin(y,[0 1],nBinsY);
0197 end
0198 
0199 % Duration for each (X,Y) sample (clipped to maxGap)
0200 dt = diff(t);dt(end+1)=dt(end);dt(dt>maxGap) = maxGap;
0201 
0202 if pointProcess,
0203     % Count occurrences for each (x,y) timestamp
0204     n = CountInIntervals(z,[t t+dt]);
0205 else
0206     % Interpolate z at (x,y) timestamps
0207     [z,discarded] = Interpolate(z,t,'maxGap',maxGap);
0208     if isempty(z), return; end
0209     if strcmp(type(end),'c'),
0210         range = isradians(z(:,2));
0211         z(:,2) = exp(j*z(:,2));
0212     end
0213     % Interpolation may have discarded values from z (due to gaps): remove them from x, y and dt too
0214     if any(discarded),
0215         x(discarded,:) = [];
0216         dt(discarded,:) = [];
0217         if ~isempty(y), y(discarded) = []; end
0218     end
0219     n = 1;
0220 end
0221 
0222 % Computations
0223 if isempty(y),
0224     % 1D (only x)
0225     map.x = linspace(0,1,nBinsX);
0226     map.count = Accumulate(x,n,nBinsX);
0227     map.time = Accumulate(x,dt,nBinsX);
0228     valid = map.time > minTime;
0229     map.count = Smooth(Interpolate1(map.x,map.count,valid,mode,maxDistance),smooth,'type',type(1))';
0230     map.time = Smooth(Interpolate1(map.x,map.time,valid,mode,maxDistance),smooth,'type',type(1))';
0231     if pointProcess,
0232         map.z = map.count./(map.time+eps);
0233     else
0234         map.z = Accumulate(x,z(:,2),nBinsX);
0235         map.z = Smooth(Interpolate1(map.x,map.z,valid,mode,maxDistance),smooth,'type',type(1))';
0236         map.z = map.z./(map.count+eps);
0237     end
0238 else
0239     % 2D (x and y)
0240     map.x = linspace(0,1,nBinsX);
0241     map.y = linspace(0,1,nBinsY);
0242     map.count = Accumulate([x y],n,nBins);
0243     map.time = Accumulate([x y],dt,nBins);
0244     valid = map.time > minTime;
0245     map.count = Smooth(Interpolate2(map.x,map.y,map.count,valid,mode,maxDistance),smooth,'type',type(1:2))';
0246     map.time = Smooth(Interpolate2(map.x,map.y,map.time,valid,mode,maxDistance),smooth,'type',type(1:2))';
0247     if pointProcess,
0248         map.z = map.count./(map.time+eps);
0249     else
0250         map.z = Accumulate([x y],z(:,2),nBins);
0251         map.z = Smooth(Interpolate2(map.x,map.y,map.z,valid,mode,maxDistance),smooth,'type',type(1:2)).';
0252         map.z = map.z./(map.count+eps);
0253     end
0254 end
0255 
0256 % Circular z
0257 if strcmp(type(end),'c'), map.z = wrap(angle(map.z),range); end
0258 
0259 
0260 % Interpolate or discard regions with insufficient sampling
0261 if strcmp(mode,'discard'),
0262     map.z(map.time<=minTime) = 0;
0263 end
0264 
0265 % ------------------------------- Helper functions -------------------------------
0266 
0267 % Interpolate if required (1D)
0268 function yint = Interpolate1(x,y,valid,mode,maxDistance)
0269 
0270 if strcmp(mode,'discard'),
0271     yint = y;
0272 else
0273     yint = interp1(x(valid),y(valid),x);
0274 end
0275 
0276 % Interpolate if required (2D)
0277 function zint = Interpolate2(x,y,z,valid,mode,maxDistance)
0278 
0279 if strcmp(mode,'discard'),
0280     % In discard mode, do nothing
0281     zint = z;
0282 else
0283     % In interpolation mode, interpolate missing points (where time < minTime) using other points
0284     d = DistanceTransform(valid);
0285     xx = repmat(x',1,length(y));
0286     yy = repmat(y,length(x),1);
0287     if exist('scatteredInterpolant') == 2,
0288         F = scatteredInterpolant(xx(d==0),yy(d==0),z(d==0));
0289         zint = F(xx,yy);
0290     else
0291         if any(imag(z(:))),
0292             Freal = TriScatteredInterp(xx(d==0),yy(d==0),real(z(d==0)));
0293             zintReal = Freal(xx,yy);
0294             Fimaginary = TriScatteredInterp(xx(d==0),yy(d==0),imag(z(d==0)));
0295             zintImaginary = Fimaginary(xx,yy);
0296             zint = zintReal + 1i.*zintImaginary;
0297         else
0298             F = TriScatteredInterp(xx(d==0),yy(d==0),z(d==0));
0299             zint = F(xx,yy);
0300         end
0301     end
0302     % (do not interpolate missing points too distant from valid points)
0303     zint(d>maxDistance) = z(d>maxDistance);
0304     zint(isnan(zint)) = z(isnan(zint));
0305 end
0306