if exist('results/H1_param.mat', 'file')
    load('results/H1_param.mat');
    return
end

%% H1 neuron's response to velocity input
% deRuyter van Steveninck 1997 data via Zhuo Wang
% it's originally from Dayan & Abbott textbook
%
% Data format
%   rho: binary spike train
%   stim: white noise stimulus
load('data/h1whitenoise.mat');

%%
Fs = 1000; % sampling rate in Hz
dt = 1/Fs; % bin size in seconds
fprintf('Mean firing rate: %f\n', mean(rho) * Fs);

%% Compute STA
lag = 120;
X = buffer(stim, lag, lag-1);
lags = -fliplr(1:lag);
STA = X * rho;

%% plot STA
fig0 = figure(2010); clf;
plot(lags, STA);
grid on;
title('Spike triggered average (STA)');
xlabel('lags (ms)');

%% project on to STA
Z = X' * STA;

%% pseudo-density of projection
plot(Z, rho + 0.2*randn(size(rho)), '.')

%% let's smooth the spike train over time
slen = 25;
y2 = filter(ones(slen, 1), 1, rho); % this is causal filtering
z2 = filter(ones(slen, 1)/slen, 1, Z);
z2 = zscore(z2);

fig1 = figure(2910); clf;
plot(z2, (y2 + 0.1*randn(size(y2)))/slen*Fs, '.');
ylabel('firing rate (spikes/sec)');
xlabel('z-scored stimulus projected onto STA');

%% (marginal) stimulus density
fig2 = figure(774); clf; hold all;
xr = linspace(-5,5);
histogram(z2, 'Normalization', 'pdf')
plot(xr, ksdensity(z2, xr), 'LineWidth', 2, 'Color', [0, 0, 0.05]);
z2density = normpdf(xr, 0, std(z2));
plot(xr, z2density, 'LineWidth', 2, 'Color', 'r');
title('stimulus distribution')
ylabel('probability density');
legend('histogram', 'kernel density', 'normal fit');

%% histogram based nonlinearity
figure(fig1); hold on;
edges = [-Inf, quantile(z2, linspace(0, 1, 100)), Inf];
my = zeros(numel(edges) - 1, 1);
for k = 1:numel(edges)-1
    bidx = edges(k) <= z2 & edges(k+1) > z2;
    my(k) = mean(y2(bidx));
end
hold all; plot(edges(2:end-2), my(1:end-2)/slen*Fs, 'k.-', 'LineWidth', 5)

%% equally binned estimator with error bars
edges = linspace(-5, 5, 101);
my = zeros(numel(edges)-1, 1);
sy = zeros(numel(edges)-1, 1);
for k = 1:numel(edges)-1
    bidx = edges(k) <= z2 & edges(k+1) > z2;
    my(k) = mean(y2(bidx));
    sy(k) = std(y2(bidx)) / sqrt(sum(bidx));
end

%%
hold on
plot(edges(2:end-2), (my(1:end-2) + 2*sy(1:end-2))/slen*Fs, 'LineWidth', 1, 'Color', 0.7 * [1,1,1]);
plot(edges(2:end-2), (my(1:end-2) - 2*sy(1:end-2))/slen*Fs, 'LineWidth', 1, 'Color', 0.7 * [1,1,1]);
plot(edges(2:end-2), my(1:end-2)/slen*Fs, 'k-', 'LineWidth', 3);


%% Let's do GP-Poisson regression
% Requires GPML toolbox: http://www.gaussianprocess.org/gpml/code/matlab/doc/

%%
mf = @meanLinear; 
lik = {@likPoisson, 'exp'};
cov = @covSEiso;
inf = @infFITC_Laplace;

%hyp = struct();
%hyp.mean = [0.7]; hyp.cov = [-0.1, -0.1];
%idx = randperm(numel(z2));
%idx = idx(1:4000);
%tic
%hyp2 = minimize(hyp, @gp, -100, inf, mf, cov, lik, z2(idx), y2(idx));
%toc

% pre-optimized for you:
hyp3.mean = 0.7849;
hyp3.cov = [-0.0886 -0.7339];

%%
idx = randperm(numel(z2));
idx = idx(1:ceil(0.02 * numel(z2)));

tic
[m1, s1, fm1, fs1] = gp(hyp3, inf, mf, cov, lik, z2(idx), y2(idx), xr(:));
toc

%%
ph = plot(xr, exp(fm1)/slen*Fs, '-r', 'LineWidth', 2);
ph1 = plot(xr, exp(fm1+2*fs1)/slen*Fs, '-r', 'LineWidth', 1);
ph2 = plot(xr, exp(fm1-2*fs1)/slen*Fs, '-r', 'LineWidth', 1);
xlim([-3, 3]);

%% use GP regression nonlinearity to get marginal response distribution
fr = exp(fm1')/slen*Fs;

%%
fig3 = figure(4381); clf;
frdensity = z2density ./ diff([fr, fr(end)+1]) .* (xr(2) - xr(1));
disp(sum(frdensity(1:end) .* diff([fr, fr(end)+1])))
plot(fr, frdensity);
yr = linspace(0, 500);
plot(yr, ksdensity(y2/slen*Fs, yr, 'Bandwidth', 20));
title('marginal response rate distribution'); xlabel('firing rate'); ylabel('probability density')

%%
sr = 0:(max(y2)+1);
pfr = zeros(numel(sr), numel(xr));
for kr = 1:numel(xr)
    pfr(:, kr) = poisspdf(sr, fr(kr)*slen/Fs);
end
frdensity(end) = 0;
y2h = histcounts(y2, [sr, Inf], 'Normalization', 'pdf');
bar(sr, [y2h', pfr * (frdensity .* diff([fr, fr(end)+1]))', poisspdf(sr, mean(y2))']);
legend('Data', 'LNPoisson', 'homogeneous Poisson')
title('marginal response rate distribution'); 
xlabel('spike count');
ylabel('probability');

%% save params
save('H1_param.mat', 'y2', 'fr', 'slen', 'Fs', 'z2', 'xr', 'sr', 'fs1', 'fm1');