load processedRecordings.mat
% processedRecordings = [{'DT1_rLS_20150729_2448um'}]
    
addpath(genpath('/zpool/Dropbox/Tools'))
addpath(genpath('/zpool/Dropbox/code'))
disp('starting assemblies analysis')
for f = length(processedRecordings):-1:1

    cd(processedRecordings{f})
    xml = dir('*xml');
    SetCurrentSession(xml.name)
    units = GetUnits();

  for n = 1:size(units,1);
        spikes{n} = GetSpikeTimes(units(n,:));
  end
  nTrials = 0;
if ~exist(['assemblies.mat'])
    lfp = GetLFP(23);
    if exist('behav.mat')
    load('behav.mat')
        for i=1:length(trials)
        for j=1:length(trials{i})
            nTrials = nTrials+1;
            all_trials{nTrials}=trials{i}{j};
        end
    end
    else
    %% make fake trials here
%     trials
tt = linspace(1,length(lfp),100);
 for t =1:100
    all_trials{t} = [0, 0, 0, tt(t),  tt(t)+6*1250];
     nTrials = nTrials+1;
 end
    end
 

    pair=[];
    for i=1:length(units)
        for j=1:length(units)
            if units(i,1) ~= units(j,1)
                if i ~= j 
                pair = [pair;i j];
                end
            end
        end
    end

    sp = single(zeros(length(units),nTrials,6000));
    for i=1:length(units)
        for j=1:nTrials
        f = find(spikes{i}>all_trials{j}(1,5));
        ff = find(spikes{i}<all_trials{j}(1,5)+6);
        sp(i,j,ceil((spikes{i}(intersect(f,ff))-all_trials{j}(1,5))*1000))=1;
        end
    end
    u = unique(pair(:,1));
    c=1;
    for i=1:length(u)
        f = find(pair(:,1)==u(i));
        tic
        disp(['running assembly analysis on ' xml.name ' on ' ...
            num2str(length(f)) ' cell pairs, of ' num2str(length(pair)) ...
            ' total pairs'])
        l = size(sp,2);
        [d dc] = cell_assembly(sp,0:150,pair(f,:));
        for j =1:length(f)
            dev(:,c,:) = squeeze(mean([mean(d(:,j,:),3)],2));
            devControl(:,c,:) = squeeze(mean([mean(dc(:,j,:),3)],2));
            c=1+c;
        end
        save('assemblies.mat','dev*','pair','units')
        toc
    end
    clear sp
    save('assemblies.mat')
    clear dev* pair f
    
%     end
end

if exist('assemblies.mat')
    load('assemblies.mat','dev*','pair','units')
    lfp = GetLFP(23);
    if exist('behav.mat')
    load('behav.mat')
    nTrials = 0;
    
    for i=1:length(trials)
        for j=1:length(trials{i})
            nTrials = nTrials+1;
            all_trials{nTrials}=trials{i}{j};
        end
    end
    pair=[];
    for i=1:length(units)
        for j=1:length(units)
            if units(i,1) ~= units(j,1)
                if i ~= j 
                pair = [pair;i j];
                end
            end
        end
    end

    sp = single(zeros(length(units),nTrials,6000));
    for i=1:length(units)
        for j=1:nTrials
        f = find(spikes{i}>all_trials{j}(1,5));
        ff = find(spikes{i}<all_trials{j}(1,5)+6);
        sp(i,j,ceil((spikes{i}(intersect(f,ff))-all_trials{j}(1,5))*1000))=1;
        end
    end
    u = unique(pair(:,1));
    finished = size(dev,2);
    p=pair(finished+1:end,:);
    uu = unique(p(:,1));
    c=size(dev,2)+1;
    
    for i=1:length(uu)
        f = find(p(:,1)==uu(i));
        tic
        disp(['running assembly analysis on ' xml.name ' on ' ...
            num2str(length(f)) ' cell pairs, of ' num2str(length(p)) ...
            ' total pairs'])
        l = size(sp,2);
        [d dc] = cell_assembly(sp,0:150,p(f,:));
        for j =1:length(f)
            dev(:,c,:) = squeeze(mean([mean(d(:,j,:),3)],2));
            devControl(:,c,:) = squeeze(mean([mean(dc(:,j,:),3)],2));
            c=1+c;
        end
        save('assemblies.mat','dev*','pair','units')
        toc
    end
end


end
cd ..
end
