% 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

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

function full_synaptic_matrix=generate_synaptic_matrix2_testsyn2(qdw,v,f,syn,flag_cond_syn,flag_poisson,voltage_fraction);
%Function to generate the transition matrix for a given distribution of synaptic weights 
%matrix, where multiple spikes are treated as powers of single-spike
%matrices. Please refer to METHODS: NUMERICAL SOLUTIONS in the paper.
%This matrix acts on p(v,t). 
%
%Inputs: qdw=synaptic weight distribution, v=geometric voltage bins,
%f=(fr*dt) is the mean number of input spikes in time dt, flag_cond_syn=1
%(for conductance-based synapses) and =0(current-based synapses), 
%flag_poisson=1(for Poisson input spikes), 
%voltage_fraction=measure of residual current after spike.

%%%%NOTE : Not fully implemented for non-Poisson inputs.

if nargin<7
    voltage_fraction=0;
end;
if nargin<6
    voltage_fraction=0;
    flag_poisson=1;
end;
if nargin<5
    voltage_fraction=0;
    flag_poisson=1;
    flag_cond_syn=0;
end;

nv=length(v);
first_plus_bin=min(v(v>0));
padding=2;                    %how far to extend binning beyond threshold
th=max(v)*padding;      
paddingbinsize=first_plus_bin;
v=[v v(end)+paddingbinsize:paddingbinsize:th];
synaptic_mat=eye(length(v),length(v));
vleft=v;
vright=[v(2:end) v(end)+1];
zero_index=find(vleft<=0 & vright>=0);
dvs=vright-vleft;

for ii=1:length(qdw(1,:))       %for each discrete weight value
    w=qdw(1,ii);
    if (w~=0)
        vr=qdw(2,ii);
        qprob=qdw(3,ii);
        f_prob=1;
        if (flag_cond_syn==1)
            new_vs=v+w*(vr-v); %new voltages as a result of synaptic activation (conductance-based)
        else
            new_vs=v+w;        %new voltages as a result of synaptic activation (current-based)
        end
        vr_min=v(1);
        for binno=1:length(new_vs)-1        %loop over all the new voltage bins
            leftedge=new_vs(binno);
            rightedge=new_vs(binno+1);
            newbinwidth=rightedge-leftedge;
            if (rightedge<th)               %no threshold crossing
                if (leftedge>=vr_min)       %no reversal potential crossing
                    withinbins=find(vleft<=leftedge & vright>=rightedge);   %shifted bin falls entirely within an existing bin
                    fullbins=find(vleft>=leftedge & vright<=rightedge);     %shifted bin entirely encompasses an existing bin
                    fullbins=setdiff(fullbins,withinbins);
                    leftoverhang=find(vleft<leftedge & vright>leftedge & vright<rightedge);
                    rightoverhang=find(vleft<rightedge & vright>rightedge & vleft>leftedge);
                    if ~isempty(withinbins) 
                        synaptic_mat(binno,withinbins)=synaptic_mat(binno,withinbins)+qprob*f_prob; 
                    end;
                    if ~isempty(fullbins) 
                        synaptic_mat(binno,fullbins)=synaptic_mat(binno,fullbins)+qprob*f_prob*dvs(fullbins)/newbinwidth; 
                    end;
                    if ~isempty(leftoverhang) 
                        synaptic_mat(binno,leftoverhang)=synaptic_mat(binno,leftoverhang)+qprob*f_prob*(vright(leftoverhang)-leftedge)/newbinwidth; 
                    end;
                    if ~isempty(rightoverhang) 
                        synaptic_mat(binno,rightoverhang)=synaptic_mat(binno,rightoverhang)+qprob*f_prob*(rightedge-vleft(rightoverhang))/newbinwidth; 
                    end;
                else                       %crossing the left (negative) boundary
                    if (rightedge<=vr_min) %entire bin crosses left boundary
                        synaptic_mat(binno,1)=synaptic_mat(binno,1)+qprob*f_prob;
                    else
                        fullbins=find(vright<=rightedge);
                        rightoverhang=find(vleft<rightedge & vright>rightedge & vleft>leftedge);
                        rightoverhang=setdiff(rightoverhang,fullbins);
                        if ~isempty(fullbins) 
                            synaptic_mat(binno,fullbins)=synaptic_mat(binno,fullbins)+qprob*f_prob*dvs(fullbins)/newbinwidth; 
                        end;
                        if ~isempty(rightoverhang) 
                            synaptic_mat(binno,rightoverhang)=synaptic_mat(binno,rightoverhang)+qprob*f_prob*(rightedge-vleft(rightoverhang))/newbinwidth; 
                        end;
                        synaptic_mat(binno,1)=synaptic_mat(binno,1)+qprob*f_prob*(max(vleft(1)-leftedge,0))/newbinwidth;
                    end;
                end
            else                          %crossing last bin, so move everything into that bin
                synaptic_mat(binno,end)=synaptic_mat(binno,end)+qprob*f_prob;   %move everything into last bin
            end
            synaptic_mat(binno,binno)=synaptic_mat(binno,binno)-qprob*f_prob;
        end;
    end;
end

%Multiply matrices to account for multiple spikes
if (syn==0) 
    full_synaptic_matrix=eye(nv,nv);
else
    combined_synaptic_matrix=zeros(size(synaptic_mat));
    if (flag_poisson==1)
        f_prob=f^syn*exp(-f)/factorial(syn);
    else
        f_prob=f;
    end;
    temp_synaptic_mat=synaptic_mat^syn;
    combined_synaptic_matrix=combined_synaptic_matrix+f_prob*temp_synaptic_mat;
    full_synaptic_matrix=combined_synaptic_matrix(1:nv,1:nv);
    %calculate firing rate as sum of everything beyond the original
    %threshold
    full_synaptic_matrix(:,end)=sum(combined_synaptic_matrix(1:nv,nv:end)');
    %%%reinject probability beyond threshold into bins near rest, depending
    %%%on voltage_frac
    if (voltage_fraction==0)
        full_synaptic_matrix(:,1)=full_synaptic_matrix(:,end);    
    else
        %shrink bins beyond the original threshold, based on voltage frac
        bins_past_thresh=(v(nv:end)-v(nv))*voltage_fraction;
        %rebin the past-threshold bins into the same bin widths as those
        %near resting potential
        for ii=1:nv
            reinject_probs=rebin_uniform_vector(bins_past_thresh,v(zero_index:nv),combined_synaptic_matrix(ii,nv:end));
            if (sum(combined_synaptic_matrix(ii,nv:end))>0)
                length(reinject_probs);
            end;
            %pad synaptic matrix with zeros if inhibitory synapses are
            %present
            if(v(1) < 0)                    
                pad_matrix = zeros(1,(nv/2-1));
            else
                pad_matrix=[];
            end
            full_synaptic_matrix(ii,:)=full_synaptic_matrix(ii,:)+[pad_matrix reinject_probs];
        end
    end
end;

%Diagonal should contain all non-transitioned probability
full_synaptic_matrix=full_synaptic_matrix-diag(diag(full_synaptic_matrix));
diagIdx = 1:nv+1:nv^2;
missing = 1-sum(full_synaptic_matrix(:,1:end-1)');
full_synaptic_matrix(diagIdx)=full_synaptic_matrix(diagIdx)+missing;
