preprocessData.m

function [rez, DATA, uproj] = preprocessData(ops)
tic;
uproj = [];
ops.nt0 	= getOr(ops, {'nt0'}, 61);


% if strcmp(ops.datatype , 'openEphys')
%    ops = convertOpenEphysToRawBInary(ops);  % convert data, only for OpenEphys
% end
% if strcmp(ops.datatype, 'h5')
%   ops = convertWillowToRawBInary(ops);      % convert willow data
% end

if ~isempty(ops.chanMap)
    if ischar(ops.chanMap)
        load(ops.chanMap);
        try
            chanMapConn = chanMap(connected>1e-6);
            xc = xcoords(connected>1e-6);
            yc = ycoords(connected>1e-6);
        catch
            chanMapConn = 1+chanNums(connected>1e-6);
            xc = zeros(numel(chanMapConn), 1);
            yc = [1:1:numel(chanMapConn)]';
            disp('WTF1!');
            keyboard
        end
        ops.Nchan    = getOr(ops, 'Nchan', sum(connected>1e-6));
        ops.NchanTOT = getOr(ops, 'NchanTOT', 1024);
        %% JPK: fs is loaded from MAT file
        if exist('fs', 'var')
            ops.fs       = getOr(ops, 'fs', fs);
        end
    else
        chanMap = ops.chanMap;
        chanMapConn = ops.chanMap;
        xc = zeros(numel(chanMapConn), 1);
        yc = [1:1:numel(chanMapConn)]';
        connected = true(numel(chanMap), 1);      
        
        ops.Nchan    = numel(connected);
        sprintf('rez.ops.Nchan: %i', ops.Nchan);
        ops.NchanTOT  = 1024;
        disp('WTF2!');
        keyboard
    end
else
    chanMap  = 1:ops.Nchan;
    connected = true(numel(chanMap), 1);
    
    chanMapConn = 1:ops.Nchan;    
    xc = zeros(numel(chanMapConn), 1);
    yc = [1:1:numel(chanMapConn)]';
    disp('WTF3!');
    keyboard
end
if exist('kcoords', 'var')
    % JPK: This should do nothing.
    kcoords = kcoords(find(connected==1));
else
    kcoords = ones(ops.Nchan, 1);
end
NchanTOT = ops.NchanTOT;
NT       = ops.NT ;

rez.ops         = ops;
rez.xc = xc;
rez.yc = yc;
rez.xcoords = xcoords;
rez.ycoords = ycoords;
rez.connected   = connected;
rez.ops.chanMap = chanMap;
rez.ops.kcoords = kcoords; 

d = dir(ops.fbinary);
ops.sampsToRead = floor(d.bytes/NchanTOT/2);

if ispc
    dmem         = memory;
    memfree      = dmem.MemAvailableAllArrays/8;
    memallocated = min(ops.ForceMaxRAMforDat, dmem.MemAvailableAllArrays) - memfree;
    memallocated = max(0, memallocated);
else
    memallocated = ops.ForceMaxRAMforDat;
end
nint16s      = memallocated/2;

NTbuff      = NT + 4*ops.ntbuff;
Nbatch      = ceil(d.bytes/2/NchanTOT /(NT-ops.ntbuff));
Nbatch_buff = floor(4/5 * nint16s/rez.ops.Nchan /(NT-ops.ntbuff)); % factor of 4/5 for storing PCs of spikes
Nbatch_buff = min(Nbatch_buff, Nbatch);

%% subtract virtual reference
% 
% % open input data file
% fid = fopen(ops.fbinary, 'r');
% buff = fread(fid, '*int16');
% fclose(fid);
% M = 1024;
% N = length(buff)/M;
% data=reshape(buff,M,N);
% % only keep good channels
% tmp = data(chanMapConn+1,:);
% data = single(data)-mean(data(chanMapConn+1,:)); 

%% load data into patches, filter, compute covariance

% compute filter weights
if isfield(ops,'fslow')&&ops.fslow<ops.fs/2
    [b1, a1] = butter(3, [ops.fshigh/ops.fs,ops.fslow/ops.fs]*2, 'bandpass');
else
    [b1, a1] = butter(3, ops.fshigh/ops.fs*2, 'high');
end

fprintf('Time %3.0fs. Loading raw data... \n', toc);

% open input data file
fid = fopen(ops.fbinary, 'r');

ibatch = 0;
Nchan = rez.ops.Nchan;
NchanTOT = rez.ops.NchanTOT;
if ops.GPU
    CC = gpuArray.zeros( Nchan,  Nchan, 'single');
