% This file is licensed under the Allen Institute Terms of Use, which can be found at:
%http://www.alleninstitute.org/Media/policies/terms_of_use_content.html

%------------------------------------------------------------------------------------

%Script to generate data for Fig. S9 panels b) and d). 
%This script convolves synaptic weight distributions q(w) matched for drift
%diffusion with a Poisson process to obtain qp(w) in order to 
%generate the tail weight numbers Qn. Please refer to eqns. (6) and (7)
%in the paper for definition of tail-weight numbers. The Poisson process
%is truncated at k = kmax. This script generates the first 150 tail weight 
%numbers for these matched synaptic weight distributions.

clear variables;

th=20*10^-3;  %V Threshold

dv=0.2*10^-3; %V defines max step in potential
v0=(0:dv:th); %V defines the grid for the membrane potential

%initial probability distribution; start at rest
p0=zeros(size(v0));
p0(1)=1;      

match_drift_only=0;          %match distribution for drift and diffusion
match_4_distributions_1_var;

qdw_hyp_dd=[0 wts_hyp; 1 ones(size(wts_hyp));0 qwts_hyp]; 
qdw_gauss_dd=[0 wts_gauss;1 ones(size(wts_gauss));0 qwts_gauss];
qdw_exp_dd=[0 wts_exp; 1 ones(size(wts_exp)); 0 qwts_exp]; 
qdw_logn_dd=[0 wts_logn; 1 ones(size(wts_logn)); 0 qwts_logn]; 

wts_bi = 0:0.1*10^-3:th;

bi_p = zeros(size(wts_bi));
bi_p(wts_bi==0.0005)=0.9967;
bi_p(wts_bi==0.015)=0.0033;

del_p = zeros(size(wts_bi));
del_p(19) = 1;

qdw_del_dd=[wts_bi;ones(size(wts_bi));del_p];
qdw_bidel_dd= [wts_bi;ones(size(wts_bi));bi_p];

kmax = 10;
lambda = 1;

%Delta function distribution convolve with Poisson 
qdw_del_probs_prime = exp(-1)*[qdw_del_dd(3,:) zeros(1,200*(kmax-1))];
    qdw_dist_fin_probs_temp = qdw_del_dd(3,:);
    for k = 2:10
        qdw_dist_fin_probs_temp = conv(qdw_dist_fin_probs_temp,qdw_del_dd(3,:));
        qdw_dist_fin_probs_padded = [qdw_dist_fin_probs_temp zeros(1,200*(kmax-k))];
        qdw_del_probs_prime = qdw_del_probs_prime+((lambda)^k*exp(-lambda)/factorial(k))*qdw_dist_fin_probs_padded;
    end
qdw_del_probs_prime(1)=exp(-1);

%Bimodal distribution convolve with Poisson 
qdw_bidel_probs_prime = exp(-1)*[qdw_bidel_dd(3,:) zeros(1,200*(kmax-1))];
    qdw_dist_fin_probs_temp = qdw_bidel_dd(3,:);
    for k = 2:10
        qdw_dist_fin_probs_temp = conv(qdw_dist_fin_probs_temp,qdw_bidel_dd(3,:));
        qdw_dist_fin_probs_padded = [qdw_dist_fin_probs_temp zeros(1,200*(kmax-k))];
        qdw_bidel_probs_prime = qdw_bidel_probs_prime+((lambda)^k*exp(-lambda)/factorial(k))*qdw_dist_fin_probs_padded;
    end
qdw_bidel_probs_prime(1)=exp(-1);

%Gaussian distribution convolve with Poisson 
qdw_gauss_probs_prime = exp(-1)*[qdw_gauss_dd(3,:) zeros(1,100*(kmax-1))];
    qdw_dist_fin_probs_temp = qdw_gauss_dd(3,:);
    for k = 2:10
        qdw_dist_fin_probs_temp = conv(qdw_dist_fin_probs_temp,qdw_gauss_dd(3,:));
        qdw_dist_fin_probs_padded = [qdw_dist_fin_probs_temp zeros(1,100*(kmax-k))];
        qdw_gauss_probs_prime = qdw_gauss_probs_prime+((lambda)^k*exp(-lambda)/factorial(k))*qdw_dist_fin_probs_padded;
    end
qdw_gauss_probs_prime(1)=exp(-1);

