close all
files = dir('../unit*.mat');
for file = files'
	a = load(strcat('../', file.name));
	spikes = a.spikes;
	angles = a.angles;
	figure()
	for k = 1:size(spikes, 1)
		subplot(3, 4, k)
		spikeraster(spikes(k,:), -0.2, 0.6);
	end
end


%% tuning curves:
close all
cosine = @(p,xdata)0.5*p(1).*(1.0-cos(2.0*pi*(xdata/180.0-p(2))));
files = dir('../unit*.mat');
phases = zeros(length(files), 1);
figure()
for j = 1:length(files)
	file = files(j);
	a = load(strcat('../', file.name));
	spikes = a.spikes;
	angles = a.angles;
	rates = zeros(size(spikes, 1), 1);
	for k = 1:size(spikes, 1)
		r = firingrate(spikes(k,:), 0.0, 0.2);
		rates(k) = r;
	end
	[mr, maxi] = max(rates);
	p0 = [mr, angles(maxi)/180.0-0.5];
	%p = p0;
	p = lsqcurvefit(cosine, p0, angles, rates');
	phase = p(2)*180.0 + 90.0;
    if phase > 180.0
        phase = phase - 180.0;
    end
    phases(j) = phase;
	subplot(2, 3, j);
	plot(angles, rates, 'b');
	hold on;
	plot(angles, cosine(p, angles), 'r');
	hold off;
	xlim([0.0 180.0])
	ylim([0.0 50.0])
	title(sprintf('unit %d', j))
end


%% read out:
a = load('../population04.mat');
spikes = a.spikes;
angle = a.angle;
unitphases = a.phases*180.0 + 90.0;
unitphases(unitphases>180.0) = unitphases(unitphases>180.0) - 180.0;
figure();
subplot(1, 3, 1);
angleestimates1 = zeros(size(spikes, 2), 1);
angleestimates2 = zeros(size(spikes, 2), 1);
for j = 1:size(spikes, 2)
    rates = zeros(size(spikes, 1), 1);
    for k = 1:size(spikes, 1)
        r = firingrate(spikes(k, j), 0.0, 0.2);
        rates(k) = r;
    end
    [x, inx] = sort(phases);
    plot(phases(inx), rates(inx), '-o');
    hold on;
    angleestimates1(j) = popvecangle(phases, rates);
    [m, i] = max(rates);
    angleestimates2(j) = phases(i);
end
hold off;
subplot(1, 3, 2);
hist(angleestimates1);
subplot(1, 3, 3);
hist(angleestimates2);
angle
mean(angleestimates1)
mean(angleestimates2)


%% read out robustness:
files = dir('../population*.mat');
angles = zeros(length(files), 1);
e1m = zeros(length(files), 1);
e1s = zeros(length(files), 1);
e2m = zeros(length(files), 1);
e2s = zeros(length(files), 1);
for i = 1:length(files)
	file = files(i);
	a = load(strcat('../', file.name));
    spikes = a.spikes;
    angle = a.angle;
    angleestimates1 = zeros(size(spikes, 2), 1);
    angleestimates2 = zeros(size(spikes, 2), 1);
    for j = 1:size(spikes, 2)
        rates = zeros(size(spikes, 1), 1);
        for k = 1:size(spikes, 1)
            r = firingrate(spikes(k, j), 0.0, 0.2);
            rates(k) = r;
        end
        angleestimates1(j) = popvecangle(phases, rates);
        [m, inx] = max(rates);
        angleestimates2(j) = phases(inx);
    end
    angles(i) = angle;
    e1m(i) = mean(angleestimates1);
    e1s(i) = std(angleestimates1);
    e2m(i) = mean(angleestimates2);
    e2s(i) = std(angleestimates2);
end
figure();
subplot(1, 2, 1);
scatter(angles, e1m);
subplot(1, 2, 2);
scatter(angles, e2m);