else
    CC = zeros( Nchan,  Nchan, 'single');
end
%if strcmp(ops.whitening, 'noSpikes')
%    if ops.GPU
%        nPairs = gpuArray.zeros( Nchan,  Nchan, 'single');
%    else
%        nPairs = zeros( Nchan,  Nchan, 'single');
%    end
%end
if ~exist('DATA', 'var')
    DATA = zeros(NT, rez.ops.Nchan, Nbatch_buff, 'int16');
end

isproc = zeros(Nbatch, 1);
while 1
    ibatch = ibatch + ops.nSkipCov;
    
    % jump to correct place in input data file
    offset = max(0, 2*NchanTOT*((NT - ops.ntbuff) * (ibatch-1) - 2*ops.ntbuff));
    if ibatch==1
        ioffset = 0;
    else
        ioffset = ops.ntbuff;
    end
    fseek(fid, offset, 'bof');

    % read from input data file in column order
    % data read is interpreted as int16 and saved as int16
    % HOW MUCH IS READ?
    % NchanTOT = 1024
    % NTbuff = (32*1024 + 64) + 4*64 = 33088 (if samples, then 1 second of data)
    buff = fread(fid, [NchanTOT NTbuff], '*int16');
    %% JPK: plot some data to make sure that it looks good

    %keyboard
    
    % if no more data, then exit while loop
    if isempty(buff)
        break;
    end

    % if remaining data is less than buffer, then fill buffer with copies of last sample
    nsampcurr = size(buff,2);
    if nsampcurr<NTbuff
        buff(:, nsampcurr+1:NTbuff) = repmat(buff(:,nsampcurr), 1, NTbuff-nsampcurr);
    end

    % save data to gpu or cpu (variable)
    if ops.GPU
        dataRAW = gpuArray(buff);
    else
        dataRAW = buff;
    end

    % transpose data
    % assumes...
    % was rows = channels, columns = time
    % now rows = time, columns = channels
    dataRAW = dataRAW';
    % convert to single precision
    dataRAW = single(dataRAW);
    % only keep subset of channels (columns)
    dataRAW = dataRAW(:, chanMapConn+1);
    
    % subtractvirtual reference
    %dataRAW1 = dataRAW-mean(dataRAW,2);
    
    virtref = mean(dataRAW,2);
    virtpower = sum(virtref.^2);
    a = size(dataRAW);
    for i=1:a(2)
        overlap = sum(virtref.*dataRAW(:,i))/virtpower;
        dataRAW(:,i) = dataRAW(:,i)-overlap*virtref;
    end
    
    % JPK: normalize data so standard deviation==1
    % since signals are strongly non-gaussian,
    % use formula for std dev based on median
    %ops.scaleproc = mean(median(abs(dataRAW)))/0.6745;
    %keyboard
    
    datr = filter(b1, a1, dataRAW);
    datr = flipud(datr);
    datr = filter(b1, a1, datr);
    datr = flipud(datr);
    %keyboard
    switch ops.whitening
        case 'noSpikes'
            smin      = my_min(datr, ops.loc_range, [1 2]);
            sd = std(datr, [], 1);
            peaks     = single(datr<smin+1e-3 & bsxfun(@lt, datr, ops.spkTh * sd));
            blankout  = 1+my_min(-peaks, ops.long_range, [1 2]);
            smin      = datr .* blankout;
            CC        = CC + (smin' * smin)/NT;
            nPairs    = nPairs + (blankout'*blankout)/NT;
        otherwise
            CC        = CC + (datr' * datr)/NT;
            %disp('This ran1.');
    end
    
    if ibatch<=Nbatch_buff
        DATA(:,:,ibatch) = gather_try(int16( datr(ioffset + (1:NT),:)));
        isproc(ibatch) = 1;
    end
end
CC = CC / ceil((Nbatch-1)/ops.nSkipCov);
%switch ops.whitening
%    case 'noSpikes'
%        nPairs = nPairs/ibatch;
%end
fclose(fid);
fprintf('Time %3.0fs. Channel-whitening filters computed. \n', toc);
switch ops.whitening
    case 'diag'
        CC = diag(diag(CC));
        %disp('This ran2.');
    case 'noSpikes'
        CC = CC ./nPairs;
end

if ops.whiteningRange<Inf
    ops.whiteningRange = min(ops.whiteningRange, Nchan);
    Wrot = whiteningLocal(gather(CC), yc, xc, ops.whiteningRange);
else
    %
    [E, D] 	= svd(CC);
    D = diag(D);
    eps 	= 1e-6;
    Wrot 	= E * diag(1./(D + eps).^.5) * E';
    disp('WTF!')
    keyboard
end
Wrot    = ops.scaleproc * Wrot;

%keyboard


fprintf('Time %3.0fs. Loading raw data and applying filters... \n', toc);

fid  = fopen(ops.fbinary, 'r');
fidW = fopen(ops.fproc, 'w');

if strcmp(ops.initialize, 'fromData')
    i0  = 0;
    ixt  = round(linspace(1, size(ops.wPCA,1), ops.nt0));
    wPCA = ops.wPCA(ixt, 1:3);
    
    rez.ops.wPCA = wPCA; % write wPCA back into the rez structure
    uproj = zeros(1e6,  size(wPCA,2) * Nchan, 'single');
end
%
for ibatch = 1:Nbatch
    if isproc(ibatch) %ibatch<=Nbatch_buff
        %disp('This ran2.');
        if ops.GPU
            datr = single(gpuArray(DATA(:,:,ibatch)));
        else
            datr = single(DATA(:,:,ibatch));
        end
    else
        offset = max(0, 2*NchanTOT*((NT - ops.ntbuff) * (ibatch-1) - 2*ops.ntbuff));
        if ibatch==1
            ioffset = 0;
        else
            ioffset = ops.ntbuff;
        end
        fseek(fid, offset, 'bof');
        
        buff = fread(fid, [NchanTOT NTbuff], '*int16');
        if isempty(buff)
            break;
        end
        nsampcurr = size(buff,2);
        if nsampcurr<NTbuff
            buff(:, nsampcurr+1:NTbuff) = repmat(buff(:,nsampcurr), 1, NTbuff-nsampcurr);
        end
        
        if ops.GPU
            dataRAW = gpuArray(buff);
        else
            dataRAW = buff;
        end
        dataRAW = dataRAW';
        dataRAW = single(dataRAW);
        % TODO why do we not need chanMapConn+1 here like we did before?
        dataRAW = dataRAW(:, chanMapConn+1);
        
        datr = filter(b1, a1, dataRAW);
        datr = flipud(datr);
        datr = filter(b1, a1, datr);
        datr = flipud(datr);
        
        datr = datr(ioffset + (1:NT),:);
    end
    
    %foo    = datr * Wrot;
    datr    = datr * Wrot;
    
    %keyboard
    
    if ops.GPU
        dataRAW = gpuArray(datr);
    else
        dataRAW = datr;
    end
    %         dataRAW = datr;
    dataRAW = single(dataRAW);
    % TODO eliminate ops.scaleproc
    dataRAW = dataRAW / ops.scaleproc;

    if strcmp(ops.initialize, 'fromData') %&& rem(ibatch, 10)==1
        % find isolated spikes
        [row, col, mu] = isolated_peaks(dataRAW, ops.loc_range, ops.long_range, ops.spkTh);
        
        % find their PC projections
        uS = get_PCproj(dataRAW, row, col, wPCA, ops.maskMaxChannels);
        
        uS = permute(uS, [2 1 3]);
        uS = reshape(uS,numel(row), Nchan * size(wPCA,2));
        
        if i0+numel(row)>size(uproj,1)
            uproj(1e6 + size(uproj,1), 1) = 0;
        end
        
        uproj(i0 + (1:numel(row)), :) = gather_try(uS);
        i0 = i0 + numel(row);
    end
    
    if ibatch<=Nbatch_buff
        DATA(:,:,ibatch) = gather_try(datr);
    else
        datcpu  = gather_try(int16(datr));
        fwrite(fidW, datcpu, 'int16');
        fprintf('Wrote whitened data to temporary file...');
    end
    
end

datcpu  = gather_try(int16(datr));
fwrite(fidW, datcpu, 'int16');
fprintf('Wrote whitened data to temporary file...');

if strcmp(ops.initialize, 'fromData')
   uproj(i0+1:end, :) = []; 
end
Wrot        = gather_try(Wrot);
rez.Wrot    = Wrot;

fclose(fidW);
fclose(fid);
if ops.verbose
    fprintf('Time %3.2f. Whitened data written to disk... \n', toc);
    fprintf('Time %3.2f. Preprocessing complete!\n', toc);
end


rez.temp.Nbatch = Nbatch;
rez.temp.Nbatch_buff = Nbatch_buff;