function [OutAx,x,sii,si_LR] = ComputeSelectivity(ev,ts,cnd,varargin)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
iod = @utils.inputordefault;
pre = iod('pre',0,varargin);
post = iod('post',1,varargin);
binsz = iod('binsz',0.01,varargin);
krn = iod('krn',0.1,varargin);
clr = iod('color','b',varargin);
corner=iod('corner',[0.1 0.1],varargin);
ax_width=iod('ax_width',0.55,varargin);
one_height=iod('one_height',0.248,varargin);
total_height=iod('total_height',0.8,varargin);
xlabel_str=iod('xlabel_str','Time from Subgoal / s',varargin);
title_str = iod('title_str','Subgoal Only',varargin);

if iscell(cnd)
    %cnd_nan = cellfun(@(x)isnan(x), cnd);
    cnd_nan = cellfun(@(x)strcmp(x,'null'), cnd);
    %cnd(cnd_nan) = 'NaN';
    cnd = categorical(cnd);
    n_cnd = categories(cnd);
else
    cnd = categorical(cnd);
    n_cnd = categories(cnd);
end


if isscalar(krn)
    dx=ceil(5*krn);
    kx=-dx:binsz:dx;
    krn=normpdf(kx,0, krn);
	if isempty(find(kx==0, 1))
		error('Your binsz needs to divide 1 second into interger # of bins');
	end
    krn(kx<0)=0;
    krn=krn/sum(krn);
end

[Y,x]=stats.spike_filter(ev,ts,krn,'pre',pre,'post',post,'kernel_bin_size',binsz);

for ci=1:numel(n_cnd)
    sampz=sum(cnd==n_cnd(1));
    ref=cnd==n_cnd(ci);
    
    y=Y(ref,:);
    
    ymn(ci,:) = nanmean(y,1);
    yst(ci,:)= stats.nanstderr(y,1);
end

dump = sort(ymn);
%sii = (dump(end,:)-(mean(dump(1:end-1,:))/(size(ymn,1)-1)))./mean(ymn);
sii = (dump(end,:)-(mean(dump(1:end-1,:))))./sum(ymn);

dump = zeros(2,length(x));
dump(1,:) = mean(ymn(1:2:end,:));%L
dump(2,:) = mean(ymn(2:2:end,:));%R
dump = sort(dump);
si_LR = (dump(2,:)-dump(1,:))./sum(dump);

ax_si_LR = draw.jaxes([corner,ax_width,one_height]);
plot(ax_si_LR,x,si_LR,'LineWidth',3);

line([0,0],[0,max(si_LR)+0.5])
set(gca,'YLim',[0 max(si_LR)+0.5])
set(gca,'XLim',[-pre post]);
title([title_str,'Preferred L/R Selectivity'])
xlabel(xlabel_str)
ylabel('SI / a.u.')

ax_sii = draw.jaxes([corner(1),corner(2)+total_height-one_height,ax_width,one_height]);
plot(ax_sii,x,sii,'LineWidth',3);
line([0,0],[0,max(sii)+0.5])
set(gca,'YLim',[0 max(sii)+0.5])
set(gca,'XLim',[-pre post]);
title([title_str,'Preferred Direction Selectivity'])
xlabel(xlabel_str)
ylabel('SI / a.u.')

OutAx(1) = ax_si_LR;
OutAx(2) = ax_sii;
end