% 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 a bounded distribution between 0 and 1 with mean of
%0.01 and 100 bins. These distributions will be mapped to synaptic weight distributions
%between 0 and threshold = 20 mV with a mean of 2mV and bin-size of 0.2 mV
%
%The algorithm essentially uses a random walk between 0 and 1 with an
%attractive term that forces the walk to return to the mean. Each time a
%bin is visited, a histogram is built up by adding 1 to that bin. 

%NOTE: The .mat file FigS10_randdist.mat consisting of the original
%distributions used for FigS10 has also been made available and can be 
%used to reproduce the figure exactly. 

clear variables;

x1=1;         %end of interval
dx=0.01;      %discretization step
x0=dx;        %beginning of interval
xm=0.1;       %mean of distribution

x=x0:dx:x1;   %discretization vector
s=length(x);  
di=zeros(s,1);%counter to build up histogram

nmax=10000;   %maximum no. of iterations to try and converge to mean
nmin=1000;    %minimum no. of iterations
tol=10^-6;    %tolerance for difference from the desired mean

tic

%Generate grid of ri and rd for multiple distributions 
ds1 = 2;      
ds2 = 1;
ratiomax = 100;
ri = 1:ds1:s/2;
for jj = 1:length(ri)
 rd = 1:ds2:ri(jj)*(1-s/2)/(ratiomax*(ri(jj)-s/2));
end;

%Initialize
dist=zeros(s,1);
cvgmat = [];
rvals = [];

%Run script
ctr=1;
while ctr<41;
    for ii = 1:length(ri)
        for jj = 1:length(rd)
            if(ri(ii)/rd(jj)>=1 && ri(ii)/rd(jj)<=ratiomax)
                r=randi(s);
                di(r)=di(r)+1;
                for i=1:nmax
                    if x*di/sum(di)>xm
                        rmin=max(1,r-ri(ii)-rd(jj));
                        rmax=min(s,r+ri(ii)-rd(jj));
                    else
                        rmin=max(1,r-ri(ii)+rd(jj));
                        rmax=min(s,r+ri(ii)+rd(jj));
                    end
                    r=rmin+randi(rmax-rmin+1)-1;
                    di(r)=di(r)+1;
                    if i>nmin && abs(x*di/sum(di)-xm)<tol
                        break
                    end
                end
                i=i+1;
                if (i>nmax)
                    converged=0;
                else
                    converged=1;
                    dist = horzcat(dist,di/sum(di));
                end
                cvgmat=[cvgmat;converged];
                racc = [ri(ii) rd(jj) converged];
                rvals = vertcat(rvals,racc);
            end
            di = zeros(s,1);
        end
    end
    ctr = ctr+1;
end
dist = dist(:,2:end);

toc

%save('FigS10_randdist.mat','dist','rvals','x');

