diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 0000000..348ceb2
Binary files /dev/null and b/.DS_Store differ
diff --git a/toefile_utils/generate_cell_groups.m b/toefile_utils/generate_cell_groups.m
new file mode 100644
index 0000000..b774966
--- /dev/null
+++ b/toefile_utils/generate_cell_groups.m
@@ -0,0 +1,49 @@
+% From population vectors, generate cell groups
+% Must activate Plex library before running
+clear
+close all
+
+load('st1215_cat_P01_S01_2ndPen_moremerged_20151102T090824_popvec.mat');
+
+%population_vectors(ncell, nwin, nstim, nbinstart);
+[ncell, nwin, nstim, n_bin_start] = size(population_vectors);
+
+%first compute average firing rate for cell:
+avg_fr = sum(sum(sum(population_vectors, 4), 3), 2)/(nwin*nstim*n_bin_start);
+
+%threshold population vectors with FR > 6*avgFR
+cell_thresh = repmat(6*avg_fr, [1, nwin, nstim, n_bin_start]);
+cell_groups = population_vectors > cell_thresh;
+
+% pull out activity during stimulus
+% windows_in_stim = (windows_starts >= 0 - bin_width & (windows_starts + bin_width) <= 1500 + bin_width)';
+% window_win = squeeze(windows_in_stim(:, 2));
+% cell_groups = cell_groups(:, window_win, :, :);
+
+% Now, build simplicial complex
+stream = api.Plex4.createExplicitSimplexStream();
+
+% Add vertices for each cell
+for i = 1:ncell
+    stream.addVertex(i);
+end
+
+%% Now go through each population vector and add simplicial complex
+for win = 1:nwin
+    for stim = 1:nstim
+        for binst = 1:n_bin_start
+            cell_group_this = squeeze(cell_groups(:, win, stim, binst));
+            cells_in_this_group = find(cell_group_this);
+            if(any(cells_in_this_group))
+                cells_in_this_group;
+                stream.addElement(cells_in_this_group);
+                stream.ensureAllFaces();
+            end
+        end
+    end
+end
+stream.finalizeStream();
+
+%% Compute homology
+persis = api.Plex4.getModularSimplicialAlgorithm(10, 7);
+persis.computeIntervals(stream)
\ No newline at end of file
diff --git a/toefile_utils/generate_population_vectors.m b/toefile_utils/generate_population_vectors.m
new file mode 100644
index 0000000..2564eb8
--- /dev/null
+++ b/toefile_utils/generate_population_vectors.m
@@ -0,0 +1,75 @@
+%cell groups
+% Generate population vectors as in Curto&Itskov 2008
+
+clear
+close all
+
+datafilename = 'st1215_cat_P01_S01_2ndPen_moremerged_20151102T090824';
+load(strcat(datafilename, '.mat'));
+
+%toedata = data_to_save.toedata;
+% parameters for analysis
+bin_width = 10; % time bin width for cell groups in ms
+n_bin_start = 5; % number of different bin start times.
+pre_stim_ms = 2000; % number of milliseconds pre-stim
+post_stim_ms = 2000; % number of milliseconds post-stim
+stim_ms = 1000; % stimuli duration in milliseconds
+
+fs = 31250;
+
+trial_duration = pre_stim_ms + post_stim_ms + stim_ms;
+nwin = floor(trial_duration / bin_width);
+
+starts_dt = floor(bin_width/n_bin_start); % equally spaced in each win.
+
+% get n_bin_start starting points for the windows within one bin time frame
+bin_starts = starts_dt.*[0:n_bin_start-1] - pre_stim_ms;
+
+windows_starts = zeros(n_bin_start, nwin);
+windows_ends = zeros(n_bin_start, nwin);
+
+windows_starts(:, 1) = bin_starts;
+windows_ends(:, 1) = bin_starts + bin_width;
+for i = 2:nwin
+    windows_starts(:, i) = windows_starts(:, i-1) + bin_width;
+    windows_ends(:, i) = windows_ends(:, i-1) + bin_width;
+end
+
+% build population vectors
+ncells = length(toedata);
+nstim = 41;
+population_vectors = zeros(ncells, nwin, nstim, n_bin_start);
+
+% For each cell, for each stim, for each trial, count the number of spikes
+% in each time bin to arrive at population vectors.
+for cell = 1:ncells
+    for stim = 1:nstim
+        spiketimes = toedata{cell,1}.stims{stim, 1}.toes;
+        ntrials = length(spiketimes);
+        for bin_start = 1:n_bin_start
+            spikecounts_thiscellstimbinst = zeros(1, 500);
+            for trial = 1:ntrials
+                spikes_this_trial = 1000*spiketimes{trial, 1}; %convert to ms
+                Nspikes = length(spikes_this_trial);
+                repspike = repmat(spikes_this_trial, [1, nwin]);
+                repbinst = repmat(windows_starts(bin_start, :), [Nspikes, 1]);
+                repbinend = repmat(windows_ends(bin_start, :), [Nspikes, 1]);
+                bool_spike = repspike >= repbinst & repspike < repbinend;
+                spikecounts_thistrial = sum(bool_spike, 1);
+                spikecounts_thiscellstimbinst = spikecounts_thiscellstimbinst + spikecounts_thistrial;
+            end
+            rate_this = spikecounts_thiscellstimbinst / bin_width;
+            population_vectors(cell, :, stim, bin_start) = rate_this;
+        end
+    end
+end
+
+% Save population vector data
+% popvec_data_to_save = rmfield(data_to_save, 'toedata');
+% popvec_data_to_save.popvecs = population_vectors;
+% popvec_data_to_save.binwidth = bin_width;
+% popvec_data_to_save.winstarts = windows_starts;
+popvecfilename = strcat(datafilename, '_popvec');
+save(popvecfilename, 'population_vectors', 'bin_width', 'windows_starts');
+
+
diff --git a/toefile_utils/kwik2mat.m b/toefile_utils/kwik2mat.m
new file mode 100755
index 0000000..403b435
--- /dev/null
+++ b/toefile_utils/kwik2mat.m
@@ -0,0 +1,177 @@
+% kwik to mat.  Get good units out of sorted kwik file with stimulus
+% information and store in one well formatted mat file
+% Brad Theilman August/September 2015
+
+clear
+close all
+
+kwikfile = './st1215/kwdfiles/pen2_latestpipeline/st1215_cat_P01_S01_2ndPen.kwik';
+
+% get information
+kwikinfo = h5info(kwikfile);
+
+% Get cluster and time information for each spike
+spike_clusters = h5read(kwikfile, '/channel_groups/0/spikes/clusters/main');
+spike_times = h5read(kwikfile, '/channel_groups/0/spikes/time_samples');
+
+% extract all the unique cluster ids:
+clusters = unique(spike_clusters);
+
+%% Build a boolean vector with 1 by a cluster id if it is a good cluster
+good_cluster_inds = zeros(1, length(clusters));
+
+for i = 1:length(clusters)
+    
+    cluster_num = clusters(i);
+    cluster_attr_path = strcat('/channel_groups/0/clusters/main/', num2str(cluster_num));
+    cluster_group = h5readatt(kwikfile, cluster_attr_path, 'cluster_group');
+    disp(cluster_group);
+    %cluster_group 2 is Good
+    % 1 is MUA
+    if cluster_group == 1
+        good_cluster_inds(i) = 1;
+        
+    end
+    
+end
+% vector containing IDs of good clusters
+good_cluster_ids = clusters(logical(good_cluster_inds));
+
+%% Now, pull all spikes that belong to good clusters
+good_spikes = cell(length(good_cluster_ids), 2);
+for i = 1:length(good_cluster_ids)
+    cluster_num = good_cluster_ids(i);
+    %pull spikes from teh good cluster
+    this_cluster_spiketimes = spike_times(spike_clusters == cluster_num);
+    good_spikes{i, 1} = cluster_num;
+    good_spikes{i, 2} = this_cluster_spiketimes;
+end
+
+%% Now get stimulus information
+
+digmark_timesamples = h5read(kwikfile, '/event_types/DigMark/time_samples');
+digmark_codes = cell2mat(h5read(kwikfile, '/event_types/DigMark/codes'));
+stim_timesamples = h5read(kwikfile, '/event_types/Stimulus/time_samples');
+stim_codes = h5read(kwikfile, '/event_types/Stimulus/codes');
+stim_text = h5read(kwikfile, '/event_types/Stimulus/text');
+
+stim_start_times = double(digmark_timesamples(digmark_codes == '<'));
+stim_end_times = double(digmark_timesamples(digmark_codes == '>'));
+intertrial_start_times = double(digmark_timesamples(digmark_codes == '('));
+intertrial_end_times = double(digmark_timesamples(digmark_codes == ')'));
+
+stim_start_filename = cell2mat(stim_text);
+stim_start_filename = stim_start_filename(2:2:end, :);
+stim_start_filename = cellstr(stim_start_filename);
+
+% stim_start_times = h5read(kwikfile, '/event_types/60/time_samples');
+% stim_start_filename = h5read(kwikfile, '/event_types/60/stim_filename');
+% stim_end_times = h5read(kwikfile, '/event_types/62/time_samples');
+% intertrial_start_times = h5read(kwikfile, '/event_types/40/time_samples');
+% intertrial_end_times = h5read(kwikfile, '/event_types/41/time_samples');
+
+stim_files_unique = unique(stim_start_filename);
+nstims = length(stim_files_unique);
+
+
+%%
+% For each stimulus, get number of trials for that stimulus
+% start building final_mat
+final_data = struct;
+final_data.nstims = nstims;
+numtrials = zeros(1, nstims);
+stim_data = cell(nstims, 1);
+
+for i = 1:nstims
+    
+    stim_entry = struct;
+    numtrials(i) = sum(strcmp(stim_files_unique(i), stim_start_filename));
+    
+    stim_start_times_this_stim = stim_start_times(strcmp(stim_files_unique(i), stim_start_filename));
+    stim_end_times_this_stim = stim_end_times(strcmp(stim_files_unique(i), stim_start_filename));
+    stim_entry.name = stim_files_unique(i);
+    stim_entry.start_times = stim_start_times_this_stim;
+    stim_entry.end_times = stim_end_times_this_stim;
+    stim_entry.ntrials = sum(strcmp(stim_files_unique(i), stim_start_filename));
+    stim_data{i, 1} = stim_entry;
+    
+end
+
+%% Go through each cell and build final data matrix
+
+fs = 31250.0; % TODO: get fs from data file itself 
+pre_stim_duration = 2; %in seconds
+post_stim_duration = 2; %seconds
+pre_stim_duration_samps = pre_stim_duration*fs;
+post_stim_duration_samps = post_stim_duration*fs;
+
+n_good_units = length(good_cluster_ids);
+toedata = cell(n_good_units, 1);
+
+
+for unit_num = 1:length(good_cluster_ids)
+    
+    unit_entry = struct;
+    unit_entry.id = good_cluster_ids(unit_num);
+    stims = cell(nstims, 1);
+    
+    % Loop through each stimulus
+    for stim_num = 1:nstims
+        
+        stim_unit_entry = struct;
+        
+        this_stim_data = stim_data{stim_num, 1};
+        stim_unit_entry.name = this_stim_data.name;
+        stim_unit_entry.ntrials = this_stim_data.ntrials;
+        
+        %create trial times
+        trial_start_times = this_stim_data.start_times - pre_stim_duration_samps;
+        trial_end_times = this_stim_data.end_times + post_stim_duration_samps;
+        stim_unit_entry.stim_start_times = this_stim_data.start_times;
+        stim_unit_entry.stim_end_times = this_stim_data.end_times;
+        stim_unit_entry.trial_start_times = trial_start_times;
+        stim_unit_entry.trial_end_times = trial_end_times;
+        
+        %Go through each trial to divy up spike times relative to stim
+        %start
+        toes = cell(this_stim_data.ntrials, 1);
+        for trialnum = 1:this_stim_data.ntrials
+            trial_start = trial_start_times(trialnum);
+            trial_end = trial_end_times(trialnum);
+            
+            spiketimes_thisunit = double(good_spikes{unit_num, 2});
+            spiketimes_samps_thistrial = spiketimes_thisunit(spiketimes_thisunit >= trial_start & spiketimes_thisunit <= trial_end);
+            spiketimes_samps_thistrial_relstimonset = spiketimes_samps_thistrial - this_stim_data.start_times(trialnum);
+            spiketimes_secs_thistrial_relstimonset = double(spiketimes_samps_thistrial_relstimonset) / fs;
+            
+            toes{trialnum, 1} = spiketimes_secs_thistrial_relstimonset;
+        end
+        stim_unit_entry.toes = toes;
+        stims{stim_num, 1} = stim_unit_entry;     
+    end
+    
+    unit_entry.stims = stims;
+    toedata{unit_num, 1} = unit_entry;
+   
+end
+
+%% Format output file name
+
+data_to_save = struct();
+data_to_save.birdID = 'st1215';
+data_to_save.penetrationID = 1;
+data_to_save.siteID = 1;
+data_to_save.siteID = 1;
+data_to_save.XYpos = [0, 0];
+data_to_save.Zpos = 0;
+data_to_save.target_structure = 'NCM';
+
+data_to_save.fs = 31250;
+data_to_save.toedata = toedata;
+
+sav_date = datestr(now, 30);
+outfile_name = strcat('./st1215/good_data/pen2/more_merged/st1215_cat_P01_S01_2ndPen_moremerged_MUA', sav_date, '.mat');
+save(outfile_name, 'data_to_save');
+
+
+
diff --git a/toefile_utils/make_neurogram.m b/toefile_utils/make_neurogram.m
new file mode 100644
index 0000000..88e384f
--- /dev/null
+++ b/toefile_utils/make_neurogram.m
@@ -0,0 +1,60 @@
+% make_neurogram.m
+% Brad Theilman September 2015
+
+clear
+close all
+
+load('st1215_cat_P01_S01_3rdPen_20150927T155818.mat');
+
+%TODO:  Get sampling frequency from data file.
+fs = 31250.0;
+
+nstims = 41;
+n_units = length(toedata);
+
+psth_winsize = 150; %ms
+psth_winsize_samps = floor((psth_winsize/1000)*fs);
+
+figure();
+for stim_index = 1:nstims
+    % unit_index is which cell you wish to plot
+    
+    spike_density = zeros(n_units, 1);
+    for unit_index = 1:n_units
+        
+        % Get the data for the chosen cell
+        unit_data = toedata{unit_index, 1};
+        
+        % Collapse across all trials.  
+        trialdata = unit_data.stims{stim_index, 1};
+        ntrials = length(trialdata.toes);
+        spiketimes_across_trials = [];
+        for trial = 1:ntrials
+            spiketimes_across_trials = [spiketimes_across_trials; trialdata.toes{trial, 1}];
+        end
+        trial_length = floor(mean(trialdata.trial_end_times - trialdata.trial_start_times));
+        stim_start = floor(mean(trialdata.stim_start_times - trialdata.trial_start_times));
+        stim_end_secs = floor(mean(trialdata.stim_end_times - trialdata.stim_start_times))/fs;
+        winnum = 1;
+        for win_start = 1:psth_winsize_samps:trial_length
+            
+            win_end = min(win_start+psth_winsize_samps, trial_length) - stim_start;
+            win_start_stim = win_start - stim_start;
+            nspikes_in_win = sum(spiketimes_across_trials <= win_end/fs & spiketimes_across_trials >= win_start_stim/fs);
+            spike_density(unit_index, winnum) = nspikes_in_win / (ntrials * (psth_winsize/1000));
+            win_starts(winnum) = win_start_stim / fs;
+            winnum = winnum+1;
+        end
+        
+    end
+    max_firing_rate = max(max(spike_density))
+    subplot(7, 6, stim_index)
+    imagesc(win_starts, 1:n_units, spike_density, [0 85] );
+    colormap('gray');
+    
+    line([0, 0], [0, n_units+1], 'Color', 'red')
+    line([stim_end_secs, stim_end_secs], [0, n_units+1], 'Color', 'red');
+end
+
+            
+        
\ No newline at end of file
diff --git a/toefile_utils/make_population_psth.m b/toefile_utils/make_population_psth.m
new file mode 100644
index 0000000..a7634fd
--- /dev/null
+++ b/toefile_utils/make_population_psth.m
@@ -0,0 +1,18 @@
+% make_population_psth
+
+clear
+close all
+
+load('st1215_cat_P01_S01_3rdPen_20150927T155818.mat');
+
+%TODO:  Get sampling frequency from data file.
+fs = 31250.0;
+
+nstims = 41;
+n_units = length(toedata);
+
+psth_winsize = 150; %ms
+psth_winsize_samps = floor((psth_winsize/1000)*fs);
+
+for stim_index = 1:nstims
+    
\ No newline at end of file
diff --git a/toefile_utils/make_raster.m b/toefile_utils/make_raster.m
new file mode 100755
index 0000000..bc372fb
--- /dev/null
+++ b/toefile_utils/make_raster.m
@@ -0,0 +1,46 @@
+%% make_raster.m 
+%  Given a toefile, plot a spike raster plot for all trials of all stimuli
+%  for a single unit
+%  Brad Theilman September 2015
+
+clear
+close all
+
+%% Make a spike raster for a given cell
+
+load('st1215_cat_P01_S01_3rdPen_20150927T155818.mat');
+
+%TODO:  Get sampling frequency from data file.
+fs = 31250.0;
+
+% unit_index is which cell you wish to plot
+unit_index = 9;
+
+% Get the data for the chosen cell
+unit_data = toedata{unit_index, 1};
+nstims = length(unit_data.stims);
+
+% Plotting parameters
+spike_color = 'blue';
+stim_color = 'red';
+
+figure();
+for stimnum = 1:nstims
+    subplot(7, 6, stimnum)
+    stim_data = unit_data.stims{stimnum, 1};
+    stim_end_secs = double(stim_data.stim_end_times - stim_data.stim_start_times)/fs;
+    ntrials = stim_data.ntrials;
+    for trialnum = 1:ntrials
+        ys = [trialnum-1, trialnum];
+        if ~isempty(stim_data.toes{trialnum, 1})
+           for spikenum = 1:length( stim_data.toes{trialnum, 1})
+           line([stim_data.toes{trialnum, 1}(spikenum), stim_data.toes{trialnum, 1}(spikenum)], ys, 'Color', spike_color);
+           end
+        end
+        xlim([-2, stim_end_secs(trialnum)+2])
+        ylim([0, ntrials]);
+        line([0, 0], [0, ntrials+1], 'Color', stim_color);
+        line([stim_end_secs(trialnum), stim_end_secs(trialnum)], [0, ntrials+1], 'Color', stim_color);
+    end
+end
+
diff --git a/toefile_utils/make_raster_all_cells.m b/toefile_utils/make_raster_all_cells.m
new file mode 100755
index 0000000..5b9897d
--- /dev/null
+++ b/toefile_utils/make_raster_all_cells.m
@@ -0,0 +1,54 @@
+%% make_raster_all_cells.m 
+%  Given a toefile, plot a spike raster plot for all trials of all stimuli
+%  Does this for every unit in the file, saving the figures to png
+%  Brad Theilman September 2015
+
+clear
+close all
+
+% Make a spike raster for a given cell
+
+datafile = './pen2/more_merged/st1215_cat_P01_S01_2ndPen_moremerged_MUA20151102T092630';
+
+load(strcat(datafile, '.mat'));
+fs = 31250.0;
+
+
+for unit_index = 1:31
+    
+    disp(unit_index)
+    % Get the data for the chose cell
+    unit_data = toedata{unit_index, 1};
+    nstims = length(unit_data.stims);
+    
+    fig = gcf;
+    fig.Visible='off'
+    for stimnum = 1:nstims
+        subplot(7, 7, stimnum)
+        
+        stim_data = unit_data.stims{stimnum, 1};
+        stim_end_secs = double(stim_data.stim_end_times - stim_data.stim_start_times)/fs;
+        ntrials = stim_data.ntrials;
+        for trialnum = 1:ntrials
+            ys = [trialnum-1, trialnum];
+            if ~isempty(stim_data.toes{trialnum, 1})
+                for spikenum = 1:length(stim_data.toes{trialnum, 1})
+                    line([stim_data.toes{trialnum, 1}(spikenum), stim_data.toes{trialnum, 1}(spikenum)], ys);
+                end
+            end
+            xlim([-2, stim_end_secs(trialnum)+2])
+            ylim([0, ntrials]);
+            line([0, 0], [0, ntrials+1], 'Color', 'red');
+            line([stim_end_secs(trialnum), stim_end_secs(trialnum)], [0, ntrials+1], 'Color', 'red');
+        end
+    end
+    
+    fig.PaperUnits = 'inches';
+    fig.PaperPosition = [0 0 11 8.5];
+    fig.PaperPositionMode = 'manual';  % TODO: Save in a large resolution
+    figfilename = strcat(datafile, '_cell', num2str(unit_index));
+    print(figfilename, '-dpng', '-r600');
+    close all
+    
+    
+end
diff --git a/toefile_utils/make_raster_singlecell_singlestim.m b/toefile_utils/make_raster_singlecell_singlestim.m
new file mode 100644
index 0000000..cc3e6ac
--- /dev/null
+++ b/toefile_utils/make_raster_singlecell_singlestim.m
@@ -0,0 +1,42 @@
+%% make_raster.m
+%  Given a toefile, plot a spike raster plot for all trials of all stimuli
+%  for a single unit
+%  Brad Theilman September 2015
+
+clear
+close all
+
+%% Make a spike raster for a given cell
+
+load('st1215_cat_P01_S01_3rdPen_20150927T155818.mat');
+
+%TODO:  Get sampling frequency from data file.
+fs = 31250.0;
+
+% unit_index is which cell you wish to plot
+unit_index = 9;
+
+% Get the data for the chosen cell
+unit_data = toedata{unit_index, 1};
+nstims = length(unit_data.stims);
+
+figure();
+stimnum = 10;
+
+stim_data = unit_data.stims{stimnum, 1};
+stim_end_secs = double(stim_data.stim_end_times - stim_data.stim_start_times)/fs;
+ntrials = stim_data.ntrials;
+for trialnum = 1:ntrials
+    ys = [trialnum-1, trialnum];
+    if ~isempty(stim_data.toes{trialnum, 1})
+        for spikenum = 1:length( stim_data.toes{trialnum, 1})
+            line([stim_data.toes{trialnum, 1}(spikenum), stim_data.toes{trialnum, 1}(spikenum)], ys);
+        end
+    end
+    xlim([-2, stim_end_secs(trialnum)+2])
+    ylim([0, ntrials]);
+    line([0, 0], [0, ntrials+1], 'Color', 'red');
+    line([stim_end_secs(trialnum), stim_end_secs(trialnum)], [0, ntrials+1], 'Color', 'red');
+end
+
+
diff --git a/toefile_utils/make_raster_stim.m b/toefile_utils/make_raster_stim.m
new file mode 100755
index 0000000..5f570c0
--- /dev/null
+++ b/toefile_utils/make_raster_stim.m
@@ -0,0 +1,39 @@
+%% make_raster.m 
+%  Given a toefile, plot a spike raster plot for all trials of one stimulus
+%  for all units
+%  Brad Theilman September 2015
+
+
+clear
+close all
+
+load('st1215_cat_P01_S01_2ndPen_fixalignment_20150924T141526.mat');
+
+fs = 31250.0; % TODO: get this from the data file
+
+stimnum = 49;
+
+nstims = 49; % TODO: get this from the data file
+
+figure();
+for unit_index= 1:31
+    subplot(8, 4, unit_index)
+
+    unit_data = toedata{unit_index, 1};
+    stim_data = unit_data.stims{stimnum, 1};
+    stim_end_secs = double(stim_data.stim_end_times - stim_data.stim_start_times)/fs;
+    ntrials = stim_data.ntrials;
+    for trialnum = 1:ntrials
+        ys = [0 + trialnum, 1+trialnum];
+        if ~isempty(stim_data.toes{trialnum, 1})
+           for spikenum = 1:length( stim_data.toes{trialnum, 1})
+           line([stim_data.toes{trialnum, 1}(spikenum), stim_data.toes{trialnum, 1}(spikenum)], ys);
+           end
+        end
+        xlim([-2, stim_end_secs(trialnum)+2])
+        ylim([0, ntrials+1]);
+        line([0, 0], [0, ntrials+1], 'Color', 'red');
+        line([stim_end_secs(trialnum), stim_end_secs(trialnum)], [0, ntrials+1], 'Color', 'red');
+    end
+end
+