From faed4994f0f1ce2fa4e597071781dde62b02e9b2 Mon Sep 17 00:00:00 2001 From: ecarmichael Date: Thu, 6 Dec 2018 16:51:01 -0500 Subject: [PATCH] added ACF sandbox which also contains gamma event detection and gamma vs non-gamma comparison. Also added req functions to shared from vdmlab codebase --- shared/io/ProcessConfig2.m | 51 ++++ shared/io/datatypes/iv/AddNActiveCellsIV.m | 73 ++++++ shared/io/datatypes/iv/AddTSDtoIV.m | 61 +++++ shared/io/datatypes/iv/CheckIV.m | 50 ++++ shared/io/datatypes/iv/DifferenceIV.m | 62 +++++ shared/io/datatypes/iv/IVcenters.m | 8 + shared/io/datatypes/iv/IntersectIV.m | 51 ++++ shared/io/datatypes/iv/InvertIV.m | 11 + shared/io/datatypes/iv/MergeIV.m | 93 ++++++++ shared/io/datatypes/iv/PlotTSDfromIV.m | 116 +++++++++ shared/io/datatypes/iv/PlotTSDfromIV_TFR.m | 98 ++++++++ shared/io/datatypes/iv/ResizeIV.m | 79 +++++++ shared/io/datatypes/iv/SelectIV.m | 170 ++++++++++++++ shared/io/datatypes/iv/TSDtoIV.m | 121 ++++++++++ shared/io/datatypes/iv/UnionIV.m | 35 +++ shared/io/datatypes/iv/iv.m | 75 ++++++ shared/io/datatypes/plot.m | 31 +++ shared/io/datatypes/tc/AnotherTCplot.m | 71 ++++++ shared/io/datatypes/tc/AnotherTCplotTwin.m | 93 ++++++++ shared/io/datatypes/tc/CheckTC.m | 88 +++++++ shared/io/datatypes/tc/PlotTC.m | 61 +++++ shared/io/datatypes/tc/TCPlot.m | 108 +++++++++ shared/io/datatypes/tc/tc.m | 54 +++++ shared/io/datatypes/tc/tccmp.m | 70 ++++++ shared/io/datatypes/ts/AddRateTS.m | 100 ++++++++ shared/io/datatypes/ts/CheckTS.m | 75 ++++++ shared/io/datatypes/ts/OrderSelectS.m | 65 ++++++ shared/io/datatypes/ts/SelectTS.m | 164 +++++++++++++ shared/io/datatypes/ts/ShuffleTS.m | 58 +++++ shared/io/datatypes/ts/concatenateTS.m | 65 ++++++ shared/io/datatypes/ts/firstSpike.m | 7 + shared/io/datatypes/ts/getISI.m | 55 +++++ shared/io/datatypes/ts/getSpikeCount.m | 42 ++++ shared/io/datatypes/ts/lastSpike.m | 7 + shared/io/datatypes/ts/removeEmptyCells.m | 13 ++ shared/io/datatypes/ts/ts.m | 71 ++++++ shared/io/datatypes/ts/tscmp.m | 49 ++++ shared/io/datatypes/tsd/CheckTSD.m | 60 +++++ shared/io/datatypes/tsd/MergeTSD.m | 102 ++++++++ shared/io/datatypes/tsd/TSD_SelectChannel.m | 25 ++ shared/io/datatypes/tsd/TSD_getidx.m | 50 ++++ shared/io/datatypes/tsd/TSD_getidx2.m | 81 +++++++ shared/io/datatypes/tsd/TSDbySpace.m | 38 +++ shared/io/datatypes/tsd/TSDbyZ.m | 33 +++ shared/io/datatypes/tsd/TSDpeth.m | 129 +++++++++++ shared/io/datatypes/tsd/TSDtoFT.m | 53 +++++ shared/io/datatypes/tsd/UnionTSD.m | 88 +++++++ shared/io/datatypes/tsd/decimate_tsd.m | 29 +++ shared/io/datatypes/tsd/tsd.m | 57 +++++ shared/io/datatypes/tsd/zscore_tsd.m | 23 ++ shared/proc/FilterLFP.m | 129 +++++++++++ workflows/sandbox_ACF_gamma_events.m | 245 ++++++++++++++++++++ workflows/sandbox_spikeSpectra.m | 29 +-- 53 files changed, 3658 insertions(+), 14 deletions(-) create mode 100644 shared/io/ProcessConfig2.m create mode 100644 shared/io/datatypes/iv/AddNActiveCellsIV.m create mode 100644 shared/io/datatypes/iv/AddTSDtoIV.m create mode 100644 shared/io/datatypes/iv/CheckIV.m create mode 100644 shared/io/datatypes/iv/DifferenceIV.m create mode 100644 shared/io/datatypes/iv/IVcenters.m create mode 100644 shared/io/datatypes/iv/IntersectIV.m create mode 100644 shared/io/datatypes/iv/InvertIV.m create mode 100644 shared/io/datatypes/iv/MergeIV.m create mode 100644 shared/io/datatypes/iv/PlotTSDfromIV.m create mode 100644 shared/io/datatypes/iv/PlotTSDfromIV_TFR.m create mode 100644 shared/io/datatypes/iv/ResizeIV.m create mode 100644 shared/io/datatypes/iv/SelectIV.m create mode 100644 shared/io/datatypes/iv/TSDtoIV.m create mode 100644 shared/io/datatypes/iv/UnionIV.m create mode 100644 shared/io/datatypes/iv/iv.m create mode 100644 shared/io/datatypes/plot.m create mode 100644 shared/io/datatypes/tc/AnotherTCplot.m create mode 100644 shared/io/datatypes/tc/AnotherTCplotTwin.m create mode 100644 shared/io/datatypes/tc/CheckTC.m create mode 100644 shared/io/datatypes/tc/PlotTC.m create mode 100644 shared/io/datatypes/tc/TCPlot.m create mode 100644 shared/io/datatypes/tc/tc.m create mode 100644 shared/io/datatypes/tc/tccmp.m create mode 100644 shared/io/datatypes/ts/AddRateTS.m create mode 100644 shared/io/datatypes/ts/CheckTS.m create mode 100644 shared/io/datatypes/ts/OrderSelectS.m create mode 100644 shared/io/datatypes/ts/SelectTS.m create mode 100644 shared/io/datatypes/ts/ShuffleTS.m create mode 100644 shared/io/datatypes/ts/concatenateTS.m create mode 100644 shared/io/datatypes/ts/firstSpike.m create mode 100644 shared/io/datatypes/ts/getISI.m create mode 100644 shared/io/datatypes/ts/getSpikeCount.m create mode 100644 shared/io/datatypes/ts/lastSpike.m create mode 100644 shared/io/datatypes/ts/removeEmptyCells.m create mode 100644 shared/io/datatypes/ts/ts.m create mode 100644 shared/io/datatypes/ts/tscmp.m create mode 100644 shared/io/datatypes/tsd/CheckTSD.m create mode 100644 shared/io/datatypes/tsd/MergeTSD.m create mode 100644 shared/io/datatypes/tsd/TSD_SelectChannel.m create mode 100644 shared/io/datatypes/tsd/TSD_getidx.m create mode 100644 shared/io/datatypes/tsd/TSD_getidx2.m create mode 100644 shared/io/datatypes/tsd/TSDbySpace.m create mode 100644 shared/io/datatypes/tsd/TSDbyZ.m create mode 100644 shared/io/datatypes/tsd/TSDpeth.m create mode 100644 shared/io/datatypes/tsd/TSDtoFT.m create mode 100644 shared/io/datatypes/tsd/UnionTSD.m create mode 100644 shared/io/datatypes/tsd/decimate_tsd.m create mode 100644 shared/io/datatypes/tsd/tsd.m create mode 100644 shared/io/datatypes/tsd/zscore_tsd.m create mode 100644 shared/proc/FilterLFP.m create mode 100755 workflows/sandbox_ACF_gamma_events.m diff --git a/shared/io/ProcessConfig2.m b/shared/io/ProcessConfig2.m new file mode 100644 index 0000000..e1e698a --- /dev/null +++ b/shared/io/ProcessConfig2.m @@ -0,0 +1,51 @@ +function [cfg_out] = ProcessConfig2(cfg_default,cfg_in) +%% PROCESSCONFIG2 +% Processes the cfg file and the cfg_in function defaults as ProcessConfig does but as a +% function. +% +% INPUTS: +% cfg_default: default parameters stored as a struct within the containing function +% cfg_in: input parameters from the input cfg file +% +% OUTPUT: +% Returns a new cfg file containing the default parameters and any updated cfg +% parameters given the input file. +% +% youkitan 2014-11-04 +% +%% Set output cfg to default cfg + +cfg_out = cfg_default; + +%% Process input cfg parameters into workspace +% NOTE - The two processes (replace and add) are separated only for clarity +% and ease of debugging. It can easily be compressed to half the amount of code (and +% therefore time). + +if ~isempty(cfg_in) + + % If there are default parameters,replace them with cfg_in input + if ~isempty(cfg_default) + default_F = fieldnames(cfg_default); + numF = length(default_F); + + for i = 1:numF + idF = default_F{i}; + if isfield(cfg_in,idF) + cfg_out.(idF) = cfg_in.(idF); + end %set cfg_out fields + end %iterate default cfg fields + end + + % Add input parameters + new_F = fieldnames(cfg_in); + numF = length(new_F); + + for i = 1:numF + inF = new_F{i}; + if ~isfield(cfg_default,inF) + cfg_out.(inF) = cfg_in.(inF); + end %set cfg_out fields + end %iterate extra cfg fields + +end diff --git a/shared/io/datatypes/iv/AddNActiveCellsIV.m b/shared/io/datatypes/iv/AddNActiveCellsIV.m new file mode 100644 index 0000000..0270459 --- /dev/null +++ b/shared/io/datatypes/iv/AddNActiveCellsIV.m @@ -0,0 +1,73 @@ +function iv_in = AddNActiveCellsIV(cfg_in,iv_in,S) +% function iv_out = AddNActiveCellsIV(cfg_in,iv_in,S) +% +% INPUT: +% cfg_in - input cfg parameters +% iv_in - input iv +% S - spikes +% +% OUTPUT: +% iv_out - iv with nActiveCells field added +% +% CFG OPTIONS: +% +% cfg_def.label = 'nActiveCells'; What to call the field that +% contains this data. +% cfg_def.dt = 0.001; % bin size for binning spikes +% +% mvdm 2015-02-01 initial version +% note: uses binning method for speed, i.e. is not exact + +% parse cfg parameters +cfg_def = []; +cfg_def.label = 'nActiveCells'; +cfg_def.dt = 0.005; +cfg_def.verbose = 1; + +mfun = mfilename; +cfg = ProcessConfig(cfg_def,cfg_in,mfun); + + +% check inputs +if ~CheckIV(iv_in) + error('Interval data must have been made with the iv constructor') +end + +if isempty(iv_in.tstart) + error('Interval data is empty') +end + +if cfg.verbose; disp([mfun,': Adding the number of active units during each interval']); end + +% make matrix of spike counts (nCells x nTimeBins) +cfg_temp.dt = cfg.dt; +cfg_temp.verbose = 0; +Q = MakeQfromS(cfg_temp,S); + +% get number of active cells for each +for iIV = length(iv_in.tstart):-1:1 + + Qr = restrict(Q,iv_in.tstart(iIV),iv_in.tend(iIV)); + + if ~isempty(Qr.tvec) + spk_counts = sum(Qr.data,2); + nC(iIV) = sum(spk_counts > 0); + else % no spikes in this iv + nC(iIV) = 0; + end + +end + +if ~isfield(iv_in,'usr') + iv_in.usr = []; +end + +if cfg.verbose && isfield(iv_in.usr,cfg.label) + disp(['WARNING in ',mfun,': iv data already includes usr.',cfg.label,', overwriting...']) +end +iv_in.usr.(cfg.label) = nC; + + +% housekeeping +iv_in.cfg.history.mfun = cat(1,iv_in.cfg.history.mfun,mfun); +iv_in.cfg.history.cfg = cat(1,iv_in.cfg.history.cfg,{cfg}); \ No newline at end of file diff --git a/shared/io/datatypes/iv/AddTSDtoIV.m b/shared/io/datatypes/iv/AddTSDtoIV.m new file mode 100644 index 0000000..0601701 --- /dev/null +++ b/shared/io/datatypes/iv/AddTSDtoIV.m @@ -0,0 +1,61 @@ +function iv = AddTSDtoIV(cfg_in,iv,tsd_in) +% function iv_out = AddTSDtoIV(cfg_in,iv_in,tsd_in) +% +% add usr field to iv based on tsd +% +% INPUTS: +% +% iv1: interval data to be selected from based on.. +% iv2: interval data to be used to select from iv1 +% +% CFG OPTIONS: +% cfg.method = 'max'; 'min', 'mean' +% cfg.target = []; % which label in tsd to use +% cfg.label = []; % what to call this in iv, i.e. usr.label +% +% OUTPUTS: +% +% iv_out: output interval data +% +% MvdM 2014-06-24 + +cfg_def.method = 'max'; + +cfg = ProcessConfig(cfg_def,cfg_in); % should take whatever is in cfg_in and put it into cfg! +mfun = mfilename; + +% check for well formed input +nData = size(tsd_in.data,1); +if nData == 1 + data_in = tsd_in.data; +else + if isempty(cfg.target) + error('Multiple data dimensions present, must specify cfg.target.'); + else + data_in = getd(tsd_in,cfg.target); + end +end + +% find indices for iv +tstart_idx = nearest_idx3(iv.tstart,tsd_in.tvec); +tend_idx = nearest_idx3(iv.tend,tsd_in.tvec); + +% collect data +data_temp = nan(size(tstart_idx)); + +switch cfg.method + case 'max' + for iI = 1:length(data_temp) + data_temp(iI) = max(data_in(tstart_idx(iI):tend_idx(iI))); + end + case 'mean' + for iI = 1:length(data_temp) + data_temp(iI) = mean(data_in(tstart_idx(iI):tend_idx(iI))); + end +end + +iv.usr.(cfg.label) = data_temp; + +% housekeeping +iv.cfg.history.mfun = cat(1,iv.cfg.history.mfun,mfun); +iv.cfg.history.cfg = cat(1,iv.cfg.history.cfg,{cfg}); \ No newline at end of file diff --git a/shared/io/datatypes/iv/CheckIV.m b/shared/io/datatypes/iv/CheckIV.m new file mode 100644 index 0000000..5294c66 --- /dev/null +++ b/shared/io/datatypes/iv/CheckIV.m @@ -0,0 +1,50 @@ +function pass_flag = CheckIV(iv_in,varargin) +%% CHECKIV Check IV for datatype violations +% pass_flag = CheckIV(iv_in,varargin) verifies that input is iv and is well formed. +% +% INPUTS: +% iv_in: iv to be checked +% varargins: +% callername: name of invoking function (for more informative warning messages) +% +% OUTPUTS: +% pass_flag: 1 if all checks pass, 0 if otherwise +% +% Checks performed: +% - .tstart and .tend fields must exist (FAIL) +% - are either .tstart or .tend empty? (WARNING) +% - both .tstart and .tend must be column vectors (FAIL) +% - are start-end pairs unidirectional? (WARNING) +% +% see also iv, CheckTS, CheckTSD, CheckTC +% +% MvdM 2014-11-12 +% aacarey edit Sept 2015, additional checks +% youkitan edit Sept 2016, additional checks +% youkitan edit Dec 2016, reformat help, add function name to output + +pass_flag = 1; + +in_mfun = ''; +if ~isempty(varargin) && ischar(varargin{1}) + in_mfun = [' in ',varargin{1}]; +end + +if isstruct(iv_in) + if ~isfield(iv_in,'tstart') || ~isfield(iv_in,'tend') + pass_flag = 0; + fprintf('FAIL%s by CheckIV: input iv must contain tstart and tend fields.\n',in_mfun); + elseif isempty(iv_in.tstart) || isempty(iv_in.tend) + fprintf('WARNING%s by CheckIV: input iv is empty.\n',in_mfun); + elseif ~iscolumn(iv_in.tstart) || ~iscolumn(iv_in.tend) + pass_flag = 0; + fprintf('FAIL%s by CheckIV: tstart and tend must be column vectors.\n',in_mfun); + elseif any(iv_in.tstart > iv_in.tend) + fprintf('WARNING%s by CheckIV: start-end pairs not unidirectional.\n',in_mfun); + end +else + pass_flag = 0; + fprintf('FAIL%s by CheckIV: input must be an iv data type.\n',in_mfun); +end + + diff --git a/shared/io/datatypes/iv/DifferenceIV.m b/shared/io/datatypes/iv/DifferenceIV.m new file mode 100644 index 0000000..17ef20b --- /dev/null +++ b/shared/io/datatypes/iv/DifferenceIV.m @@ -0,0 +1,62 @@ +function iv_out = DifferenceIV(cfg_in,iv1,iv2) +% function iv_out = DifferenceIV(cfg,iv1,iv2) +% +% keep only those iv1's that DO NOT include a piece of iv2 +% +% iv1 ____________ _____ ______ ____________ ______ +% iv2 ________ _____ _______ ______ +% +% iv_out ______ ______ +% +% INPUTS: +% +% iv1: interval data to be selected from based on.. +% iv2: interval data to be used to select from iv1 +% +% CFG OPTIONS: +% cfg.verbose = 1; If 1, displays text in command window; if 0, doesn't. +% +% OUTPUTS: +% +% iv_out: output interval data +% +% MvdM 2014-06-28 + +cfg_def.verbose = 1; +mfun = mfilename; + +cfg = ProcessConfig(cfg_def,cfg_in,mfun); % should take whatever is in cfg_in and put it into cfg! + +if ~CheckIV(iv1) || ~CheckIV(iv2); error('Inputs must be iv data type'); end + +keep = ones(length(iv1.tstart),1); +for iI = 1:length(iv1.tstart) + + % first, check if any start or end is within this interval + temp1 = find(iv2.tstart >= iv1.tstart(iI) & iv2.tstart <= iv1.tend(iI)); + temp2 = find(iv2.tend >= iv1.tstart(iI) & iv2.tend <= iv1.tend(iI)); + + if ~isempty(temp1) | ~isempty(temp2) + keep(iI) = 0; + continue; + end + + % check if interval is enveloped by anything + for iJ = 1:length(iv2.tstart) + if iv2.tstart(iJ) < iv1.tstart(iI) & iv2.tend(iJ) > iv1.tend(iI) + keep(iI) = 0; + break; + end + end + +end + +iv_out = iv1; +iv_out.tstart = iv_out.tstart(logical(keep)); +iv_out.tend = iv_out.tend(logical(keep)); + +if cfg.verbose; fprintf('%s: %d intervals in, %d intervals out\n',mfun,length(iv1.tstart) + length(iv2.tstart),length(iv_out.tstart)); end + +% housekeeping +iv1.cfg.history.mfun = cat(1,iv1.cfg.history.mfun,mfun); +iv1.cfg.history.cfg = cat(1,iv1.cfg.history.cfg,{cfg}); \ No newline at end of file diff --git a/shared/io/datatypes/iv/IVcenters.m b/shared/io/datatypes/iv/IVcenters.m new file mode 100644 index 0000000..8faad97 --- /dev/null +++ b/shared/io/datatypes/iv/IVcenters.m @@ -0,0 +1,8 @@ +function c = IVcenters(iv_in) +% returns nIV x 1 vector with center times of iv_in + +if ~CheckIV(iv_in) + return; +end + +c = nanmean(cat(2,iv_in.tstart,iv_in.tend),2); \ No newline at end of file diff --git a/shared/io/datatypes/iv/IntersectIV.m b/shared/io/datatypes/iv/IntersectIV.m new file mode 100644 index 0000000..02c1e01 --- /dev/null +++ b/shared/io/datatypes/iv/IntersectIV.m @@ -0,0 +1,51 @@ +function iv1 = IntersectIV(cfg_in,iv1,iv2) +% function iv = IntersectIV(cfg,iv1,iv2) +% +% keep only those iv1's that include a piece of iv2 +% +% INPUTS: +% +% iv1: interval data to be selected from based on.. +% iv2: interval data to be used to select from iv1 +% +% CFG OPTIONS: +% +% OUTPUTS: +% +% iv_out: output interval data +% +% MvdM 2014-06-24 + +cfg_def = []; + +mfun = mfilename; +cfg = ProcessConfig(cfg_def,cfg_in,mfun); % should take whatever is in cfg_in and put it into cfg! + +keep = zeros(length(iv1.tstart),1); +for iI = 1:length(iv1.tstart) + + % first, check if any start or end is within this interval + temp1 = find(iv2.tstart > iv1.tstart(iI) & iv2.tstart < iv1.tend(iI)); + temp2 = find(iv2.tend > iv1.tstart(iI) & iv2.tend < iv1.tend(iI)); + + if ~isempty(temp1) || ~isempty(temp2) + keep(iI) = 1; + continue; + end + + % check if interval is enveloped by anything + for iJ = 1:length(iv2.tstart) + if iv2.tstart(iJ) < iv1.tstart(iI) && iv2.tend(iJ) > iv1.tend(iI) + keep(iI) = 1; + break; + end + end + +end + +iv1.tstart = iv1.tstart(logical(keep)); +iv1.tend = iv1.tend(logical(keep)); + +% housekeeping +iv1.cfg.history.mfun = cat(1,iv1.cfg.history.mfun,mfun); +iv1.cfg.history.cfg = cat(1,iv1.cfg.history.cfg,{cfg}); \ No newline at end of file diff --git a/shared/io/datatypes/iv/InvertIV.m b/shared/io/datatypes/iv/InvertIV.m new file mode 100644 index 0000000..610d6bb --- /dev/null +++ b/shared/io/datatypes/iv/InvertIV.m @@ -0,0 +1,11 @@ +function iv_out = InvertIV(iv_in,t0,t1) +% function iv_out = InvertIV(iv_in,t0,t1) +% +% "flips" included intervals to non-included, and vice versa +% +% MvdM 2016-04-27 + +iv_out = iv; +iv_out.tstart = cat(1,t0,iv_in.tend); +iv_out.tend = cat(1,iv_in.tstart,t1); + diff --git a/shared/io/datatypes/iv/MergeIV.m b/shared/io/datatypes/iv/MergeIV.m new file mode 100644 index 0000000..13aad40 --- /dev/null +++ b/shared/io/datatypes/iv/MergeIV.m @@ -0,0 +1,93 @@ +function iv_out = MergeIV2(cfg_in,iv_in) +%MERGEIV Merge touching, overlapping, or nearby intervals within an iv struct +% iv_out = MERGEIV(cfg,iv_in) +% +% If an interval's end time is less than or equal to the next interval's +% start time, the two intervals are merged into a single interval. +% +% tstart(1) ______________________tend(1) +% iv_in +% tstart(2)______________________tend(2) +% +% +% iv_out tstart(1) ___________________________________tend(1) +% +% If intervals are less than a user-defined distance away from each other, +% the intervals are merged (see cfg.gap). +% +% tstart(1) ___________tend(1) +% +% iv_in |<---cfg.gap--->| +% +% tstart(2)_____________tend(2) +% +% +% iv_out tstart(1) ___________________________________tend(1) +% +% INPUTS: +% cfg: config struct with fields controlling function behavior +% iv_in: iv struct with intervals in ascending order +% +% CFG OPTIONS: +% cfg.gap = 0; Merge events closer than this +% cfg.verbose = 1; Tell me how many intervals came in, and how many +% went out. +% +% aacarey oct 2015 +% +% see also UnionIV, IntersectIV, ResizeIV + +%% +% set cfg defaults +cfg_def.gap = 0; +cfg_def.verbose = 1; + +if ~CheckIV(iv_in); + error('iv_in must be an iv data type.') +end + +% make sure the intervals are ordered +if ~issorted(iv_in.tstart); + error('Intervals must be in ascending order.') +end + +mfun = mfilename; + +% parse cfg parameters +cfg = ProcessConfig(cfg_def,cfg_in,mfun); + +% add 1/2 cfg.gap to the intervals +cfg_temp.d = [-cfg.gap/2 cfg.gap/2]; +cfg_temp.allowOverlap = 1; +iv_temp = ResizeIV(cfg_temp,iv_in); + +% merge the resulting iv +remove = zeros(size(iv_temp.tstart)); +for iInterval = 1:length(iv_temp.tstart)-1; + if iv_temp.tstart(iInterval+1) <= iv_temp.tend(iInterval) + remove(iInterval) = iInterval; + iv_temp.tstart(iInterval+1) = iv_temp.tstart(iInterval); + end +end + +remove = remove > 0; % kill will be 1 for intervals to remove, and zero for the ones to be kept +iv_temp.tstart(remove) = []; +iv_temp.tend(remove) = []; + +% remove excess from outer flanks (this happens because we're using ResizeIV +% to grow the inner flanks) +cfg_temp =[]; +cfg_temp.d = [cfg.gap/2 -cfg.gap/2]; +iv_out = ResizeIV(cfg_temp,iv_temp); + +% talk to me +if cfg.verbose + disp([mfun,': ',num2str(length(iv_in.tstart)),' intervals in, ',num2str(length(iv_out.tstart)),' intervals out.']) +end + +% housekeeping +iv_out.cfg.history.mfun = cat(1,iv_in.cfg.history.mfun,mfun); +iv_out.cfg.history.cfg = cat(1,iv_in.cfg.history.cfg,{cfg}); + +end + diff --git a/shared/io/datatypes/iv/PlotTSDfromIV.m b/shared/io/datatypes/iv/PlotTSDfromIV.m new file mode 100644 index 0000000..063516a --- /dev/null +++ b/shared/io/datatypes/iv/PlotTSDfromIV.m @@ -0,0 +1,116 @@ +function h = PlotTSDfromIV(cfg_in,iv_in,tsd_in) +% function PlotTSDfromIV(cfg,iv_in,tsd_in) +% +% display intervals defined relative to tsd +% +% INPUTS: +% +% iv_in: input iv +% tsd_in: input tsd +% +% OUTPUTS: +% +% h: struct with handles to LFP and highlighted iv's +% +% CFG OPTIONS: +% cfg_def.display = 'tsd'; % {'tsd','iv'}, if 'tsd' then plot tsd and +% highlight iv's within it; if 'iv' then only plot iv's +% cfg_def.bgcol = 'k'; % tsd color outside iv's +% cfg_def.fgcol = 'r'; % tsd color within iv's +% cfg_def.target = []; % which tsd channel to use if tsd has multiple data dimensions +% +% ONLY USED IN 'IV' DISPLAY MODE: +% cfg_def.mode = 'edges'; % {'edges','center'}: plot each iv either using a +% specific time window ('center') or relative to each iv's edges ('edges') +% cfg_def.width = 0.2; % size of time window in s (for 'center') or time +% window to be added to each edge (for 'edges') +% cfg_def.subplotdim = [10 8]; % for 'iv' display mode only, specifies +% subplots to use in single figure +% cfg_def.title = []; % for each iv, add contents of specified usr +% field as title (SHOULD GENERALIZE TO ARBITRARY STRING WITH MULTIPLE +% FIELDS, FORMATTING...) +% +% MvdM 2014-06-24, edit 2016-01-0y to add title option + +cfg_def = []; +cfg_def.verbose = 1; +cfg_def.display = 'tsd'; % 'iv' +cfg_def.mode = 'edges'; % 'center' +cfg_def.width = 0.2; % in s +cfg_def.subplotdim = [10 8]; +cfg_def.bgcol = 'k'; +cfg_def.fgcol = 'r'; +cfg_def.MarkerSize = 10; +cfg_def.iv_only = 0; % if 1, don't plot tsd +cfg_def.title = []; + +cfg = ProcessConfig(cfg_def,cfg_in); % should take whatever is in cfg_in and put it into cfg! +mfun = mfilename; + +% check if conditions are in place +nData = size(tsd_in.data,1); +if nData > 1 + if ~isempty(cfg.target) + temp_data = getd(tsd_in,cfg.target); + else + error('Multiple data dimensions exist but no label is specified.'); + end +else + temp_data = tsd_in.data; +end + +% find indices for iv +tstart_idx = nearest_idx3(iv_in.tstart,tsd_in.tvec); +tend_idx = nearest_idx3(iv_in.tend,tsd_in.tvec); + +% if iv display mode, also find indices for background +switch cfg.mode + case 'edges' + bg_tstart_idx = nearest_idx3(iv_in.tstart-cfg.width,tsd_in.tvec); + bg_tend_idx = nearest_idx3(iv_in.tend+cfg.width,tsd_in.tvec); + case 'center' + ctr = mean(cat(2,iv_in.tstart,iv_in.tend),2); + bg_tstart_idx = nearest_idx3(ctr-cfg.width/2,tsd_in.tvec); + bg_tend_idx = nearest_idx3(ctr+cfg.width/2,tsd_in.tvec); +end + +switch cfg.display + case 'tsd' % plot tsd with highlighted iv + + if ~cfg.iv_only + h.LFP = plot(tsd_in.tvec,temp_data,cfg.bgcol,'MarkerSize',cfg.MarkerSize); + end + hold on; + + h.LFP_iv = nan(size(tstart_idx)); + for iI = 1:length(tstart_idx) + + h.LFP_iv(iI) = plot(tsd_in.tvec(tstart_idx(iI):tend_idx(iI)),temp_data(tstart_idx(iI):tend_idx(iI)),cfg.fgcol,'MarkerSize',cfg.MarkerSize); + + end + + case 'iv' + + nppf = prod(cfg.subplotdim); + for iI = 1:length(tstart_idx) + + figno = ceil(iI./nppf); + plotno = mod(iI-1,nppf) + 1; + maximize; + + figure(figno); + subtightplot(cfg.subplotdim(1),cfg.subplotdim(2),plotno); + + h.LFP(iI) = plot(tsd_in.tvec(bg_tstart_idx(iI):bg_tend_idx(iI)),temp_data(bg_tstart_idx(iI):bg_tend_idx(iI)),cfg.bgcol,'MarkerSize',1); + hold on; + h.LFP_iv(iI) = plot(tsd_in.tvec(tstart_idx(iI):tend_idx(iI)),temp_data(tstart_idx(iI):tend_idx(iI)),cfg.fgcol,'MarkerSize',1); + + axis off; axis tight; + if ~isempty(cfg.title) + all_usr = iv_in.usr.(cfg.title); + title(all_usr(iI)); + end + + end + +end % switch cfg.display diff --git a/shared/io/datatypes/iv/PlotTSDfromIV_TFR.m b/shared/io/datatypes/iv/PlotTSDfromIV_TFR.m new file mode 100644 index 0000000..8f9dd72 --- /dev/null +++ b/shared/io/datatypes/iv/PlotTSDfromIV_TFR.m @@ -0,0 +1,98 @@ +function PlotTSDfromIV_TFR(cfg_in,iv_in,ft_in) +% function PlotTSDfromIV_TFR(cfg,iv_in,ft_in) +% +% plot ft_in data spectrograms and raw data as defined by iv_in +% +% INPUTS: +% +% ft_in: input fieldtrip data to plot from +% iv: intervals to plot data for +% +% CFG OPTIONS: +% +% cfg.target = []; % if ft_in has multiple data dimensions (labels) +% cfg.twin = [-0.5 0.5]; +% cfg.clim = [lower_lim upper_lim] used to make the caxis in the imagesc plots even. default is [0 500] +% ...and a bunch more +% +% MvdM 2014-11-12 + +cfg_def = []; +cfg_def.twin = [-0.5 0.5]; +cfg_def.dt = 0.01; +cfg_def.method = 'mtmconvol'; +cfg_def.taper = 'hanning'; +cfg_def.foi = 1:1:100; % frequencies of interest +cfg_def.subplotdim = [4 5]; +cfg_def.clim = [0 500]; % sets the caxis for the imagesc plots. default [0 500] + +cfg = ProcessConfig2(cfg_def,cfg_in); % should take whatever is in cfg_in and put it into cfg! + +cfg.t_ftimwin = 15*(1./cfg.foi); +cfg.toi = cfg.twin(1):cfg.dt:cfg.twin(2); + +% check if conditions are in place +nTrials = length(ft_in.trial); +if nTrials > 1 + error('Multiple trials not yet implemented.'); + % MvdM: ..and may not even make sense +end + +nData = size(ft_in.trial{1},1); +if nData > 1 + if ~isempty(cfg.target) + temp_cfg = []; + temp_cfg.channel = cfg.target; + ft_in = ft_selectdata(temp_cfg,ft_in); + else + error('Multiple data dimensions exist but no label is specified.'); + end +end + +% create trials from iv +trl_cfg = []; +trl_cfg.t = IVcenters(iv_in); +trl_cfg.mode = 'ft'; +trl_cfg.hdr = ft_in.hdr; + +trl_cfg.twin = cfg.twin + [-2 2]; % padding for spectrogram +trl_cfg.hdr.FirstTimeStamp = ft_in.time{1}(1); +trl_cfg.tvec = ft_in.time{1}; % For neuralynx data conversion + +trl = ft_maketrl(trl_cfg); + +temp_cfg = []; temp_cfg.trl = trl; +ft_in = ft_redefinetrial(temp_cfg,ft_in); + +% specgram +cfg.output = 'pow'; +cfg.keeptrials = 'yes'; % need this for stats later + +TFR = ft_freqanalysis(cfg, ft_in); + +% plot +temp_cfg = []; +temp_cfg.latency = cfg.twin; +TFR = ft_selectdata(temp_cfg,TFR); +ft_in = ft_selectdata(temp_cfg,ft_in); + +nppf = prod(cfg.subplotdim); +for iI = 1:length(ft_in.trial) + + figno = ceil(iI./nppf); + plotno = mod(iI-1,nppf) + 1; + maximize; + + figure(figno); + subtightplot(cfg.subplotdim(1),cfg.subplotdim(2),plotno); + + imagesc(TFR.time,TFR.freq,sq(TFR.powspctrm(iI,1,:,:))); axis xy + caxis(cfg.clim) + hold on; + + temp_data = rescale(ft_in.trial{iI},cfg.foi(1),cfg.foi(end)); + % note, may need alternative option that doesn't rescale + plot(ft_in.time{iI},temp_data,'LineWidth',1,'Color',[1 1 1]); + + axis off; axis tight; +end \ No newline at end of file diff --git a/shared/io/datatypes/iv/ResizeIV.m b/shared/io/datatypes/iv/ResizeIV.m new file mode 100644 index 0000000..7a96863 --- /dev/null +++ b/shared/io/datatypes/iv/ResizeIV.m @@ -0,0 +1,79 @@ +function iv = ResizeIV(cfg_in,iv) +% ResizeIV expand or contract intervals by specified amount +% function iv = ResizeIV(cfg_in,iv) +% +% INPUTS +% iv - iv struct containing interval data +% +% CONFIG options/defaults: +% +% cfg.d = [0 0]; amount to expand/contract; must be two-element array +% which will be added to tstart and tend iv fields +% respectively. +% ex: [-1 1] expand +% [1 -1] contract +% +% note that intervals can be shifted in cases like these: +% [-1 -1] shift left +% [1 1] shift right +% +% cfg.allowOverlap = 1; if 1, IV output may overlap if cfg.d is large; +% if 0, IVs will touch but not overlap +% +% if cfg.allowOverlap is set to 1 the intervals might overlap: +% tstart(1) ______________________tend(1) +% tstart(2)______________________tend(2) +% +% if cfg.allowOverlap is set to 0 the intervals will touch: +% +% tstart(1) __________________tend(1) +% tstart(2)_________________tend(2) +% +% +% MvdM 2014-08-28 initial version (originally called addIV) +% aacarey edit Mar 2015, added allowOverlap and merge +% aacarey edit Sept 2015, renamed addIV to ResizeIV and removed merge option + +% make sure the input was created using the iv constructor +is_iv = CheckIV(iv); + +if is_iv % proceed + + cfg_def.d = [0 0]; + cfg_def.allowOverlap = 1; + cfg = ProcessConfig2(cfg_def,cfg_in); + + mfun = mfilename; + + if length(cfg.d) == 2 + iv.tstart = iv.tstart + cfg.d(1); + iv.tend = iv.tend + cfg.d(2); + else + error('cfg.d must contain 2 elements.'); + end + + % remove "bad" intervals + iv_diffs = iv.tend-iv.tstart; + keep = find(iv_diffs > 0); + + if length(keep) < length(iv.tend) + warning('%d/%d intervals discarded',length(iv.tend)-length(keep),length(iv.tend)); + iv.tstart = iv.tstart(keep); + iv.tend = iv.tend(keep); + end + + % check/correct overlap + if ~cfg.allowOverlap + for ii = 1:length(iv.tstart)-1 + if iv.tend(ii) > iv.tstart(ii+1) + avg = mean([iv.tend(ii) iv.tstart(ii+1)]); + iv.tend(ii) = avg; + iv.tstart(ii+1) = avg; + end + end + end + + % housekeeping + iv.cfg.history.mfun = cat(1,iv.cfg.history.mfun,mfun); + iv.cfg.history.cfg = cat(1,iv.cfg.history.cfg,{cfg}); +end \ No newline at end of file diff --git a/shared/io/datatypes/iv/SelectIV.m b/shared/io/datatypes/iv/SelectIV.m new file mode 100644 index 0000000..21bab3e --- /dev/null +++ b/shared/io/datatypes/iv/SelectIV.m @@ -0,0 +1,170 @@ +function [iv_out,idx] = SelectIV(cfg_in,iv_in,selectspec) +%SELECTIV Specify intervals to keep. +% [iv_out,idx] = SELECTIV(cfg_in,iv_in,selectspec) +% +% INPUTS: +% cfg: config struct with fields controlling function behavior +% iv_in: iv struct +% selectspec: selection specifics, either: +% - [nx1] double: logical array or indices specifying which +% intervals to keep. +% - string: string specifying which usr field to work with. If +% selectspec is a string, the config options cfg.operation, +% cfg.threshold, and cfg.str apply. +% +% +% OUTPUTS +% iv_out: iv struct with specified intervals selected and all +% corresponding same-length usr trimmed accordingly +% +% CFG OPTIONS +% cfg.operation = '>='; How to perform numerical selection, see +% cfg.threshold. +% '>' - usr data > threshold +% '>=' - usr data >= threshold +% '<' - usr data < threshold +% '<=' - usr data <= threshold +% '=' - usr data = threshold +% cfg.threshold = 0; Set a numerical threshold for keeping intervals. +% This works on numerical usr contents, but can also be +% applied to strings as long as the first character is +% number-convertible: +% If your field contains strings and cfg.str is +% empty, SelectIV assumes that the first character is a +% number (i.e. a rating) and thresholds based on this +% number. An example would be '1, very good', for which +% SelectIV considers the 1 only. +% cfg.str = ''; If your target usr field contains strings that +% are NOT number-convertible, input the string you +% want to select by. If this is not empty, it overrides +% numerical selection. Examples of non-number-convertible +% strings might be 'good' or 'maybe' or 'poor'. +% cfg.verbose = 1; Tell me how many intervals came in, and how many +% went out. +% aacarey Nov 2015 +% +% see also RemoveIV, restrict + +cfg_def.operation = '>='; +cfg_def.threshold = 0; +cfg_def.str = ''; % if this is not empty, it overrides numerical selection +cfg_def.verbose = 1; +mfun = mfilename; + +if ~CheckIV(iv_in,mfun) + error('iv_in must be an iv data type.') +end + +% parse cfg parameters +cfg = ProcessConfig(cfg_def,cfg_in,mfun); + +if isempty(iv_in.tstart) + if cfg.verbose + fprintf('%s: iv_in is empty, returning iv_in\n',mfun) + iv_out = History(iv_in,mfun,cfg); + return + end +end + +% choose which thing to do +if islogical(selectspec) || isnumeric(selectspec) + spec_type = 'log_or_num'; % specifying intervals to keep using a logical or numerical array +elseif ischar(selectspec) + spec_type = 'string'; % specifying a string which corresponds to a usr field name +else + error('selectspec must be a logical array, numeric array of indices, or a string specifying a usr field name.') +end + +switch spec_type + case 'string' + % make sure usr exists + if ~isfield(iv_in,'usr') + error('iv_in requires usr for this type of selection.') + end + % check that the field actually exists and that it's the right length + if ischar(selectspec) && ~isfield(iv_in.usr,selectspec) + error([selectspec,' does not exist.']) + elseif ischar(selectspec) && length(iv_in.usr.(selectspec)) ~= length(iv_in.tstart) + error(['iv_in.usr.',selectspec,' must have the same dimensions as iv_in.tstart.']) + end + + % if the field contains strings, get ratings in numerical form + if isempty(cfg.str) && ~isnumeric(iv_in.usr.(selectspec)(1)) + str_type = 'rating'; % something like '5, or delete' + temp = nan(size(iv_in.usr.(selectspec))); + for ii = 1:length(temp) + temp(ii,1) = str2double(iv_in.usr.(selectspec){ii,1}(1)); % we assume the rating is the first character in the string + end + elseif isempty(cfg.str) && isnumeric(iv_in.usr.(selectspec)) + str_type = 'rating'; % something like '5, or delete' + temp = iv_in.usr.(selectspec); + elseif ~isempty(cfg.str) + str_type = 'description'; % something like 'good' + temp = iv_in.usr.(selectspec); + end + + % do the thing + switch str_type + case 'rating' + switch cfg.operation + case '>' + keep = temp > cfg.threshold; + case '>=' + keep = temp >= cfg.threshold; + case '<' + keep = temp < cfg.threshold; + case '<=' + keep = temp <= cfg.threshold; + case '=' + keep = temp == cfg.threshold; + otherwise + error('Unrecognized cfg.operation') + end + case 'description' + keep = nan(size(temp)); + for iStr = 1:length(temp) + keep(iStr) = strcmp(cfg.str,temp(iStr)); + end + end + + keep = logical(keep); + + case 'log_or_num' + keep = selectspec; + % these config options do not apply in this case, so don't give + % them a value in history since they did not affect the output + cfg.operation = ''; + cfg.threshold = []; + cfg.str = ''; +end + +iv_out = iv_in; +iv_out.tstart = iv_out.tstart(keep); +iv_out.tend = iv_out.tend(keep); + +% also select data from other same-length usr fields +if isfield(iv_out,'usr') && ~isempty(iv_out.usr) + ivfields = fieldnames(iv_out.usr); + for iField = 1:length(ivfields) + iv_out.usr.(ivfields{iField}) = iv_out.usr.(ivfields{iField})(keep); + end +end + +% make idx output +if islogical(keep) + idx = find(keep); +elseif isnumeric(keep) + idx = keep; +end + +% talk to me +if cfg.verbose + disp([mfun,': ',num2str(length(iv_in.tstart)),' intervals in, ',num2str(length(iv_out.tstart)),' intervals out.']) +end + +% housekeeping +iv_out.cfg.history.mfun = cat(1,iv_in.cfg.history.mfun,mfun); +iv_out.cfg.history.cfg = cat(1,iv_in.cfg.history.cfg,{cfg}); + +end % of function + diff --git a/shared/io/datatypes/iv/TSDtoIV.m b/shared/io/datatypes/iv/TSDtoIV.m new file mode 100644 index 0000000..9758660 --- /dev/null +++ b/shared/io/datatypes/iv/TSDtoIV.m @@ -0,0 +1,121 @@ +function [iv_out,thr_out] = TSDtoIV(cfg_in,tsd_in) +% function [iv_out,thr_out] = TSDtoIV(cfg,tsd) +% +% create interval data from tsd by tresholding +% +% INPUTS: +% +% tsd_in: input tsd +% +% CFG OPTIONS: +% cfg.method = 'zscore'; +% cfg.threshold = 5; +% cfg.operation = '>'; % '<', '>' +% cfg.merge_thr = 0.05; % merge events closer than this +% cfg.target = []; % which data (label) to use +% cfg.minlen = 0.05; % minimum interval length +% cfg.verbose = 1; 1 display command window text, 0 don't +% +% OUTPUTS: +% +% iv_out: output interval data +% thr_out: threshold used (in same units as input; helpful if z-scored +% threshold is to be applied later to other data) +% +% MvdM 2014-06-24 + +cfg_def.method = 'zscore'; +cfg_def.threshold = 5; +cfg_def.operation = '>'; % return intervals where threshold is exceeded +cfg_def.merge_thr = 0.05; % merge events closer than this +cfg_def.target = []; +cfg_def.minlen = 0.05; % minimum interval length +cfg_def.verbose = 1; + +mfun = mfilename; +cfg = ProcessConfig(cfg_def,cfg_in,mfun); % should take whatever is in cfg_in and put it into cfg! + +iv_out = iv; % initialize new iv struct + +% check if conditions are in place +nData = size(tsd_in.data,1); +if nData > 1 + if ~isempty(cfg.target) + temp_data = getd(tsd_in,cfg.target); + else + error('Multiple data dimensions exist but no label is specified.'); + end +else + temp_data = tsd_in.data; +end + +% apply transform to data if requested +thr_out = cfg.threshold; +switch cfg.method + case 'zscore' + [temp_data,mu,sigma] = zscore(temp_data); + switch cfg.dcn + case '>' + thr_out = mu + cfg.threshold.*sigma; % store threshold for output argument + case '<' + thr_out = mu - cfg.threshold.*sigma; + end + case 'percentile' + temp_data_sorted = sort(temp_data,'ascend'); + thr_out_idx = round(cfg.threshold.*length(temp_data)); + thr_out = temp_data_sorted(thr_out_idx); + + temp_data = tiedrank(temp_data)./length(temp_data); % assign percentile to each data point + +end + +% detect crossings +switch cfg.operation + case '>' + detec = temp_data > cfg.threshold; + case '<' + detec = temp_data < cfg.threshold; + case 'range' + detec = temp_data > cfg.threshold(1) & temp_data < cfg.threshold(2); +end + +% pad the detection so we can deal with cases where first or last samples are detects +detec = cat(2,0,detec,0); + +dfs = diff(detec); +up_idx = find(dfs == 1); +down_idx = find(dfs == -1) - 1; + +up_t = tsd_in.tvec(up_idx); +down_t = tsd_in.tvec(down_idx); + +% remove doubles +d = up_t(2:end)-down_t(1:end-1); + +merge_idx = find(d < cfg.merge_thr); +down_t(merge_idx) = []; +up_t(merge_idx+1) = []; + +% remove too short +iv_len = down_t-up_t; +keep_idx = iv_len > cfg.minlen; + +iv_out.tstart = up_t(keep_idx); +iv_out.tend = down_t(keep_idx); + +% ensure column vectors +if ~iscolumn(iv_out.tstart) + iv_out.tstart = iv_out.tstart'; +end + +if ~iscolumn(iv_out.tend) + iv_out.tend = iv_out.tend'; +end + +if cfg.verbose + disp([mfun,': ',num2str(length(iv_out.tstart)),' intervals found.']) +end + +% housekeeping +iv_out.cfg.history.mfun = cat(1,iv_out.cfg.history.mfun,mfun); +iv_out.cfg.history.cfg = cat(1,iv_out.cfg.history.cfg,{cfg}); \ No newline at end of file diff --git a/shared/io/datatypes/iv/UnionIV.m b/shared/io/datatypes/iv/UnionIV.m new file mode 100644 index 0000000..4a8b4b5 --- /dev/null +++ b/shared/io/datatypes/iv/UnionIV.m @@ -0,0 +1,35 @@ +function iv1 = UnionIV(cfg_in,iv1,iv2) +% function iv = UnionIV(cfg,iv1,iv2) +% +% union of iv objects +% +% output is resorted based on start times (ascending) +% +% MvdM 2014-08-28 initial version + +cfg_def = []; +cfg = ProcessConfig(cfg_def,cfg_in); + +mfun = mfilename; + +if isempty(iv1) % function should work for empty arguments + iv1 = iv([],[]); +end + +if isempty(iv2) + iv2 = iv([],[]); +end + +if ~CheckIV(iv1) | ~CheckIV(iv2) + error('Malformed IV.'); +end + +iv1.tstart = cat(1,iv1.tstart,iv2.tstart); % note that constructor guarantees column vectors +iv1.tend = cat(1,iv1.tend,iv2.tend); + +[iv1.tstart,sort_idx] = sort(iv1.tstart,'ascend'); +iv1.tend = iv1.tend(sort_idx); + +% housekeeping +iv1.cfg.history.mfun = cat(1,iv1.cfg.history.mfun,mfun); +iv1.cfg.history.cfg = cat(1,iv1.cfg.history.cfg,{cfg}); \ No newline at end of file diff --git a/shared/io/datatypes/iv/iv.m b/shared/io/datatypes/iv/iv.m new file mode 100644 index 0000000..359f0d9 --- /dev/null +++ b/shared/io/datatypes/iv/iv.m @@ -0,0 +1,75 @@ +function iv_out = iv(varargin) +% IV Interval datatype constructor +% Interval, or iv, data is one of the main datatypes. Each iv struct +% contains a set of intervals and accompanying information (in usr field). +% Many functions require or assume that the intervals are non-overlapping. +% +% function iv_out = IV([tstart tend]) +% +% function iv_out = IV(tstart,tend) +% +% INPUTS: +% +% [tstart end]: [nx2] double containing interval start times and end times +% +% OR +% +% tstart: [nx1] double or [1xn] double containing interval start times +% tend: [nx1] double or [1xn] double containing interval end times +% +% OUTPUTS +% iv_out: iv struct with fields: +% .type - 'iv'; datatype identification +% .tstart - [nx1] double containing interval start times +% .tend - [nx1] double containing interval end times +% .usr - [] initialized field that can be used to store data +% or information pertaining to individual intervals +% within the iv struct. For example, signal power can be +% stored as iv1.usr.power = [nx1] double (NOTE: usr data +% must have the same dimensions as the tstart and tend +% fields). +% .cfg - record of the data's history including config +% parameters and functions visited +% +% see also CheckIV +% +% MvdM 2014-06-24 +% aacarey edit, Nov 2015 + +iv_out.type = 'iv'; +iv_out.tstart = []; +iv_out.tend = []; +iv_out.usr = []; + +if nargin == 1 + + if size(varargin{1},2) == 2 + iv_out.tstart = varargin{1}(:,1); + iv_out.tend = varargin{1}(:,2); + else + error('Single input argument must have length 2 (tstart, tend)'); + end + +elseif nargin == 2 + + if numel(varargin{1}) == numel(varargin{2}) + iv_out.tstart = varargin{1}; + iv_out.tend = varargin{2}; + else + error('Input arguments must have same length (tstart, tend)'); + end + +end + +% ensure column vectors +if ~iscolumn(iv_out.tstart) + iv_out.tstart = iv_out.tstart'; +end + +if ~iscolumn(iv_out.tend) + iv_out.tend = iv_out.tend'; +end + +% housekeeping +iv_out.cfg.history.mfun{1} = mfilename; +iv_out.cfg.history.cfg{1} = []; \ No newline at end of file diff --git a/shared/io/datatypes/plot.m b/shared/io/datatypes/plot.m new file mode 100644 index 0000000..40ea2aa --- /dev/null +++ b/shared/io/datatypes/plot.m @@ -0,0 +1,31 @@ +function h_out = plot(varargin) +% overload plot function to accept ts, tsd, iv structs +% +% MvdM 2016-01-20 + +if isstruct(varargin{1}) + if isfield(varargin{1},'type') + switch varargin{1}.type + case 'tsd' + if size(varargin{1}.data,1) == 1 + h = builtin('plot',varargin{1}.tvec,varargin{1}.data,varargin{2:end}); + elseif size(varargin{1}.data,1) == 2 + h = builtin('plot',varargin{1}.data(1,:),varargin{1}.data(2,:),varargin{2:end}); + else + error('plotting for tsds with more than two dimensions not implemented yet.') + end + case 'ts' + h = MultiRaster([],varargin{1}); + otherwise + error('%s plot not implemented yet.',varargin{1}.type); + end + else + error('Unknown or undefined struct type.'); + end +else + h = builtin('plot',varargin{:}); +end + +if nargout == 1 + h_out = h; +end \ No newline at end of file diff --git a/shared/io/datatypes/tc/AnotherTCplot.m b/shared/io/datatypes/tc/AnotherTCplot.m new file mode 100644 index 0000000..7416c68 --- /dev/null +++ b/shared/io/datatypes/tc/AnotherTCplot.m @@ -0,0 +1,71 @@ +function AnotherTCplot(cfg_in,S,TC,pos) +% function AnotherTCplot(cfg_in,S,TC,pos) +% +% plot some linearized tuning tuning curves with their scatterfields +% +% configs: +% +% cfg_def.npf = 6; % number of subplots per figure +% cfg_def.write_output = 0; + +cfg_def.npf = 6; % number of subplots per figure +cfg_def.write_output = 0; + +cfg = ProcessConfig2(cfg_def,cfg_in); + +% should also create a version that plots left and right versions side by +% side in a 4-column, 2-row layout + +for iC = 1:length(S.t) + + % find correct place to plot + fno = floor((iC-1)/cfg.npf) + 1; + sno = mod(iC,cfg.npf); if sno == 0, sno = cfg.npf; end + + fh = figure(fno); set(fh,'Color',[0 0 0]); + + % plot tuning curve + subplot(2,6,1+(sno-1)*2); + + smoo_pfR{iC} = TC.tc(:,iC); + + ah = area(smoo_pfR{iC}); set(ah,'FaceColor',[0 0 0.5],'EdgeColor',[0 0 0.5]); + axis off; set(gca,'YLim',[0 max(30,max(smoo_pfR{iC}))]); hold on; + th = text(1,max(smoo_pfR{iC}),num2str(round(max(smoo_pfR{iC})))); set(th,'Color',[1 1 1]); + + % see if any fields were found, if so plot peaks also + if ~isempty(TC.peak_loc{iC}) + pks.loc = TC.peak_loc{iC}; + pks.val = smoo_pfR{iC}(pks.loc); + + plot(pks.loc,pks.val,'.','MarkerSize',10,'Color',[0 0.7 0]); + else + set(ah,'FaceColor',[0.3 0.3 0.3],'EdgeColor',[0.3 0.3 0.3]); + end + + % plot scatterfield + subplot(2,6,sno*2); + + plot(getd(pos,'y'),getd(pos,'x'),'.','MarkerSize',1,'Color',[0.5 0.5 0.5]); + hold on; axis off + + spk = S.t{iC}; + spk_x = interp1(pos.tvec,getd(pos,'x'),spk,'nearest'); + spk_y = interp1(pos.tvec,getd(pos,'y'),spk,'nearest'); + + h = plot(spk_y,spk_x,'.r','MarkerSize',10); + + [~,fn,fe] = fileparts(S.label{iC}); + fs = cat(2,fn,fe); + + th = title(sprintf('c%d %s',iC,fs)); set(th,'Interpreter','none','Color',[1 1 1],'FontSize',8); + maximize; + + if cfg.write_output & sno == cfg.npf + fn_base = cat(2,S.cfg.SessionID,'-fields'); + fn = cat(2,fn_base,num2str(fno),'.png'); + set(gcf, 'InvertHardCopy', 'off'); + print(gcf,'-r300','-dpng',fn); + end + +end \ No newline at end of file diff --git a/shared/io/datatypes/tc/AnotherTCplotTwin.m b/shared/io/datatypes/tc/AnotherTCplotTwin.m new file mode 100644 index 0000000..c3970d6 --- /dev/null +++ b/shared/io/datatypes/tc/AnotherTCplotTwin.m @@ -0,0 +1,93 @@ +function fh = AnotherTCplotTwin(cfg_in,S,TC1,TC2,pos) +% function fh = AnotherTCplotTwin(cfg_in,S,TC1,TC2,pos) +% +% plot some linearized tuning tuning curves with their scatterfields +% +% INPUTS: +% +% self-explanatory +% +% OUTPUT: +% +% fh: array of figure handles +% +% configs: +% +% cfg_def.npf = 6; % number of subplots per figure +% cfg_def.write_output = 0; +% +% MvdM in a hurry Feb 2015 + +cfg_def.npf = 6; % number of subplots per figure +cfg_def.write_output = 0; +cfg_def.c1 = [1 0 0]; +cfg_def.c2 = [0 1 0]; + +cfg = ProcessConfig2(cfg_def,cfg_in); + +% should also create a version that plots left and right versions side by +% side in a 4-column, 2-row layout + +for iC = 1:length(S.t) + + % find correct place to plot + fno = floor((iC-1)/cfg.npf) + 1; + sno = mod(iC,cfg.npf); if sno == 0, sno = cfg.npf; end + + fh(fno) = figure(fno); set(fh(fno),'Color',[0 0 0]); + + % plot tuning curve + subplot(2,6,1+(sno-1)*2); + + for iT = 1:2 + + if iT == 1, TC = TC1; c = cfg.c1; else TC = TC2; c = cfg.c2; end + + smoo_pfR{iC} = TC.tc(:,iC); + + plot(smoo_pfR{iC},'Color',c,'LineWidth',1); hold on; + set(gca,'Color',[0 0 0],'XColor',[0.3 0.3 0.3],'YColor',[0.3 0.3 0.3],'XTick',[],'YTick',[],'XLim',[1 length(smoo_pfR{iC})]); box off; + + th = text(1,max(smoo_pfR{iC}),num2str(round(max(smoo_pfR{iC})))); set(th,'Color',c); + + % see if any fields were found, if so plot peaks also + f_idx = find(TC.template_idx == iC); + if ~isempty(f_idx) + pks.loc = TC.peak_loc{f_idx}; + pks.val = smoo_pfR{iC}(pks.loc); + + plot(pks.loc,pks.val,'.','MarkerSize',10,'Color',c); + plot(pks.loc,pks.val,'o','MarkerSize',10,'Color',c); + else + %set(ah,'FaceColor',[0.3 0.3 0.3],'EdgeColor',[0.3 0.3 0.3]); + end + + end + yl = ylim; ylim([0 ceil(yl(2))+1]); % make sure we can see peak highlights + + % plot scatterfield + subplot(2,6,sno*2); + + plot(getd(pos,'y'),getd(pos,'x'),'.','MarkerSize',1,'Color',[0.5 0.5 0.5]); + hold on; axis off + + spk = S.t{iC}; + spk_x = interp1(pos.tvec,getd(pos,'x'),spk,'nearest'); + spk_y = interp1(pos.tvec,getd(pos,'y'),spk,'nearest'); + + h = plot(spk_y,spk_x,'.r','MarkerSize',10); + + [~,fn,fe] = fileparts(S.label{iC}); + fs = cat(2,fn,fe); + + th = title(sprintf('c%d %s',iC,fs)); set(th,'Interpreter','none','Color',[1 1 1],'FontSize',8); + maximize; + + if cfg.write_output & sno == cfg.npf + fn_base = cat(2,S.cfg.SessionID,'-fields'); + fn = cat(2,fn_base,num2str(fno),'.png'); + set(gcf, 'InvertHardCopy', 'off'); + print(gcf,'-r300','-dpng',fn); + end + +end \ No newline at end of file diff --git a/shared/io/datatypes/tc/CheckTC.m b/shared/io/datatypes/tc/CheckTC.m new file mode 100644 index 0000000..8fc86d8 --- /dev/null +++ b/shared/io/datatypes/tc/CheckTC.m @@ -0,0 +1,88 @@ +function pass_flag = CheckTC(tc_in,varargin) +%% CHECKTC Check TC for datatype violations +% pass_flag = CheckTC(tc_in,varargin) verifies that input is tc and is well formed. +% +% INPUTS: +% tc_in: tc to be checked +% varargins: +% callername: name of invoking function (for more informative warning messages) +% +% OUTPUTS: +% pass_flag: 1 if all checks pass, 0 if otherwise +% +% Checks performed: +% - .tc or .tc2D field must exist (FAIL) +% - .tc must be an nCells x nBins matrix (FAIL) +% - .tc2D must be an nCells x nXBins x nYBins matrix (FAIL) +% - nBins must agree across variables (FAIL) +% - number of cells must agree across variables (FAIL) +% +% see also: tc, TuningCurves, CheckTSD, CheckTS, CheckIV +% +% youkitan 2016-11-29 initial +% youkitan edit Dec 2016, added some tc2D functionality, reformat help + +pass_flag = 1; + +in_mfun = ''; +if ~isempty(varargin) && ischar(varargin{1}) + in_mfun = [' in ',varargin{1}]; +end + +if ~isstruct(tc_in) + pass_flag = 0; + fprintf('FAIL%s by CheckTC: input must be a tc struct.\n',in_mfun); +else + % basic checks + if isfield(tc_in,'type') + if ~isequal(tc_in.type,'tc') + fprintf('FAIL%s by CheckTC: input must be a tc type.\n',in_mfun); + pass_flag = 0; + + % ~~~~~~~~~~~~~~~~~~~~ checks for all required fields ~~~~~~~~~~~~~~~~~~~~ + elseif ~isfield(tc_in,'occ_hist') + fprintf('FAIL%s by CheckTC: input must have a occupancy field.\n',in_mfun); + pass_flag = 0; + + elseif ~isfield(tc_in,'spk_hist') + fprintf('FAIL%s by CheckTC: input must have a spike count field.\n',in_mfun); + pass_flag = 0; + + elseif ~isfield(tc_in,'binEdges') + fprintf('FAIL%s by CheckTC: input must have bin edge field.\n',in_mfun); + pass_flag = 0; + + % ~~~~~~~~~~~~~~~~~~~~ checks for 1D case ~~~~~~~~~~~~~~~~~~~~ + elseif isfield(tc_in,'tc') + if length(size(tc_in.tc)) ~= 2 + fprintf('FAIL%s by CheckTC: tc variable should be an nCells x nBins matrix.\n',in_mfun); + pass_flag = 0; + + elseif size(tc_in.tc,2) ~= size(tc_in.occ_hist,2) + pass_flag = 0; + fprintf('FAIL%s by CheckTC: number of bins differ between tc and occ_hist.\n',in_mfun); + end + + % ~~~~~~~~~~~~~~~~~~~~ checks for 2D case ~~~~~~~~~~~~~~~~~~~~ + elseif isfield(tc_in,'tc2D') + if length(size(tc_in.tc2D)) ~= 3 + fprintf('FAIL%s by CheckTC: tc variable should be an nCells x nXBins x nYBins matrix.\n',in_mfun); + pass_flag = 0; + end + fprintf('WARNING by CheckTC: Checks for 2D tuning curves not implemented yet!!!') + + % ~~~~~~~~~~~~~~~~~~~~ no tuning variable ~~~~~~~~~~~~~~~~~~~~ + else + pass_flag = 0; + fprintf('FAIL%s by CheckTC: input tc must contain a tc field.\n',in_mfun); + end + + % ~~~~~~~~~~~~~~~~~~~~ no type parameter ~~~~~~~~~~~~~~~~~~~~ + else + fprintf('FAIL%s by CheckTC: input has no type parameter.\n',in_mfun) + pass_flag = 0; + end + +end + +end \ No newline at end of file diff --git a/shared/io/datatypes/tc/PlotTC.m b/shared/io/datatypes/tc/PlotTC.m new file mode 100644 index 0000000..e3b65e5 --- /dev/null +++ b/shared/io/datatypes/tc/PlotTC.m @@ -0,0 +1,61 @@ +function [] = tcplot(cfg_in,tc_in) +%% TCPLOT Tuning Curve Plot +% +% Plots the tuning curves from a tc object +% +% INPUTS: +% cfg_in: input cfg file +% tc_in: input tc +% +% CFG OPTIONS: +% cfg.display = 'separate' % 'together' +% +% NOTES: +% Currently only plots 1D tuning curves. +% +% youkitan 2014-11-05 +%% set defaults and process cfg +cfg_def.display = 'separate'; + +mfun = mfilename; +cfg = ProcessConfig(cfg_def,cfg_in,mfun); + +%% create plot layout +numPF = length(tc_in.firing_rates); +[nsp,~] = numSubplots(numPF); + +%% plot +cases = {'separate','together'}; +assert(any(strcmp(cfg.display,cases)),... + '"%s" is an invalid plot type!',cfg.display) + +switch cfg.display + case 'separate' + indPF = 1; + hold on; + for i = 1:(nsp(1)*nsp(2)) + if indPF <= numPF + subplot(nsp(1),nsp(2),indPF); + plot(tc_in.stimulus{indPF},tc_in.firing_rates{indPF}); + title(['Cell ',num2str(i)],'Fontsize',8); + xlabel('Position (cm)','Fontsize',8); + ylabel('Firing Rate (Hz)','Fontsize',8); + set(gca,'Fontsize',6) + indPF = indPF + 1; + end + end%iterate rows + + case 'together' + plotcolors = {'c','r','g','b','k'}; + hold on; + for i = 1:numPF + cval = mod(i,5) + 1; + plot(tc_in.stimulus{i},tc_in.firing_rates{i},'Color',plotcolors{cval}); + xlabel('Position (cm)','Fontsize',8); + ylabel('Firing Rate (Hz)','Fontsize',8); + set(gca,'Fontsize',6) + end %iterate rows +end %end switch + + +end \ No newline at end of file diff --git a/shared/io/datatypes/tc/TCPlot.m b/shared/io/datatypes/tc/TCPlot.m new file mode 100644 index 0000000..821360e --- /dev/null +++ b/shared/io/datatypes/tc/TCPlot.m @@ -0,0 +1,108 @@ +function tch = TCPlot(cfg_in,TC) +% function TCPlot(cfg,TC) +% +% plot tuning curves, ordered by peak location +% +% INPUTS +% +% TC: tuning curves object (as returned by MakeTC()) +% +% OUTPUTS +% +% tch: vector of handles to plotted tuning curves (so can do things like, +% set(tch,'FaceColor',[1 0 0]) +% +% CONFIGS +% +% cfg_def.ax = []; % if specified, plot here; otherwise create new figure +% cfg_def.order = 1; % if 1, order by template_idx; 2, order by field_template_idx +% cfg_def.mode = 'area'; % 'line', 'area', 'heat' +% cfg_def.color = [0 0 0]; +% cfg_def.alpha = 0; +% cfg_def.binsize = 1; % how many cm in each bin +% cfg_def.cp = []; % location of choice point +% cfg_def.YOrd = 'normal'; % order of Y axis plotting, 'normal' in +% decreasing order or 'reverse' in increasing order. +% cfg_def.XDir = 'normal'; % 'normal', 'reverse' to flip horizontal +% +% disclaimer: +% +% MvdM rough version copied from youkitan workflow chunk; many improvements +% to make especially options to plot all tuning curves, only those with +% fields, etc.. +% aacarey edit sept 2015 for XDir and YOrd + +cfg_def = []; +cfg_def.ax = []; % if specified, plot here; otherwise create new figure +cfg_def.order = 1; % if 1, order by template_idx; 2, order by field_template_idx +cfg_def.mode = 'area'; % 'line' +cfg_def.color = [0 0 0]; +cfg_def.alpha = 0; +cfg_def.binsize = 1; +cfg_def.cp = []; +cfg_def.YOrd = 'normal'; % order of plotting along Y, 'normal' in decreasing order or 'reverse' in increasing order +cfg_def.XDir = 'normal'; % 'normal' 'reverse' + +cfg = ProcessConfig2(cfg_def,cfg_in); + +% +if isempty(cfg.ax) + figure; +else + axes(cfg.ax); +end + +% +switch cfg.order + case 1 + tc_temp = TC.tc(:,TC.template_idx)'; + case 2 + tc_temp = TC.tc(:,TC.field_template_idx)'; + otherwise + tc_temp = TC.tc'; +end + +switch cfg.YOrd % handle reverse Y axis plotting of tuning curves, if specified (aacarey) + case 'normal' + cellNum = 1:size(tc_temp,1); + cellNum = strtrim(cellstr(num2str(cellNum'))'); % convert to cell array of strings cellNum = {'3' '2' '1'} + case 'reverse' + cellNum = size(tc_temp,1):-1:1; % we flip the cell numbering and plot order for cell 1 to be at bottom + cellNum = strtrim(cellstr(num2str(cellNum'))'); % convert to cell array of strings cellNum = {'1' '2' '3'} + tc_temp = flipud(tc_temp); % flip the data up-down for reverse plot order along y axis +end + +x_end = size(tc_temp,2); + +clear tch; +for iC = 1:size(tc_temp,1) + + subaxis(size(tc_temp,1),1,iC,'SpacingVert',0); + + switch cfg.mode + case 'area' + tch(iC) = area(tc_temp(iC,:)); set(tch(iC),'FaceColor',cfg.color,'EdgeColor',[1 1 1]); + if cfg.alpha > 0 % attempt at transparency doesn't work + patchobjs = findobj(get(tch(iC),'children'), '-property', 'AlphaData'); + set(patchobjs, 'AlphaData', cfg.alpha); + end + case 'line' + tch(iC) = plot(tc_temp(iC,:),'Color',cfg.color,'LineWidth',1); + end + hold on; + + %ylims = get(gca,'Ylim'); + %plot([TC.peak_idx(iC) TC.peak_idx(iC)],[ylims(1) ylims(2)],'r:'); + + set(gca, 'XTick', [],'YTick',[],'XLim',[1 x_end],'LineWidth',0.5,'XColor',[0.5 0.5 0.5],'YColor',[0.5 0.5 0.5],'XDir',cfg.XDir); + box off; + + yl = ylabel(cellNum(iC)); + set(yl,'Rotation',0,'Fontsize',8); + %if iC == 1; title(ENC_data(iT).trial_type); end + %axis off; +end + +set(gca,'XTick',[1 cfg.cp x_end],'XTicklabel',[1 cfg.cp*cfg.binsize x_end*cfg.binsize],'Ticklength', [0 0],'Fontsize',8,'YTickLabel',cellNum); + + diff --git a/shared/io/datatypes/tc/tc.m b/shared/io/datatypes/tc/tc.m new file mode 100644 index 0000000..1f5c042 --- /dev/null +++ b/shared/io/datatypes/tc/tc.m @@ -0,0 +1,54 @@ +function tc_out = tc(varargin) +%% TC TUNING CURVE DATATYPE +% The tc dataype is a container for tuning curves typically used in experimental +% analysis. +% +% tc_out = TC(varargin) is a constructor for the tc (tuning curve) struct. +% +% INPUTS: +% tc/tc2D: firing rate map for each cell (spike count normalized by sampling of +% tuning variable) +% occ_hist: total sampling of tuning variable (occupancy for +% position variable) spk_hist: binned raw firing rates (optionally smoothed) +% +% OUTPUTS +% tc_out: tc struct with fields: +% .tc - [NxM] array, where N is the number of cells and M is the number of +% tuning variable bins. +% .tc2D - [NxMxP] array, where N is the number of cells and M and P are the number +% of tuning variable bins(for 2D tuning_var only) +% .occ_hist - [1xM] array, where M is the number of tuning variable bins +% .spk_hist - [1xM] array, where M is the number of tuning variable bins +% +% See also, TuningCurves, MakeTC, DetectPlaceCells1D +% Workflow example WORKFLOW_PlotOrderedRaster +% +% youkitan init 2014-10-27 +% youkitan edit 2015-09-16 +% youkitan edit Feb 2017 overhaul + +%% Set fields +tc_out.tc = []; %nCells x nBins +tc_out.occ_hist = []; % +tc_out.spk_hist = []; + +switch nargin + case 2 + tc_out.tc = varargin{1}; + tc_out.occ_hist = varargin{2}; + + case 3 + tc_out.tc = varargin{1}; + tc_out.occ_hist = varargin{2}; + tc_out.spk_hist = varargin{3}; + + case 4 + tc_out.tc = varargin{1}; + tc_out.occ_hist = varargin{2}; + tc_out.spk_hist = varargin{3}; + tc_out.usr.binEdges = varargin{4}; +end + +% housekeeping +tc_out.cfg.history.mfun{1} = mfilename; +tc_out.cfg.history.cfg{1} = []; \ No newline at end of file diff --git a/shared/io/datatypes/tc/tccmp.m b/shared/io/datatypes/tc/tccmp.m new file mode 100644 index 0000000..9427272 --- /dev/null +++ b/shared/io/datatypes/tc/tccmp.m @@ -0,0 +1,70 @@ +function pass_flag = tccmp(TC1,TC2) +%% TCCMP Compare Tuning Curve objects +% TF = tccmp(TS1,TS2) compares the tc structs tc1 and tc2 and returns logical 1 (true) +% if they are identical, and returns logical 0 (false) otherwise. +% +% See also STRCMP +% +% youkitan 2016-12-02 initial version + +%% check input +mfun = mfilename; + +if ~CheckTC(TC1,mfun) + error('tc1 is not a correctly formed tc.'); +end + +if ~CheckTC(TC2,mfun) + error('tc2 is not a correctly formed tc.'); +end + +%% check overall structure +if ~isequal(size(TC1.tc),size(TC2.tc)) + fprintf('The sizes of the tuning curves are unequal.\n') + pass_flag = 0; + return +end + +if ~isequal(size(TC1.pos_idx),size(TC2.pos_idx)) + fprintf('The sizes of the binned tuning variables are unequal.\n') + pass_flag = 0; + return +end + +if isfield(TC1,'usr') ~= isfield(TC2,'usr') + fprintf('<<< WARNING: Only one of the tc objects has usr fields! >>>\n') +end + +% IMPLEMENT MORE TESTS... + +%% check content +if ~isequal(TC1.tc,TC2.tc) + fprintf('The contens of the tuning curves are not equal.\n') + pass_flag = 0; + return +end + +if ~isequal(TC1.spk_hist,TC2.spk_hist) + fprintf('The contents of the spike counts are not equal.\n') + pass_flag = 0; + return +end + +if ~isequal(TC1.occ_hist,TC2.occ_hist) + fprintf('The contents of the occupancies are not equal.\n') + pass_flag = 0; + return +end + +if ~isequal(TC1.pos_idx,TC2.pos_idx) + fprintf('The contents of the binned tuning variables are not equal.\n') + pass_flag = 0; + return +end + +% IMPLEMENT MORE TESTS... + +%% return true if all checks pass +pass_flag = true; + +end \ No newline at end of file diff --git a/shared/io/datatypes/ts/AddRateTS.m b/shared/io/datatypes/ts/AddRateTS.m new file mode 100644 index 0000000..616924e --- /dev/null +++ b/shared/io/datatypes/ts/AddRateTS.m @@ -0,0 +1,100 @@ +function [TS,gap_idx] = AddRateTS(cfg,TS,tvec) +%ADDRATETS Add usr field to TS that contains the number of occurrances per +% unit time (for S, this is firing rate) +% [TS,gap_idx] = ADDRATETS(cfg,TS,tvec) +% +% INPUTS +% cfg: Struct with field controlling function behavior (see config +% options) +% TS: TS struct containing timestamps such as spike times +% tvec: Time vector (can be obtained from LFP or position data) +% +% OUTPUTS +% TS: TS struct with usr rate field added. +% gap_idx: Indices where gaps were found (is empty if ~cfg.DetectGaps) +% +% CONFIG OPTIONS +% +% cfg.FieldName = 'rate'; +% What to call the usr field that will contain rate data. +% +% cfg.DetectGaps = 0; +% If 1, account for large gaps found in the tvec that might make +% firing rate incorrectly appear lower, if 0 don't. If you don't +% have massive gaps, or lots of them, this will not impact the +% rate very much. +% +% cfg.method = 'mean'; +% How to determine whether there is a gap. Applies if cfg.DetectGaps +% 'raw' - if the gap is greater than cfg.threshold (in tvec +% units) +% 'mean' - if the gap is greater than cfg.threshold times +% above the mean tvec difference +% 'zscore' - if the gap is greater than cfg.threshold standard +% deviations above the mean tvec difference +% +% cfg.threshold = 2; +% Threshold for identifying a gap. Applies only if cfg.DetectGaps +% +% cfg.verbose = 1; +% +% aacarey Jan 2016 + +% Set config defaults +cfg_def.FieldName = 'rate'; +cfg_def.DetectGaps = 0; +cfg_def.method = 'mean'; +cfg_def.threshold = 2; +cfg_def.verbose = 1; +mfun = mfilename; + +cfg = ProcessConfig(cfg_def,cfg,mfun); + +% check inputs +if ~CheckTS(TS) + error('TS data must have been made with the TS constructor. See ts()') +end + +if isempty(TS.t) + error('TS data is empty. Cannot add rate.') +end + +if cfg.verbose; disp([mfun,': Adding rate to usr data']); end + +% Get total time +ElapsedTime = tvec(end)-tvec(1); +gap_idx = []; + +% Handle gaps in time, if requested +if cfg.DetectGaps + TimeDifference = diff(tvec); + + switch cfg.method + case 'raw' + gap_idx = find(TimeDifference >= cfg.threshold); + + case 'mean' + gap_idx = find(TimeDifference >= mean(TimeDifference)*cfg.threshold); + + case 'zscore' + gap_idx = find(zscore(TimeDifference) >= cfg.threshold); + + otherwise + error('Unrecognized option specified in cfg.method') + end + + LostTime = sum(tvec(gap_idx + 1) - tvec(gap_idx)); + ElapsedTime = ElapsedTime - LostTime; +end + +% Calculate rate and add usr field +for iTS = 1:length(TS.t) + nOccurrances = length(TS.t{iTS}); + TS.usr.(cfg.FieldName)(iTS) = nOccurrances/ElapsedTime; +end + +% Record history +TS = History(TS,cfg,mfun); + +end + diff --git a/shared/io/datatypes/ts/CheckTS.m b/shared/io/datatypes/ts/CheckTS.m new file mode 100644 index 0000000..75fff99 --- /dev/null +++ b/shared/io/datatypes/ts/CheckTS.m @@ -0,0 +1,75 @@ +function pass_flag = CheckTS(ts_in,varargin) +%% CHECKTS Check TS for datatype violations +% pass_flag = CheckTS(ts_in,varargin) verifies that input is ts and is well formed. +% +% INPUTS: +% ts_in: ts to be checked +% varargins: +% callername: name of invoking function (for more informative warning messages) +% +% OUTPUTS: +% pass_flag: 1 if all checks pass, 0 if otherwise +% +% Checks performed: +% - .t field must exist (FAIL) +% - all contents of .t must be column vectors (FAIL) +% - is there a .label field? (WARNING) +% - is .label the same length as .t? (FAIL) +% - is there a .type field? (WARNING) +% - does the .type field correctly identify as ts? (WARNING) +% +% see also ts, CheckTSD, CheckIV, CheckTC +% +% aacarey Nov 2015 +% youkitan edit Dec 2016, reformat help, add function name to output + +pass_flag = 1; + +in_mfun = ''; +if ~isempty(varargin) && ischar(varargin{1}) + in_mfun = [' in ',varargin{1}]; +end + +if ~isstruct(ts_in) + pass_flag = 0; + fprintf('FAIL%s by CheckTS: input must be a ts datatype.\n',in_mfun); +else + % check for .t field + if ~isfield(ts_in,'t') + pass_flag = 0; + fprintf('FAIL%s by CheckTS: input ts must contain t field.\n',in_mfun); + + else + % check that all t are column vectors + for iT = 1:length(ts_in.t) + if ~iscolumn(ts_in.t{iT}) + pass_flag = 0; + fprintf('FAIL%s by CheckTS: all contents of ts_in.t must be column vectors.\n',in_mfun); + end + end + end + + % check for label field + if ~isfield(ts_in,'label') + %pass_flag = 0; + fprintf('WARNING%s by CheckTS: ts_in lacks a .label field.\n',in_mfun); + else + % check that label is the same length as t + if length(ts_in.label) ~= length(ts_in.t) + pass_flag = 0; + fprintf('FAIL%s by CheckTS: ts_in.t and ts_in.label must have the same length.\n',in_mfun); + end + end + + % check for type field + if pass_flag && ~isfield(ts_in,'type') + %pass_flag = 0; + fprintf('WARNING%s by CheckTS: input appears to be ts datatype but lacks the .type field.\n',in_mfun); + elseif pass_flag && isfield(ts_in,'type') && ~strcmp(ts_in.type,'ts') + %pass_flag = 0; + fprintf('WARNING%s by CheckTS: input appears to be ts datatype but is identified as %s in the .type field.\n',in_mfun,ts_in.type); + end + +end + +end \ No newline at end of file diff --git a/shared/io/datatypes/ts/OrderSelectS.m b/shared/io/datatypes/ts/OrderSelectS.m new file mode 100644 index 0000000..16b9f7f --- /dev/null +++ b/shared/io/datatypes/ts/OrderSelectS.m @@ -0,0 +1,65 @@ +function S_out = OrderSelectS(cfg_in,S,idx) +%ORDERSELECTS Order and/or select spiketrains +% Return spiketrains and all related information (labels, usr fields) +% according to the specified indices in idx. +% +% S = OrderSelectS(cfg,S,idx) +% +% INPUTS +% S - spiketrains (output from LoadSpikes) +% idx - indices specifying which spiketrains to keep and/or how to +% redorder them +% ex: if idx = [54 5 17] then S_out.t(1) corresponds to S.t(54) +% and so on +% +% OUTPUTS +% S_out - spiketrains ordered or selected according to idx +% +% CONFIG OPTIONS +% cfg.verbose = 1; % if 1, tell me how many spiketrains went in and +% how many went out; if 0, don't +% +% NOTE: this function makes a copy of S and then reorders or selects S.t, +% S.label, and S.usr.data. It also updates cfg history. This function does +% not know about any additional data fields that may be present. +% +% aacarey Oct 2015, edit Nov 2015 + +%% Parse cfg parameters + +cfg_def.verbose = 1; + +mfun = mfilename; % get name of function as string +cfg = ProcessConfig(cfg_def,cfg_in,mfun); + +% do the thing +S_out = S; +S_out.t = S_out.t(idx); +S_out.label = S_out.label(idx); +%S_out.usr.data = S_out.usr.data(idx); + +% also select data from other same-length usr fields +if isfield(S_out,'usr') && ~isempty(S_out.usr) + Sfields = fieldnames(S_out.usr); + for iField = 1:length(Sfields) + S_out.usr.(Sfields{iField}) = S_out.usr.(Sfields{iField})(idx); + end +end + +if cfg.verbose % talk to me + disp([mfun,': ',num2str(length(S.t)),' spiketrains in, ',num2str(length(S_out.t)),' spiketrains out']) +end + +% keep a record of cfg history + +if isfield(S,'cfg') + S_out.cfg = S.cfg; + S_out.cfg.history.mfun = cat(1,S.cfg.history.mfun,mfilename); + S_out.cfg.history.cfg = cat(1,S.cfg.history.cfg,{cfg}); +else + S_out.cfg.history.mfun = mfun; + S_out.cfg.history.cfg = {cfg}; +end + +end + diff --git a/shared/io/datatypes/ts/SelectTS.m b/shared/io/datatypes/ts/SelectTS.m new file mode 100644 index 0000000..6b9f995 --- /dev/null +++ b/shared/io/datatypes/ts/SelectTS.m @@ -0,0 +1,164 @@ +function [ts_out,idx] = SelectTS(cfg_in,ts_in,selectspec) +%SELECTTS Specify ts data to keep. (Also performs reordering). +% This function may be of particular interest for ordering and selecting +% spiketrains (S output from LoadSpikes) but is applicable to any ts data. +% +% ts_out = SELECTTS(cfg,ts_in,selectspec) +% +% INPUTS: +% cfg: config struct with fields controlling function behavior +% ts_in: ts struct +% selectspec: selection specifics, either: +% - [nx1] double: logical array or indices specifying which +% intervals to keep. +% - string: string specifying which usr field to work with. If +% selectspec is a string, the config options cfg.operation, +% cfg.threshold, and cfg.str apply. +% +% OUTPUTS +% ts_out - ts data ordered or selected according to selectspec +% +% CONFIG OPTIONS +% cfg.operation = '>='; How to perform numerical selection, see +% cfg.threshold. +% cfg.threshold = 0; Set a numerical threshold for keeping ts. +% This works on numerical usr contents, but can also be +% applied to strings as long as the first character is +% number-convertible: +% If your field contains strings and cfg.str is +% empty, SelectTS assumes that the first character is a +% number (i.e. a rating) and thresholds based on this +% number. An example would be '1, very good', for which +% SelectTS considers the 1 only. +% cfg.str = ''; If your target usr field contents contain strings that +% are NOT number-convertible, input the string you +% want to select by. If this is not empty, it overrides +% numerical selection. Examples of non-number-convertible +% strings might be 'good' or 'maybe' or 'poor'. +% cfg.verbose = 1; % if 1, tell me how many ts went in and how many +% went out; if 0, don't +% +% aacarey Oct 2015, edit Nov 2015 + +%% + + +mfun = mfilename; +if ~CheckTS(ts_in) + error('ts_in is either not a ts datatype or is poorly formed.') +end + +% Parse cfg parameters +cfg_def.operation = '>='; +cfg_def.threshold = 0; +cfg_def.str = ''; % if this is not empty, it overrides numerical selection +cfg_def.verbose = 1; + +cfg = ProcessConfig(cfg_def,cfg_in,mfun); + +% choose which thing to do +if islogical(selectspec) || isnumeric(selectspec) + goto = 'there'; % -_- +elseif ischar(selectspec) + goto = 'here'; % D: +else + error('selectspec must be a logical array, numeric array of indices, or a string specifying a usr field name.') +end + +switch goto + case 'here' + % make sure usr exists + if ~isfield(ts_in,'usr') + error([mfun,': ts_in requires usr for this type of selection.']) + end + % check that the field actually exists and that it's the right length + if ischar(selectspec) && ~isfield(ts_in.usr,selectspec) + error([selectspec,' does not exist.']) + elseif ischar(selectspec) && length(ts_in.usr.(selectspec)) ~= length(ts_in.t) + error(['ts_in.usr.',selectspec,' must have the same dimensions as ts_in.t.']) + end + + % if the field contains strings, get ratings in numerical form + if isempty(cfg.str) && ~isnumeric(ts_in.usr.(selectspec)(1)) + type = 1; + temp = nan(size(ts_in.usr.(selectspec))); + for ii = 1:length(temp) + temp(ii,1) = str2double(ts_in.usr.(selectspec){ii,1}(1)); % we assume the rating is the first character in the string + end + elseif isempty(cfg.str) && isnumeric(ts_in.usr.(selectspec)) + type = 1; + temp = ts_in.usr.(selectspec); + elseif ~isempty(cfg.str) + type = 2; + temp = ts_in.usr.(selectspec); + end + + assignin('base','temp',temp) + + % do the thing + switch type + case 1 + switch cfg.operation + case '>' + keep = temp > cfg.threshold; + case '>=' + keep = temp >= cfg.threshold; + case '<' + keep = temp < cfg.threshold; + case '<=' + keep = temp <= cfg.threshold; + case '=' + keep = temp == cfg.threshold; + otherwise + error('Unrecognized cfg.operation') + end + case 2 + keep = nan(size(temp)); + for iStr = 1:length(temp) + keep(iStr) = strcmp(cfg.str,temp(iStr)); + end + end + + keep = logical(keep); + + case 'there' + keep = selectspec; +end + +% select +ts_out = ts_in; +ts_out.t = ts_out.t(keep); +ts_out.label = ts_out.label(keep); + +% also select data from other same-length usr fields +if isfield(ts_out,'usr') && ~isempty(ts_out.usr) + TSfields = fieldnames(ts_out.usr); + for iField = 1:length(TSfields) + ts_out.usr.(TSfields{iField}) = ts_out.usr.(TSfields{iField})(keep); + end +end + +% make idx output +if islogical(keep) + idx = find(keep); +elseif isnumeric(keep) + idx = keep; +end + +if cfg.verbose % talk to me + disp([mfun,': ',num2str(length(ts_in.t)),' ts in, ',num2str(length(ts_out.t)),' ts out']) +end + +% keep a record of cfg history + +if isfield(ts_in,'cfg') + ts_out.cfg = ts_in.cfg; + ts_out.cfg.history.mfun = cat(1,ts_in.cfg.history.mfun,mfilename); + ts_out.cfg.history.cfg = cat(1,ts_in.cfg.history.cfg,{cfg}); +else + ts_out.cfg.history.mfun = mfun; + ts_out.cfg.history.cfg = {cfg}; +end + +end + diff --git a/shared/io/datatypes/ts/ShuffleTS.m b/shared/io/datatypes/ts/ShuffleTS.m new file mode 100644 index 0000000..f5aeab4 --- /dev/null +++ b/shared/io/datatypes/ts/ShuffleTS.m @@ -0,0 +1,58 @@ +function ts_in = ShuffleTS(cfg_in,ts_in) +% function ts_out = ShuffleTS(cfg,ts_in) +% +% shuffle spikes in ts struct +% +% INPUTS: +% +% ts_in: input ts +% +% OUTPUT: +% +% ts_out: shuffled ts +% +% CONFIGS: +% +% cfg_def.mode = 1; % 1: shuffle labels, 2: shuffle ISIs, 3: uniform random +% with same number of spikes +% cfg_def.t0 = 0; % should specify this for mode 2, 3 to indicate start of +% window +% cfg.def.t1 = 0; % should specify for mode 2, 3 to indicate end of window +% +% MvdM 2015-02-09 initial version + +cfg_def.mode = 1; % 1: shuffle labels, 2: shuffle ISIs + +cfg = ProcessConfig2(cfg_def,cfg_in); + +switch cfg.mode + + case 1 % shuffle labels + + r = randperm(length(ts_in.t)); + ts_in.t = ts_in.t(r); + + case 2 % shuffle ISIs + + ts_in = restrict2(ts_in,cfg.t0,cfg.t1); + + for iC = 1:length(ts_in.t) + + isi = diff(cat(1,cfg.t0,ts_in.t{iC})); + r = randperm(length(isi)); + ts_in.t{iC} = cfg.t0+cumsum(isi(r)); + + end + + case 3 % random uniform + + ts_in = restrict2(ts_in,cfg.t0,cfg.t1); + + for iC = 1:length(ts_in.t) + + N = length(ts_in.t{iC}); + + ts_in.t{iC} = cfg.t0 + (cfg.t1-cfg.t0)*rand(N,1); + + end +end \ No newline at end of file diff --git a/shared/io/datatypes/ts/concatenateTS.m b/shared/io/datatypes/ts/concatenateTS.m new file mode 100644 index 0000000..5f6f9d7 --- /dev/null +++ b/shared/io/datatypes/ts/concatenateTS.m @@ -0,0 +1,65 @@ +function ts_out = concatenateTS(ts1,ts2) +%concatenateTS concatenate two ts objects +% Concatenates .t fields and all related information (labels, usr fields) +% +% ts.out = concatenateTS(ts1,ts2) +% +% INPUTS +% ts1 - first ts (second ts will be appended to this) +% ts2 - second ts +% +% OUTPUTS +% ts_out - concatenated ts +% +% MvdM 2015-11-18 initial version + +mfun = mfilename; % get name of function as string + +% append t and label +ts_out = ts1; + +ts_out.t = cat(2,ts_out.t,ts2.t); +ts_out.label = cat(2,ts_out.label,ts2.label); + +% prepare to append data from usr fields + +% first check if have same usr fields +if ~isfield(ts_out,'usr') + nFields1 = 0; +else + ts1_fieldnames = fieldnames(ts_out.usr); + nFields1 = length(ts1_fieldnames); +end + +if ~isfield(ts2,'usr') + nFields2 = 0; +else + ts2_fieldnames = fieldnames(ts2.usr); + nFields2 = length(ts2_fieldnames); +end + +if nFields1 ~= nFields2 + error('Inputs have non-matching usr field counts.'); +elseif sum(strcmp(sort(ts1_fieldnames),sort(ts2_fieldnames))) ~= nFields1 % some name mismatch + error('Inputs have non-matching usr field names.'); +end + +% append +if isfield(ts_out,'usr') && ~isempty(ts_out.usr) + ts1_fieldnames = fieldnames(ts_out.usr); + for iField = 1:length(ts1_fieldnames) + ts_out.usr.(ts1_fieldnames{iField}) = cat(2,ts_out.usr.(ts1_fieldnames{iField}),ts2.usr.(ts1_fieldnames{iField})); + end +end + +% keep a record of cfg history -- note this is a case with two parents, not +% handled well +if isfield(ts_out,'cfg') + ts_out.cfg.history.mfun = cat(1,ts_out.cfg.history.mfun,mfilename); + ts_out.cfg.history.cfg = cat(1,ts_out.cfg.history.cfg,{}); +else + ts_out.cfg.history.mfun = mfun; + ts_out.cfg.history.cfg = {}; +end + + diff --git a/shared/io/datatypes/ts/firstSpike.m b/shared/io/datatypes/ts/firstSpike.m new file mode 100644 index 0000000..fd14d58 --- /dev/null +++ b/shared/io/datatypes/ts/firstSpike.m @@ -0,0 +1,7 @@ +function spk_t = firstSpike(S) +% function spk_t = firstSpike(S) +% +% returns time of first spike + +spk = vertcat(S.t{:}); +spk_t = min(spk); \ No newline at end of file diff --git a/shared/io/datatypes/ts/getISI.m b/shared/io/datatypes/ts/getISI.m new file mode 100644 index 0000000..a509607 --- /dev/null +++ b/shared/io/datatypes/ts/getISI.m @@ -0,0 +1,55 @@ +function isi = getISI(cfg_in,S) +%% GETISI Calculate interspike intervals +% isi = getISI(cfg_in,S) returns the ISI values for each cell in a TS variable +% +% INPUTS: +% cfg_in: config field +% S: ts struct +% +% OUTPUTS: +% isi: {NxI} cell array where N is the number of cells and I is the number of +% intervals set by the cfg.iv parameter (default 1). Each cell contains the +% interspike (event) intervals as a [Dx1] double. +% +% CONFIG OPTIONS +% cfg.iv = []; If specified, return ISIs for each interval (output is nCells x +% nIntervals). Otherwise the interval is the whole TS time span. +% +% See also GETSPIKECOUNT +% +% youkitan 2016-07-28 initial version + +%% Parse cfg parameters and error check + +cfg_def = []; +cfg_def.iv = []; % if specified, return spike counts for each interval (output is nCells x nIntervals) + +cfg = ProcessConfig(cfg_def,cfg_in); + +if ~CheckTS(S) + error('Input is not a correctly formed ts.'); +end + +%% do the things + +% create dummy iv to avoid repetition later +if isempty(cfg.iv) + cfg.iv = iv(); + cfg.iv.tstart = firstSpike(S)-eps; + cfg.iv.tend = lastSpike(S)+eps; +end + +% create array for efficiency +nCells = length(S.t); +nIV = length(cfg.iv.tstart); +isi = cell(nCells,nIV); + +% iterate and get ISIs +for iI = 1:nIV % for each interval, get spike counts + this_S = restrict(S,cfg.iv.tstart(iI),cfg.iv.tend(iI)); + + for iC = 1:nCells + isi{iC,iI} = diff(this_S.t{iC}); + end % of cells + +end % of intervals diff --git a/shared/io/datatypes/ts/getSpikeCount.m b/shared/io/datatypes/ts/getSpikeCount.m new file mode 100644 index 0000000..d2c4509 --- /dev/null +++ b/shared/io/datatypes/ts/getSpikeCount.m @@ -0,0 +1,42 @@ +function spk_c = getSpikeCount(cfg_in,S) +% function spk_c = getSpikeCount(cfg_in,S) +% +% returns spike counts (number of spikes) in a nCells x 1 vector of S input +% +% if optional cfg_in.iv argument is specified, return spike counts for each +% interval (output is nCells x nIntervals) +% +% MvdM 2016-01-11 initial version + +cfg_def = []; +cfg_def.iv = []; % if specified, return spike counts for each interval (output is nCells x nIntervals) +cfg_def.convertToRate = 0; % if 1, divide spike count by length of interval +cfg_def.verbose = 1; + +cfg = ProcessConfig(cfg_def,cfg_in); + +if ~CheckTS(S) + error('Input is not a correctly formed ts.'); +end + +if isempty(cfg.iv) % create dummy iv to avoid repetition later + cfg.iv = iv(); + cfg.iv.tstart = firstSpike(S)-eps; + cfg.iv.tend = lastSpike(S)+eps; +end + +nCells = length(S.t); + +for iI = length(cfg.iv.tstart):-1:1 % for each interval, get spike counts + + this_S = restrict(S,cfg.iv.tstart(iI),cfg.iv.tend(iI)); + + for iC = nCells:-1:1 + + spk_c(iC,iI) = length(this_S.t{iC}); + + if cfg.convertToRate, spk_c(iC,iI) = spk_c(iC,iI) ./ (cfg.iv.tend(iI) - cfg.iv.tstart(iI)); end + + end % of cells + +end % of intervals \ No newline at end of file diff --git a/shared/io/datatypes/ts/lastSpike.m b/shared/io/datatypes/ts/lastSpike.m new file mode 100644 index 0000000..20acc38 --- /dev/null +++ b/shared/io/datatypes/ts/lastSpike.m @@ -0,0 +1,7 @@ +function spk_t = lastSpike(S) +% function spk_t = lastSpike(S) +% +% returns time of last spike + +spk = vertcat(S.t{:}); +spk_t = max(spk); \ No newline at end of file diff --git a/shared/io/datatypes/ts/removeEmptyCells.m b/shared/io/datatypes/ts/removeEmptyCells.m new file mode 100644 index 0000000..bca3163 --- /dev/null +++ b/shared/io/datatypes/ts/removeEmptyCells.m @@ -0,0 +1,13 @@ +function [S,keep_idx] = removeEmptyCells(S) +% function [S,keep_idx] = removeEmptyCells(S) + +cfg = []; +mfun = mfilename; + +keep_idx = ~cellfun(@isempty,S.t); +S.t = S.t(keep_idx); +S.label = S.label(keep_idx); + +% housekeeping +S.cfg.history.mfun = cat(1,S.cfg.history.mfun,mfun); +S.cfg.history.cfg = cat(1,S.cfg.history.cfg,{cfg}); diff --git a/shared/io/datatypes/ts/ts.m b/shared/io/datatypes/ts/ts.m new file mode 100644 index 0000000..2ad2960 --- /dev/null +++ b/shared/io/datatypes/ts/ts.m @@ -0,0 +1,71 @@ +function ts_out = ts(varargin) +% TS Timestamp datatype constructor +% The ts datatype differs from tsd in that it contains only the timestamps of events +% with no additional accompanying data. +% +% function ts_out = TS(t,label) +% +% function ts_out = TS(t) +% +% function ts_out = TS(N) returns an ts struct with N specified empty cells +% +% OUTPUTS: +% ts_out: [1x1] ts struct with fields: +% .type - 'ts'; datatype identification +% .t - {1xn} cell, each containing a [nTimestamps x 1] double +% .label - {1xn} cell containing unique identifiers for the +% timestamps contained within the .t field, such as a +% record of the source (in the case of S, the label +% corresponds to the tetrode that recorded the spike +% trains) +% .cfg - record of the data's history including config +% parameters and functions visited +% +% see also tsd, CheckTS +% +% MvdM 2014-06-17 +% aacarey edit Nov 2015 +% youkitan edit 2016-11-22 added initialization with size +% youkitan edit Jan 2017 added intitialization with data + +ts_out.type = 'ts'; +ts_out.t = {}; +ts_out.label = {}; + +% initialization of a ts struct with specified data +if nargin == 2 + timestamps = varargin{1}; + labels = varargin{2}; + + if ~iscell(timestamps) + error('Input .t data are not in cell format.') + elseif ~iscell(labels) + error('Input .label data are not in cell format.') + end + +% initialization of an empty ts struct with prespecified number of data bins +elseif nargin == 1 + if isnumeric(varargin{1}) && size(varargin{1},1) == 1 && size(varargin{1},2) == 1 + ts_out.t = cell(1,varargin{1}); + ts_out.label = cell(1,varargin{1}); + + for iC = 1:varargin{1} + ts_out.t{iC} = zeros(0,1); + ts_out.label{iC} = iC; + end + elseif iscell(varargin{1}) + ts_out.t = varargin{1}; + else + error('Input is neither a numeric or cell type.') + end + +% initialization of an empty ts struct +elseif nargin == 0 + %doesn't have to do aything +else + error('Invalid number of arguments.') +end + +% housekeeping +ts_out.cfg.history.mfun{1} = mfilename; +ts_out.cfg.history.cfg{1} = []; \ No newline at end of file diff --git a/shared/io/datatypes/ts/tscmp.m b/shared/io/datatypes/ts/tscmp.m new file mode 100644 index 0000000..0d2d554 --- /dev/null +++ b/shared/io/datatypes/ts/tscmp.m @@ -0,0 +1,49 @@ +function pass_flag = tscmp(TS1,TS2) +%% tscmp Compare TS objects +% TF = strcmp(TS1,TS2) compares the ts structs ts1 and ts2 and returns logical 1 (true) +% if they are identical, and returns logical 0 (false) otherwise. +% +% See also STRCMP +% +% youkitan 2016-03-22 initial version +% youkitan 2016-08-05 input check + +%% check input +if ~CheckTS(TS1) + error('ts1 is not a correctly formed ts.'); +end + +if ~CheckTS(TS2) + error('ts2 is not a correctly formed ts.'); +end + +%% check overall structure +if length(TS1.t) ~= length(TS2.t); + fprintf('There are unequal numbers of cells') + pass_flag = false; + return; +end + +if length(TS1.label) ~= length(TS2.label); + fprintf('There are unequal numbers of labels') + pass_flag = false; + return; +end + +%% check content + +if any(~cellfun(@isequal,TS1.t,TS2.t)); + fprintf('The spike times are not the same across all cells... \n') + + idx = find(~cellfun(@isequal,TS1.t,TS2.t)); + fprintf(['The following cells are not the same between files: ',num2str(idx),'\n']) + + pass_flag = false; + return; +end + + +%% return true if all checks pass +pass_flag = true; + +end \ No newline at end of file diff --git a/shared/io/datatypes/tsd/CheckTSD.m b/shared/io/datatypes/tsd/CheckTSD.m new file mode 100644 index 0000000..de3ebf8 --- /dev/null +++ b/shared/io/datatypes/tsd/CheckTSD.m @@ -0,0 +1,60 @@ +function pass_flag = CheckTSD(tsd_in,varargin) +%% CHECKTSD Check TSD for datatype violations +% pass_flag = CheckTSD(tsd_in,varargin) verifies that input is TSD and well formed. +% +% INPUTS: +% tsd_in: tsd to be checked +% varargins: +% callername: name of invoking function (for more informative warning messages) +% +% OUTPUTS: +% pass_flag: 1 if all checks pass, 0 if otherwise +% +% Checks performed: +% - .units field exists (WARNING) +% - .data field cannot be a cell (FAIL) +% - is the number of signals greater than the number of samples? (WARNING) +% - number of samples in .data and .tvec must be equal (FAIL) +% +% see also tsd, CheckTS, CheckIV, CheckTC +% +% MvdM 2014-06-24 +% youkitan edit Dec 2016, reformat help, add function name to output +% youkitan edit Feb 2017, included check for units + +pass_flag = 1; + +in_mfun = ''; +if ~isempty(varargin) && ischar(varargin{1}) + in_mfun = [' in ',varargin{1}]; +end + +if ~isstruct(tsd_in) + pass_flag = 0; + fprintf('FAIL%s by CheckTSD: data is not a struct.\n',in_mfun); + return; +end + +if ~isfield(tsd_in,'units') + fprintf('WARNING%s by CheckTSD: Missing units. TSD was generated with an outdated version of the constructor.\n',in_mfun); +end + +if iscell(tsd_in.data) + pass_flag = 0; + fprintf('FAIL%s by CheckTSD: data is a cell.\n',in_mfun); + return; +end + +nSignals = size(tsd_in.data,1); +nSamples_data = size(tsd_in.data,2); + +if nSignals > nSamples_data + %pass_flag = 0; + fprintf('WARNING%s by CheckTSD: more signals (%d) than samples (%d) in data.\n',in_mfun,nSignals,nSamples_data); +end + +nSamples_tvec = length(tsd_in.tvec); +if nSamples_data ~= nSamples_tvec + pass_flag = 0; + fprintf('FAIL in %s by CheckTSD: samples in data (%d) does not match samples in tvec.\n',in_mfun,nSamples_data,nSamples_tvec); +end \ No newline at end of file diff --git a/shared/io/datatypes/tsd/MergeTSD.m b/shared/io/datatypes/tsd/MergeTSD.m new file mode 100644 index 0000000..c393708 --- /dev/null +++ b/shared/io/datatypes/tsd/MergeTSD.m @@ -0,0 +1,102 @@ +function tsd_out = MergeTSD(cfg_in,varargin) +%MERGETSD Merge data from multiple TSDs into a single, new TSD +% Corresponding data values from each TSD.data are combined according to +% cfg.method. So, tsd1.data(1) is combined with tsd2.data(1) and so on. +% If data values in one TSD are much larger than in another, you may +% consider rescaling your TSDs beforehand so that one is not weighted +% higher than another: see rescale() and rescmean(). +% +% tsd_out = MergeTSD(cfg,tsd1,tsd2,...) +% +% cfg.method = 'mean'; +% 'median' - the mean of the middle two data values in sorted order +% 'mean' - the average of the data values +% 'geometricmean' - the nth root of the mean of the data values +% note: cfg.geometricmean is not compatible with TSDs that contain +% negative values if the nth root is even +% 'addition' - not implemented +% 'subtraction' - not implemented +% +% cfg.verbose = 1; % 1 - tell me what you are doing. 0 - don't +% +% Why use MergeTSD? +% 1. To combine score vectors from different detectors +% 2. In some cases, you may want to merge your CSCs to get a general LFP. +% +% aacarey Sept 2015 + +cfg_def.method = 'mean'; +cfg_def.verbose = 1; + +mfun = mfilename; +cfg = ProcessConfig(cfg_def,cfg_in,mfun); + +numArgs = length(varargin); +if numArgs < 2 + error('Require at least two input TSDs') +end + +% check that the inputs are TSDs +for iArg = 1:numArgs + pass = CheckTSD(varargin{iArg}); + + if ~ pass + error([inputname(iArg+1),' must be a TSD object.']) + end +end + +% check that the inputs all have the same number of samples +checktvecs = nan(1,numArgs); +for iArg = 1:numArgs + checktvecs(iArg) = length(varargin{iArg}.tvec); +end +checktvecs = unique(checktvecs); +if length(checktvecs) > 1 + error('All TSDs must have the same number of samples.') +end + +% tell me what you are doing +if cfg.verbose + names = []; + for iArg = 1:numArgs + if iArg < numArgs + names = [names,' ',inputname(iArg+1)]; + else + names = [names,' and ', inputname(iArg+1)]; + end + end + disp(['MergeTSD: outputting the ',cfg.method,' of',names]) +end + +% time to do the thing + +% collect all the data fields into an array +allTSDs = nan(numArgs,length(varargin{iArg}.tvec)); +for iArg = 1:numArgs + allTSDs(iArg,:) = varargin{iArg}.data; +end + +% merge the data, as specified by user +switch cfg.method + case 'median' + data = median(allTSDs,1); + case 'mean' + data = mean(allTSDs,1); + case 'geometricmean' + data = nthroot(prod(allTSDs,1),numArgs); + case 'sum' + data = sum(allTSDs,1); + otherwise + error('Unrecognized cfg.method. Better check that spelling ^_^') +end + +% form the output TSD + +tsd_out = tsd(varargin{1}.tvec,data); + +% need to add label if csc??? + +tsd_out = History(tsd_out,mfun,cfg); + +end + diff --git a/shared/io/datatypes/tsd/TSD_SelectChannel.m b/shared/io/datatypes/tsd/TSD_SelectChannel.m new file mode 100644 index 0000000..28cea51 --- /dev/null +++ b/shared/io/datatypes/tsd/TSD_SelectChannel.m @@ -0,0 +1,25 @@ +function tsd_out = TSD_SelectChannel(tsd_in,channel_label) +% function tsd_out = TSD_SelectChannel(tsd_in,channel_label) +% +% returns tsd containing only channels with labels in channel_label +% +% MvdM 2015-11-03 initial version + +if ~iscell(channel_label) + keep_idx = strmatch(channel_label,tsd_in.label,'exact'); +else + keep_idx = []; + for iCh = 1:length(channel_label) + + keep_idx = cat(1,keep_idx,strmatch(channel_label{iCh},tsd_in.label,'exact')); + + end +end + +if isempty(keep_idx) + error('No channels selected.'); +end + +tsd_out = tsd_in; +tsd_out.data = tsd_out.data(keep_idx,:); +tsd_out.label = tsd_out.label(keep_idx); \ No newline at end of file diff --git a/shared/io/datatypes/tsd/TSD_getidx.m b/shared/io/datatypes/tsd/TSD_getidx.m new file mode 100644 index 0000000..998cdf6 --- /dev/null +++ b/shared/io/datatypes/tsd/TSD_getidx.m @@ -0,0 +1,50 @@ +function idx = TSD_getidx(in,varargin) +% function idx = TSD_getidx(in,varargin) +% +% returns idxs into ts or tsd, delimited by input interval(s) +% +% if a data time is equal to an interval edge, this data point is INCLUDED +% +% usage: +% idx = TSD_getidx(in,iv) +% idx = TSD_getidx(in,tstart,tend) % tstart and tend can be vectors +% +% MvdM 2015-04-26 initial version: SLOW for large numbers of inputs, fix + +% convert input arguments to iv if not already done +if nargin == 2 + if ~strcmp(varargin{1}.cfg.history.mfun{1},'iv') + error('Single input argument must be iv (interval) type.'); + end + + iv_use = varargin{1}; + +elseif nargin == 3 + + iv_use = iv(varargin{1},varargin{2}); + +else + error('Unsupported number of input arguments.'); +end + +% get indices to keep +if isfield(in,'tvec') % tsd + keep = false(size(in.tvec)); + + for iT = 1:length(iv_use.tstart) + keep = keep | (in.tvec >= iv_use.tstart(iT) & in.tvec <= iv_use.tend(iT)); + end + +elseif isfield(in,'t') % ts + + for iC = length(in.t):-1:1 + keep{iC} = false(size(in.t{iC})); + + for iT = 1:length(iv_use.tstart) + keep{iC} = keep{iC} | (in.t{iC} >= iv_use.tstart(iT) & in.t{iC} <= iv_use.tend(iT)); + end + end + +end + +idx = find(keep); \ No newline at end of file diff --git a/shared/io/datatypes/tsd/TSD_getidx2.m b/shared/io/datatypes/tsd/TSD_getidx2.m new file mode 100644 index 0000000..946c85b --- /dev/null +++ b/shared/io/datatypes/tsd/TSD_getidx2.m @@ -0,0 +1,81 @@ +function idx = TSD_getidx(in,varargin) +% function idx = TSD_getidx(in,varargin) +% +% returns idxs into ts or tsd, delimited by input interval(s) +% +% if a data time is equal to an interval edge, this data point is INCLUDED +% +% usage: +% idx = TSD_getidx(in,iv) +% idx = TSD_getidx(in,tstart,tend) % tstart and tend can be vectors +% +% MvdM 2015-04-26 initial version: SLOW for large numbers of inputs, fix +% MvdM 2016-01-07 faster version based on nearest_idx3 (taken from +% restrict2()) + +% convert input arguments to iv if not already done +if nargin == 2 + if ~strcmp(varargin{1}.cfg.history.mfun{1},'iv') + error('Single input argument must be iv (interval) type.'); + end + + iv_use = varargin{1}; + +elseif nargin == 3 + + iv_use = iv(varargin{1},varargin{2}); + +else + error('Unsupported number of input arguments.'); +end + +% get indices to keep +if isfield(in,'tvec') % tsd + + keep = []; + + for iT = 1:length(iv_use.tstart) + + ind_start = nearest_idx3(iv_use.tstart(iT),in.tvec,1); + ind_end = nearest_idx3(iv_use.tend(iT),in.tvec,-1); + keep = [keep ind_start:ind_end]; + + end + keep = unique(keep); + +elseif isfield(in,'t') % ts + + for iC = length(in.t):-1:1 + keep{iC} = []; + + lookups = in.t{iC}; + if isempty(lookups) + continue + end + + for iT = 1:length(iv_use.tstart) %for each input interval... + + iv_start = iv_use.tstart(iT); + iv_end = iv_use.tend(iT); + + ind_start = nearest_idx3(iv_start,lookups,1); + ind_end = nearest_idx3(iv_end,lookups,-1); + + keepvals_new = ind_start:ind_end; + + if (ind_start == 1 && ind_end == 1) ||... + (ind_start == length(lookups) && ind_end == length(lookups)) + keepvals_new = find(in.t{iC} >= iv_use.tstart(iT) & in.t{iC} <= iv_use.tend(iT))'; + end + + if ind_start <= ind_end + keep{iC} = [keep{iC} keepvals_new]; + end + + end % iterate events + + end % iterate cells + +end + +idx = keep; \ No newline at end of file diff --git a/shared/io/datatypes/tsd/TSDbySpace.m b/shared/io/datatypes/tsd/TSDbySpace.m new file mode 100644 index 0000000..74995a9 --- /dev/null +++ b/shared/io/datatypes/tsd/TSDbySpace.m @@ -0,0 +1,38 @@ +function space_av = TSDbySpace(cfg_in,pos,tsd_in) +% function space_av = TSDbySpace(cfg,pos,tsd_in) +% +% construct 2D matrix with spatial averages of the values in tsd_in +% +% inputs: +% +% outputs: +% +% cfg options: +% +% + +cfg_def = []; +cfg_def.x_edges = []; +cfg_def.y_edges = []; + +cfg = ProcessConfig(cfg_def,cfg_in); + +pos_mat(:,1) = getd(pos,'y'); +pos_mat(:,2) = getd(pos,'x'); + +[occ_hist,edges,mid,pos_idx] = histcn(pos_mat,cfg.y_edges,cfg.x_edges); + +tsd_pos = interp1(tsd_in.tvec,tsd_in.data,pos.tvec,'nearest'); % tsd value for each pos sample + +space_av = nan(length(cfg.y_edges)-1,length(cfg.x_edges)-1); % need to fix points that fall on last edge though + +temp = sub2ind(size(space_av),pos_idx(:,1),pos_idx(:,2)); % idx into binned space for each pos sample + +idx_list = unique(temp); +for iI = 1:length(idx_list) + curr_idx = idx_list(iI); + space_av(curr_idx) = nanmean(tsd_pos(find(temp == curr_idx))); +end + + + \ No newline at end of file diff --git a/shared/io/datatypes/tsd/TSDbyZ.m b/shared/io/datatypes/tsd/TSDbyZ.m new file mode 100644 index 0000000..fc0e0db --- /dev/null +++ b/shared/io/datatypes/tsd/TSDbyZ.m @@ -0,0 +1,33 @@ +function z_av = TSDbyZ(cfg_in,z,tsd_in) +% function z_av = TSDbyZ(cfg,z,tsd_in) +% +% compute averages of the values in tsd_in binned by z's cfg.edges +% +% inputs: +% +% outputs: +% +% cfg options: +% +% + +cfg_def = []; +cfg_def.edges = []; + +cfg = ProcessConfig(cfg_def,cfg_in); + +[z_hist,z_idx] = histc(z.data(1,:),cfg.edges); % need to fix points that fall on last edge +[~,z_idx] = trim_histc(z_hist,z_idx); + +tsd_z = interp1(tsd_in.tvec,tsd_in.data(1,:),z.tvec,'nearest'); % tsd value for each z sample + +z_av = nan(length(cfg.edges)-1,1); + +idx_list = unique(z_idx); +for iI = 1:length(idx_list) + curr_idx = idx_list(iI); + z_av(curr_idx) = nanmean(tsd_z(z_idx == curr_idx)); +end + + + \ No newline at end of file diff --git a/shared/io/datatypes/tsd/TSDpeth.m b/shared/io/datatypes/tsd/TSDpeth.m new file mode 100644 index 0000000..a4aa6d8 --- /dev/null +++ b/shared/io/datatypes/tsd/TSDpeth.m @@ -0,0 +1,129 @@ +function out = tsdPETH(cfg_in,tsd_in,t_in) +% function out = tsdPETH(cfg_in,tsd_in,t_in) +% +% compute peri-event time histogram (average) +% +% INPUTS: +% +% tsd_in: tsd to compute average over +% t_in: times to average: can be raw timestamps, ts or iv data. If iv, cfg.window is ignored +% +% OUTPUTS: +% +% out: tsd with PETH +% +% cfg options: +% +% cfg_def.window = [-2 2]; % start and end times of window (in s) +% cfg_def.dt = []; % time step, specify this for 'interp' mode +% cfg_def.mode = 'raw'; % 'raw' or 'interp'; if 'interp', need to specify cfg_def.dt +% cfg_def.interp_mode = 'linear'; +% +% EXAMPLE USAGE: +% +% cd('R042-2013-08-18'); LoadExpKeys; LoadMetadata; +% cfg = []; cfg.fc = ExpKeys.goodSWR(1); +% lfp = LoadCSC(cfg); +% peth = TSDpeth([],lfp,metadata.taskvars.trial_iv.tstart(2:end)); +% plot(peth); +% +% NOTES: +% +% 'raw' mode collects tsd samples that fall within the specified intervals. +% Because it does not interpolate, it can easily fail if there are gaps in +% the data or the diffs between samples are unequal. +% +% 'interp' interpolates on a fixed timebase, so guarantees the number of +% samples will be constant. +% +% multidimensional TSDs are not yet supported! +% +% MvdM 2017-08-16 initial version + +cfg_def = []; +cfg_def.window = [-2 2]; % start and end times of window (in s) +cfg_def.dt = []; % time step, specify this for 'interp' mode +cfg_def.mode = 'raw'; % 'raw' or 'interp'; if 'interp', need to specify cfg_def.dt +cfg_def.interp_mode = 'linear'; + +cfg = ProcessConfig2(cfg_def,cfg_in); + +if ~CheckTSD(tsd_in) + error('Incorrectly formed tsd input.'); +end + +if size(tsd_in.data,1) > 1 + error('More than 1 dimension in TSD input not yet supported.'); +end + +if ~isfield(t_in,'type') % assume raw times + this_iv = iv(t_in + cfg.window(1), t_in + cfg.window(2)); +else + switch t_in.type + case 'iv' + + if ~CheckIV(t_in) % this should check that there are no t_end before t_start, etc... + error('Incorrectly formed iv input.'); + end + + this_iv = t_in; + + case ts + + if length(t_in.t) ~= 1 + error('ts input must have exactly one .t cell'); + end + + this_iv = iv(t_in.t{1} + cfg.window(1), t_in.t{1} + cfg_window(2)); + + otherwise + + error('t input must be raw times, iv or ts data') + + end +end + +nT = length(this_iv.tstart); + +switch cfg.mode + + case 'raw' + + % make big idx matrix + start_idx = nearest_idx3(this_iv.tstart,tsd_in.tvec); + end_idx = nearest_idx3(this_iv.tend,tsd_in.tvec); + + if length(unique(end_idx-start_idx)) ~= 1 % unequal length trials + error('Raw mode requires equal tsd samples for each trial.'); + end + + for iT = nT:-1:1 % slow! could be vectorized + + out_data(iT,:) = tsd_in.data(start_idx(iT):end_idx(iT)); % need to generalize to 2-D + + end + out_tvec = tsd_in.tvec(start_idx(1):end_idx(1)); + out_tvec = out_tvec-nanmean(out_tvec); % this is an approximation -- depends on exact spacing of input tsd + + case 'interp' + + if isempty(cfg.dt) + error('interp mode requires cfg.dt to be specified.'); + end + + for iT = nT:-1:1 % slow! should be vectorized to do one interp1 first, then reshape + + out_data(iT,:) = interp1(tsd_in.tvec,tsd_in.data,this_iv.tstart(iT):cfg.dt:this_iv.tend(iT),cfg.interp_mode); % need to generalize to 2-D + + end + out_tvec = cfg.window(1):cfg.dt:cfg.window(2); + + otherwise + + error('unknown cfg.mode %s',cfg.mode); +end + +% average and package +out = tsd; +out.data = nanmean(out_data); +out.tvec = out_tvec; \ No newline at end of file diff --git a/shared/io/datatypes/tsd/TSDtoFT.m b/shared/io/datatypes/tsd/TSDtoFT.m new file mode 100644 index 0000000..eeb7377 --- /dev/null +++ b/shared/io/datatypes/tsd/TSDtoFT.m @@ -0,0 +1,53 @@ +function data_ft = TSDtoFT(cfg_in,data) +% function data_ft = TSDtoFT(cfg_in,data_tsd) +% +% converts tsd into fieldtrip (ft) data structure +% +% cfg_def.mode = 'as-is'; % {'as-is','resample'}, defines how to deal with gaps in data +% +% MvdM 2014-11-12 initial version +% NOTE multiple channels not yet implemented! + +cfg_def = []; +cfg_def.mode = 'as-is'; % {'as-is','resample'}, defines how to deal with gaps in data +cfg = ProcessConfig(cfg_def,cfg_in); + +if ~CheckTSD(data) + return; +end + +dts = unique(diff(data.tvec)); +if length(dts) > 1 + fprintf('\nTSDtoFT.m: WARNING: tvec diffs are not constant, cannot determine Fs.'); + % could approximate with a median if matches average closely enough + Fs = 1./median(dts); + fprintf('\nTSDtoFT.m: Fs %.2f estimated.\n',Fs); +else + Fs = 1./dts; + fprintf('\nTSDtoFT.m: Fs %.2f detected.\n',Fs); +end + +% +data_ft = []; +switch cfg.mode + + case 'as-is' + nSamples = size(data.data,2); + + data_ft.trial{1} = data.data; + data_ft.time{1} = data.tvec; + + case 'resample' % currently does interpolation of data, could be improved with options like inserting NaNs or zeros if no sample nearby + + data_ft.time{1} = data.tvec(1):1./Fs:data.tvec(end); + data_ft.trial{1} = interp1(data.tvec,data.data,data_ft.time{1},'nearest'); + + nSamples = length(data_ft.time{1}); + otherwise + error('Unknown mode.'); +end + +data_ft.hdr.Fs = Fs; +data_ft.hdr.nSamples = nSamples; +data_ft.label = data.label; +data_ft.sampleinfo = [1 data_ft.hdr.nSamples]; \ No newline at end of file diff --git a/shared/io/datatypes/tsd/UnionTSD.m b/shared/io/datatypes/tsd/UnionTSD.m new file mode 100644 index 0000000..a87592d --- /dev/null +++ b/shared/io/datatypes/tsd/UnionTSD.m @@ -0,0 +1,88 @@ +function tsd_out = UnionTSD(cfg_in,tsd1,tsd2) +% function tsd_out = UnionTSD(cfg_in,tsd1,tsd2) +% +% returns the union (merge) of two TSD inputs +% +% example: +% +% tsd1.tvec = [1 3 5]; tsd1.data = [1 3 5]; +% +% tsd2.tvec = [2 4]; tsd1.data = [22 44]; +% +% tsd_out.tvec = [1 2 3 4 5]; tsd_out.data = [1 22 3 44 5]; +% +% - throws error if any of the tvec values are identical +% - output tvec will be sorted +% +% MvdM 2016-04-15 initial version +% NOTE merging of usr fields doesn't work well if one input is empty + +cfg_def = []; +cfg = ProcessConfig2(cfg_def,cfg_in); + +mfun = mfilename; + +%%%%%%%%%%%%%%%%%%%%%% +%%% INPUT CHECKING %%% +%%%%%%%%%%%%%%%%%%%%%% + +if ~CheckTSD(tsd1) | ~CheckTSD(tsd2) + error('Malformed TSD'); +end + +% check if any timestamps are identical; if so, give up +common_values = intersect(tsd1.tvec,tsd2.tvec); +if ~isempty(common_values) + error('tsd1 and tsd2 have common tvec value(s)'); +end + +% check if any usr field names are different; if so, give up +fn1 = []; fn2 = []; +if xor(isfield(tsd1,'usr'),isfield(tsd2,'usr')) & (~isempty(tsd1.tvec) & ~isempty(tsd2.tvec)) % OK if one is empty... + error('tsd1/tsd2 usr field presence mismatch'); +elseif isfield(tsd1,'usr') & isfield(tsd2,'usr') + fn1 = fieldnames(tsd1.usr); fn2 = fieldnames(tsd2.usr); + fnx = setxor(fn1,fn2); % names that are in one, but not the other + if ~isempty(fnx) + error('tsd1 and tsd2 have different usr field names'); + end +end + +%%%%%%%%%%%%% +%%% MERGE %%% +%%%%%%%%%%%%% + +tsd_out = tsd; % create blank tsd + +% construct new row tvec +if ~isrow(tsd1.tvec), tsd1.tvec = tsd1.tvec'; end +if ~isrow(tsd2.tvec), tsd2.tvec = tsd2.tvec'; end + +tsd_out.tvec = cat(2,tsd1.tvec,tsd2.tvec); + +% construct new data -- first check shape +if size(tsd1.data,1) > size(tsd1.data,2) + fprintf('WARNING: UnionTSD.m: tsd1 has more nSignals than nSamples!\n'); +end + +if size(tsd2.data,1) > size(tsd2.data,2) + fprintf('WARNING: UnionTSD.m: tsd2 has more nSignals than nSamples!\n') +end + +tsd_out.data = cat(2,tsd1.data,tsd2.data); + +[tsd_out.tvec,sort_idx] = sort(tsd_out.tvec,'ascend'); +tsd_out.data = tsd_out.data(:,sort_idx); + +% handle usr fields +for iF = 1:length(fn1) + + tsd_out.usr.(fn1{iF}) = cat(2,tsd1.usr.(fn1{iF}),tsd2.usr.(fn1{iF})); + tsd_out.usr.(fn1{iF}) = tsd_out.usr.(fn1{iF})(:,sort_idx); + +end + +% housekeeping +tsd_out.cfg.history.mfun = cat(1,tsd_out.cfg.history.mfun,mfun); +tsd_out.cfg.history.cfg = cat(1,tsd_out.cfg.history.cfg,{cfg}); + diff --git a/shared/io/datatypes/tsd/decimate_tsd.m b/shared/io/datatypes/tsd/decimate_tsd.m new file mode 100644 index 0000000..24759b3 --- /dev/null +++ b/shared/io/datatypes/tsd/decimate_tsd.m @@ -0,0 +1,29 @@ +function tsd_in = decimate_tsd(cfg_in,tsd_in) +% function tsd_in = decimate_tsd(cfg,tsd_in) +% +% +% DOES NOT YET WORK FOR MULTIDIMENSIONAL DATA + +cfg = []; +cfg.decimateFactor = 4; +cfg = ProcessConfig(cfg, cfg_in); + +mfun = mfilename; + +if ~CheckTSD(tsd_in) + error('Input data is not correctly formed.'); +end + +for iData = size(tsd_in.data,1):-1:1 + + data_temp(iData,:) = decimate(tsd_in.data(iData,:),cfg.decimateFactor); + tsd_in.cfg.hdr{iData}.SamplingFrequency = tsd_in.cfg.hdr{iData}.SamplingFrequency ./ cfg.decimateFactor; + +end + +tsd_in.data = data_temp; +tsd_in.tvec = downsample(tsd_in.tvec,cfg.decimateFactor); + +% housekeeping +tsd_in.cfg.history.mfun = cat(1,tsd_in.cfg.history.mfun,mfun); +tsd_in.cfg.history.cfg = cat(1,tsd_in.cfg.history.cfg,{cfg}); \ No newline at end of file diff --git a/shared/io/datatypes/tsd/tsd.m b/shared/io/datatypes/tsd/tsd.m new file mode 100644 index 0000000..6f24d30 --- /dev/null +++ b/shared/io/datatypes/tsd/tsd.m @@ -0,0 +1,57 @@ +function tsd_out = tsd(varargin) +% TSD Time-stamped data datatype constructor +% Time-stamped data, or tsd, is one of the main data types. Each tsd +% struct contains a time vector and corresponding data (such as voltages +% or position). +% +% function tsd_out = TSD(tvec,data) +% +% function tsd_out = TSD(tvec,data,'label') +% +% INPUTS: +% tvec: time vector (time stamps in ascending order) +% data: measurements taken at each time stamp +% label: unique identifier +% units: measurement units +% +% OUTPUTS +% iv_out: iv struct with fields: +% .type - 'tsd'; datatype identification +% .units - 'au'; unit marker, e.g., 'cm','px','mV' +% .tvec - [nx1] double, where n is the number of time samples taken +% .data - [1xn] double containing measures taken at each sample +% point +% .label - Unique identifier for the data contained with the tsd, +% such as a record of the data source (in the case of a +% CSC, the label corresponds to the tetrode that recorded +% the signal) +% .cfg - record of the data's history including config +% parameters and functions visited +% +% MvdM 2014-06-17 +% aacarey edit, Nov 2015 +% youkitan edit Feb 2017, unit marker + +tsd_out.type = 'tsd'; +tsd_out.units = 'au'; +tsd_out.tvec = []; +tsd_out.data = []; +tsd_out.label = {}; + +if nargin == 2 + tsd_out.tvec = varargin{1}; + tsd_out.data = varargin{2}; +elseif nargin == 3 + tsd_out.tvec = varargin{1}; + tsd_out.data = varargin{2}; + tsd_out.label = varargin{3}; +elseif nargin == 4 + tsd_out.tvec = varargin{1}; + tsd_out.data = varargin{2}; + tsd_out.label = varargin{3}; + tsd_out.units = varargin{4}; +end + +% housekeeping +tsd_out.cfg.history.mfun{1} = mfilename; +tsd_out.cfg.history.cfg{1} = []; \ No newline at end of file diff --git a/shared/io/datatypes/tsd/zscore_tsd.m b/shared/io/datatypes/tsd/zscore_tsd.m new file mode 100644 index 0000000..94ab3ca --- /dev/null +++ b/shared/io/datatypes/tsd/zscore_tsd.m @@ -0,0 +1,23 @@ +function tsd_in = zscore_tsd(tsd_in) +% function tsd_out = zscore_tsd(tsd_in) +% +% z-scores data +% +% MvdM 2015-11-03 added functionality for multiple data channels + +cfg = []; +mfun = mfilename; + +nDim = size(tsd_in.data,1); + +for iDim = 1:nDim + + keep_idx = ~isnan(tsd_in.data(iDim,:)); + + tsd_in.data(iDim,keep_idx) = zscore(tsd_in.data(iDim,keep_idx)); + +end + +% housekeeping +tsd_in.cfg.history.mfun = cat(1,tsd_in.cfg.history.mfun,mfun); +tsd_in.cfg.history.cfg = cat(1,tsd_in.cfg.history.cfg,{cfg}); \ No newline at end of file diff --git a/shared/proc/FilterLFP.m b/shared/proc/FilterLFP.m new file mode 100644 index 0000000..fdebfc3 --- /dev/null +++ b/shared/proc/FilterLFP.m @@ -0,0 +1,129 @@ +function lfp_tsd = FilterLFP(cfg_in,lfp_tsd) +% function lfp_tsd = FilterLFP(cfg,lfp_tsd) +% +% INPUT: +% tsd with LFP data +% +% OUTPUT: +% tsd with filtered LFP data +% +% CFG OPTIONS with defaults: +% +% cfg.type = 'butter'; % {'cheby1','butter','fdesign'} -- type of filter to be used +% cfg.order = 4; % filter order +% cfg.display_filter = 0; % show output of fvtool on filter +% cfg.bandtype = 'bandpass'; % 'highpass', 'lowpass' +% cfg.R = 0.5; % passband ripple (in dB) for Chebyshev filters only +% cfg.f = [6 10]; filter range to use (in Hz) +% cfg.verbose = 1; If 1 display helpful text in command window, if 0 don't +% +% MvdM 2014-06-24, 25 (update cfg_in) + +cfg_def.type = 'butter'; +cfg_def.order = 4; +cfg_def.display_filter = 0; +cfg_def.band = 'bandpass'; % not used +cfg_def.R = 0.5; % passband ripple (in dB) for Chebyshev filters only +cfg_def.f = [6 10]; +cfg_def.verbose = 1; + +mfun = mfilename; +cfg = ProcessConfig(cfg_def,cfg_in,mfun); % this takes fields from cfg_in and puts them into cfg + +% do some checks on the data +if ~CheckTSD(lfp_tsd) + error('FilterLFP.m: CheckTSD failed.'); +end + +% check reported Fs in headers +nSignals = length(lfp_tsd.cfg.hdr); +for iS = nSignals:-1:1 + + if isfield(lfp_tsd.cfg.hdr{iS},'SamplingFrequency') + reported_Fs(iS) = lfp_tsd.cfg.hdr{iS}.SamplingFrequency; + elseif isfield(lfp_tsd.cfg.hdr{iS},'Fs') + reported_Fs(iS) = lfp_tsd.cfg.hdr{iS}.Fs; + else + error('Unknown Fs.'); + end + + +end + +reported_Fs = unique(reported_Fs); +if length(reported_Fs) > 1 + error('FilterLFP.m: multiple sampling frequencies in header.'); +end + +tvec_diffs = diff(lfp_tsd.tvec); +median_Fs = 1./median(tvec_diffs); + +if cfg.verbose; fprintf('FilterLFP.m: reported Fs %.2f, median tvec Fs %.2f.\n',reported_Fs,median_Fs); end + +Fs = reported_Fs; + +% construct the filter +Wn = cfg.f ./ (Fs/2); % convert to units of half sample rate +switch cfg.type + + case 'butter' + + %[z,p,k] = butter(cfg.order,Wn); + [b,a] = butter(cfg.order,Wn); + + case 'cheby1' + + %[z,p,k] = cheby1(cfg.order,cfg.R,Wn); + [b,a] = cheby1(cfg.order,cfg.R,Wn); + + case 'fdesign' + + d = fdesign.bandpass('N,F3dB1,F3dB2',cfg.order,cfg.f(1),cfg.f(2),Fs); + Hd = design(d,'butter'); + b = Hd.sosMatrix; a = Hd.scaleValues; +end + +% convert to SOS format +%[sos,g] = zp2sos(z,p,k); + +% display if requested +if cfg.display_filter + %h = dfilt.df2sos(sos,g); + %fvtool(h); + if exist('Hd','var') + fvtool(Hd); + else + fvtool(b,a); + end + fprintf('FilterLFP.m: paused, press key to continue...\n'); + pause(1); +end + +% process signals +for iS = 1:nSignals + + if cfg.verbose; fprintf('FilterLFP.m: filtering signal %d/%d...\n',iS,nSignals); end + + % check for NaNs in the data; if there are, issue a warning and replace by + % zeros + temp_sig = lfp_tsd.data(iS,:); + + nan_idx = find(isnan(temp_sig)); + + if ~isempty(nan_idx) + fprintf('WARNING: FilterLFP.m: signal %d contains NaNs (%d).\n',iS,length(nan_idx)); + temp_sig(nan_idx) = 0; + end + + % filter + %temp_sig = filtfilt(sos,g,temp_sig); + temp_sig = filtfilt(b,a,temp_sig); + + % reinstate NaNs and put signal back into tsd + temp_sig(nan_idx) = NaN; + lfp_tsd.data(iS,:) = temp_sig; +end + +% housekeeping +lfp_tsd.cfg.history.mfun = cat(1,lfp_tsd.cfg.history.mfun,mfun); +lfp_tsd.cfg.history.cfg = cat(1,lfp_tsd.cfg.history.cfg,{cfg}); \ No newline at end of file diff --git a/workflows/sandbox_ACF_gamma_events.m b/workflows/sandbox_ACF_gamma_events.m new file mode 100755 index 0000000..ee5beb8 --- /dev/null +++ b/workflows/sandbox_ACF_gamma_events.m @@ -0,0 +1,245 @@ +%% script exploring spike autocorrelations, spike spectra, and spike-triggered averages +% remember to set path +clear +addpath(genpath('/Users/jericcarmichael/Documents/GitHub/EC_Multisite')) +addpath(genpath('/Users/jericcarmichael/Documents/GitHub/striatal-spike-rhythms'))%chronux_2_12\spectral_analysis')); +% fd = {'D:\data\adrlab\R117\R117-2007-06-20'}; +% fd = {'C:\data\adrlab\R117-2007-06-20'}; +% fd = {'C:\data\R042\R042-2013-08-17'}; +fd = {'/Volumes/Fenrir/Str_rhythms/R117-2007-06-20'}; % EC dta location +cd(fd{1}); + +LoadExpKeys; +cfg = []; cfg.fc = ExpKeys.goodGamma_vStr; +% cfg = []; cfg.fc = ExpKeys.goodTheta; +csc = LoadCSC(cfg); csc.data = csc.data-nanmean(csc.data); % may need to locdetrend + +%% load data +S = LoadSpikes([]); + +%% artificial "gamma cell" which is firing at 55 Hz +S.label{end+1} = 'artificial'; +S.t{end+1} = [1:(csc.cfg.hdr{1}.SamplingFrequency/55)*0.001:length(csc.tvec)/4]'; +CheckTS(S) +%% detect gamma events +cfg_gamma_detect = []; +cfg_gamma_detect.detect_thr = [0.7 .7 ]; % threshold for event detection: 95th percentile of (amplitude) envelope + +gamma_iv = MS_DetectEvents(cfg_gamma_detect, csc, ExpKeys); + +%% detect NON-gamma events + +cfg_other = []; +cfg_other.epoch = 'all'; % chewing occurs during task mostly +cfg_other.minlen = 0.005; +cfg_other.filter_cfg.f = [45 65]; % default [25 35] +cfg_other.threshold = .5; % 0.25 for session 2, 0.5 for session 1? +cfg_other.smooth = 0.05; % convolve with Gaussian of this SD +cfg_other.dcn = '<'; +non_gamma_iv = DetectEvents(cfg_other, csc, ExpKeys); + +non_gamma_iv.usr.evt_len = (non_gamma_iv.tend - non_gamma_iv.tstart)'; + +% remove long events +cfg_max_len = []; +cfg_max_len.operation = '<'; +cfg_max_len.threshold = max(gamma_iv.low.tend - gamma_iv.low.tstart); +non_gamma_iv = SelectIV(cfg_max_len,non_gamma_iv,'evt_len'); + +% remove events that are too short +cfg_max_len = []; +cfg_max_len.operation = '>'; +cfg_max_len.threshold = min(gamma_iv.low.tend - gamma_iv.low.tstart); +evt = SelectIV(cfg_max_len,non_gamma_iv,'evt_len'); + +% select a random assortment of the non-gamma events to get the same number +% of events +keep_idx = datasample(1:length(non_gamma_iv.tend),length(gamma_iv.low.tend)); % get nSamples from gamma events chosen randomly from the number of samples in non gamma +non_gamma_iv = SelectIV([],non_gamma_iv,keep_idx); + +% PlotTSDfromIV([], non_gamma_iv, csc) + +%% try multiraster +% cfg = []; +% cfg.evt = gamma_iv.low; +% cfg.lfp = csc; +% cfg.spkColor = 'jet'; +% MultiRaster(cfg,S) + +%% get ACF for each gamma event + +this_band = 'low'; +% +cfg_acf = []; +cfg_acf.binsize = .001; +cfg_acf.maxlag = 200; +all_acf = cell(1,length(S.t)); +for iEvt = 1:500%length(gamma_iv.(this_band).tstart) % loop events + disp(num2str(iEvt)) + this_event = restrict(S, gamma_iv.(this_band).tstart(iEvt), gamma_iv.(this_band).tend(iEvt)); + % create bins for this gamma event. + tbin_edges = gamma_iv.(this_band).tstart(iEvt):cfg_acf.binsize:gamma_iv.(this_band).tend(iEvt); % vector of time bin edges (for histogram) + tbin_centers = tbin_edges(1:end-1)+cfg_acf.binsize/2; % vector of time bin centers (for plotting) + + % compute + for iC = 1:length(S.t) % loop cells + if ~isempty(this_event.t{iC}) + if length(this_event.t{iC})>1 + % bin the data for this cell. + spk_count = histc(S.t{iC},tbin_edges); % get spike counts for each bin + spk_count = spk_count(1:end-1); + + [all_acf{iC}(iEvt,:) ,tvec_acf] = ComputeACF(cfg_acf, spk_count); + + end + end + end +end + +%% same thing for NON-GAMMA epochs +cfg_acf = []; +cfg_acf.binsize = .001; +cfg_acf.maxlag = 200; +all_acf_non = cell(1,length(S.t)); +for iEvt = 1:500%length(non_gamma_iv.tstart) % loop events + this_event = restrict(S, non_gamma_iv.tstart(iEvt), non_gamma_iv.tend(iEvt)); + % create bins for this gamma event. + tbin_edges = non_gamma_iv.tstart(iEvt):cfg_acf.binsize:non_gamma_iv.tend(iEvt); % vector of time bin edges (for histogram) + tbin_centers = tbin_edges(1:end-1)+cfg_acf.binsize/2; % vector of time bin centers (for plotting) + + % compute + for iC = 1:length(S.t) % loop cells + if ~isempty(this_event.t{iC}) + if length(this_event.t{iC})>1 + % bin the data for this cell. + spk_count = histc(S.t{iC},tbin_edges); % get spike counts for each bin + spk_count = spk_count(1:end-1); + + [all_acf_non{iC}(iEvt,:) ,tvec] = ComputeACF(cfg_acf, spk_count); + end + end + end +end + +%% remove cells that have very few data points +min_points = 10; +cells_to_keep =[]; +for iC = 1:length(all_acf) + if max(sum(all_acf{iC})) > min_points && max(sum(all_acf_non{iC})) > min_points + cells_to_keep =[cells_to_keep iC]; + else + fprintf('\nCell %.0f had < %.0f number of data points\n', iC, min_points) + end +end + +%% plot all the gamma events vs non events. + +figure(222) +x_plot = 1; +for iC = 1:length(cells_to_keep) + subplot(4,4,x_plot) + +plot(tvec, nanmean(all_acf{cells_to_keep(iC)}),'r', tvec, nanmean(all_acf_non{cells_to_keep(iC)}), 'k') +title(num2str(cells_to_keep(iC))) +x_plot = x_plot+1; +end + +%% plots? +figure(100) +plot_n = 0; +Fs = 1./cfg_acf.binsize; +S_gamma = restrict(S, gamma_iv.(this_band).tstart, gamma_iv.(this_band).tend); +S_non = restrict(S, non_gamma_iv.tstart, non_gamma_iv.tend); + +%% +figure(100) +plot_n = 0; +for iC =cells_to_keep(1:ceil(length(cells_to_keep)/2)) + + plot_n = plot_n+1; + subplot(4,length(cells_to_keep), plot_n) + hold on + plot(tvec_acf, sum(all_acf_non{iC}), 'r') + plot(tvec_acf, sum(all_acf{iC}), 'b') + xlim([-0.05 0.05]) + ylim([0 inf]) + + % plot the PSD of the ACORR + subplot(4, length(cells_to_keep), plot_n+length(cells_to_keep)) + + [P,F] = pwelch(all_acf{iC}(1:floor(length(all_acf{iC})/2)),2^7,2^6,[],Fs); + + plot(F,10*log10(P)); xlim([0 100]); + title('acorr spectrum'); + + % get the sta in gamma events + this_spk_binned = zeros(size(csc.tvec)); + idx = nearest_idx3(S_gamma.t{iC},csc.tvec); + this_spk_binned(idx) = 1; + [xc_gamma,tvec] = xcorr(csc.data,this_spk_binned,1000); + tvec = tvec.*median(diff(csc.tvec)); + % same for non gamma + this_spk_binned = zeros(size(csc.tvec)); + idx = nearest_idx3(S_non.t{iC},csc.tvec); + this_spk_binned(idx) = 1; + [xc_non,tvec] = xcorr(csc.data,this_spk_binned,1000); + + subplot(4, length(cells_to_keep), plot_n+length(cells_to_keep)*2) + + plot(tvec,xc_non,'r', tvec, xc_gamma, 'b'); + + title('raw STA'); + +end +tightfig +maximize +%% +figure(200) +plot_n = 0; +for iC =cells_to_keep(ceil(length(cells_to_keep)/2):end) + + plot_n = plot_n+1; + subplot(4,ceil(length(cells_to_keep)/2), plot_n) + hold on + plot(tvec_acf, sum(all_acf_non{iC}), 'r') + plot(tvec_acf, sum(all_acf{iC}), 'b') + xlim([-0.05 0.05]) + ylim([0 inf]) + + % plot the PSD of the XCORR + subplot(4, ceil(length(cells_to_keep)/2), plot_n+ceil(length(cells_to_keep)/2)) + + [P,F] = pwelch(all_acf{iC}(1:floor(length(all_acf{iC})/2)),2^10,2^9,[],Fs); + + plot(F,10*log10(P)); xlim([0 100]); + title('acorr spectrum'); + + % get the sta in gamma events + this_spk_binned = zeros(size(csc.tvec)); + idx = nearest_idx3(S_gamma.t{iC},csc.tvec); + this_spk_binned(idx) = 1; + [xc_gamma,tvec] = xcorr(csc.data,this_spk_binned,1000); + tvec = tvec.*median(diff(csc.tvec)); + % same for non gamma + this_spk_binned = zeros(size(csc.tvec)); + idx = nearest_idx3(S_non.t{iC},csc.tvec); + this_spk_binned(idx) = 1; + [xc_non,tvec] = xcorr(csc.data,this_spk_binned,1000); + + subplot(4, ceil(length(cells_to_keep)/2), plot_n+(ceil(length(cells_to_keep)/2))*2) + + plot(tvec,xc_non,'r', tvec, xc_gamma, 'b'); + + title('raw STA'); + +end +tightfig +maximize + +%% get the psd for each STA across all events. + + + + + + diff --git a/workflows/sandbox_spikeSpectra.m b/workflows/sandbox_spikeSpectra.m index 0a9969e..a240919 100644 --- a/workflows/sandbox_spikeSpectra.m +++ b/workflows/sandbox_spikeSpectra.m @@ -3,14 +3,15 @@ clear %addpath(genpath('D:\My_Documents\GitHub\striatal-spike-rhythms\chronux_2_12\spectral_analysis')); -%fd = {'D:\data\adrlab\R117\R117-2007-06-20'}; -fd = {'C:\data\adrlab\R117-2007-06-20'}; -fd = {'C:\data\R042\R042-2013-08-17'}; +% fd = {'D:\data\adrlab\R117\R117-2007-06-20'}; +% fd = {'C:\data\adrlab\R117-2007-06-20'}; +% fd = {'C:\data\R042\R042-2013-08-17'}; +fd = {'/Volumes/Fenrir/Str_rhythms/R117-2007-06-20'}; % EC dta location cd(fd{1}); LoadExpKeys; -%cfg = []; cfg.fc = ExpKeys.goodGamma_vStr; -cfg = []; cfg.fc = ExpKeys.goodTheta; +cfg = []; cfg.fc = ExpKeys.goodGamma_vStr; +% cfg = []; cfg.fc = ExpKeys.goodTheta; csc = LoadCSC(cfg); csc.data = csc.data-nanmean(csc.data); % may need to locdetrend %% load data @@ -161,12 +162,12 @@ keep = pval < 0.05; hold on; plot(xval(keep),cc(keep),'.b','MarkerSIze',20); %% -function [acf,tvec] = ComputeACF(cfg,spk_binned) -if isfield(cfg,'maxlag') - [acf,tvec] = xcorr(spk_binned,spk_binned,cfg.maxlag); -else - [acf,tvec] = xcorr(spk_binned,spk_binned); -end -tvec = tvec.*cfg.binsize; -acf(ceil(length(acf)/2)) = 0; -end \ No newline at end of file +% function [acf,tvec] = ComputeACF(cfg,spk_binned) +% if isfield(cfg,'maxlag') +% [acf,tvec] = xcorr(spk_binned,spk_binned,cfg.maxlag); +% else +% [acf,tvec] = xcorr(spk_binned,spk_binned); +% end +% tvec = tvec.*cfg.binsize; +% acf(ceil(length(acf)/2)) = 0; +% end \ No newline at end of file