%Exponential distribution convolve with Poisson 
qdw_exp_probs_prime = exp(-1)*[qdw_exp_dd(3,:) zeros(1,100*(kmax-1))];
    qdw_dist_fin_probs_temp = qdw_exp_dd(3,:);
    for k = 2:10
        qdw_dist_fin_probs_temp = conv(qdw_dist_fin_probs_temp,qdw_exp_dd(3,:));
        qdw_dist_fin_probs_padded = [qdw_dist_fin_probs_temp zeros(1,100*(kmax-k))];
        qdw_exp_probs_prime = qdw_exp_probs_prime+((lambda)^k*exp(-lambda)/factorial(k))*qdw_dist_fin_probs_padded;
    end
qdw_exp_probs_prime(1)=exp(-1);

%Lognormal distribution convolve with Poisson  
qdw_logn_probs_prime = exp(-1)*[qdw_logn_dd(3,:) zeros(1,100*(kmax-1))];
    qdw_dist_fin_probs_temp = qdw_logn_dd(3,:);
    for k = 2:10
        qdw_dist_fin_probs_temp = conv(qdw_dist_fin_probs_temp,qdw_logn_dd(3,:));
        qdw_dist_fin_probs_padded = [qdw_dist_fin_probs_temp zeros(1,100*(kmax-k))];
        qdw_logn_probs_prime = qdw_logn_probs_prime+((lambda)^k*exp(-lambda)/factorial(k))*qdw_dist_fin_probs_padded;
    end
qdw_logn_probs_prime(1)=exp(-1);

%Power law distribution convolve with Poisson
qdw_hyp_probs_prime = exp(-1)*[qdw_hyp_dd(3,:) zeros(1,100*(kmax-1))];
    qdw_dist_fin_probs_temp = qdw_hyp_dd(3,:);
    for k = 2:10
        qdw_dist_fin_probs_temp = conv(qdw_dist_fin_probs_temp,qdw_hyp_dd(3,:));
        qdw_dist_fin_probs_padded = [qdw_dist_fin_probs_temp zeros(1,100*(kmax-k))];
        qdw_hyp_probs_prime = qdw_hyp_probs_prime+((lambda)^k*exp(-lambda)/factorial(k))*qdw_dist_fin_probs_padded;
    end
qdw_hyp_probs_prime(1)=exp(-1);

%convolve q(w) with initial probability distribution
c_del = conv(p0,qdw_del_probs_prime);
c_gauss = conv(p0,qdw_gauss_probs_prime);
c_exp = conv(p0,qdw_exp_probs_prime);
c_logn = conv(p0,qdw_logn_probs_prime);
c_hyp = conv(p0,qdw_hyp_probs_prime);
c_bidel = conv(p0,qdw_bidel_probs_prime);

%threshold indices
th1 = length(wts_hyp)+1;%101;
th_bi = length(wts_bi);%201;

%calculate tail weight numbers Qn
r_del_dd = [sum(c_del(th_bi:end))];
r_gauss_dd = [sum(c_gauss(th1:end))];
r_exp_dd = [sum(c_exp(th1:end))];
r_logn_dd = [sum(c_logn(th1:end))];
r_hyp_dd = [sum(c_hyp(th1:end))];
r_bidel_dd = [sum(c_bidel(th_bi:end))];

for i = 1:150
    c_del = conv(c_del,qdw_del_probs_prime);
    r_del_dd = [r_del_dd sum(c_del(th_bi:end))];
    c_gauss = conv(c_gauss,qdw_gauss_probs_prime);
    r_gauss_dd = [r_gauss_dd sum(c_gauss(th1:end))];
    c_exp = conv(c_exp,qdw_exp_probs_prime);
    r_exp_dd = [r_exp_dd sum(c_exp(th1:end))];
    c_logn = conv(c_logn,qdw_logn_probs_prime);
    r_logn_dd = [r_logn_dd sum(c_logn(th1:end))];
    c_hyp = conv(c_hyp,qdw_hyp_probs_prime);
    r_hyp_dd = [r_hyp_dd sum(c_hyp(th1:end))];
    c_bidel = conv(c_bidel,qdw_bidel_probs_prime);
    r_bidel_dd = [r_bidel_dd sum(c_bidel(th_bi:end))];
end

save('drdiffmatched_resp_nos_poisson.mat','r_del_dd','r_gauss_dd','r_exp_dd','r_logn_dd','r_hyp_dd','r_bidel_dd','qdw_del_dd','qdw_gauss_dd','qdw_exp_dd','qdw_logn_dd','qdw_hyp_dd','qdw_bidel_dd');
