% 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 to generate synaptic distributions that match a given drift, where 
%drift=sum(wts.*probs)*fr
function [wts qwts xopt]=match_distribution_1_var(opt_drift,wts,dist_type,fr);

switch dist_type
    case 'gauss'
        disp('Matching Gaussian distribution, searching over mean=stdev')
        wts_mean=wts(floor(length(wts)/5));
        [xopt fval exitflag] = fminsearch(@(x) drift_err_gauss(x,wts,opt_drift,fr),[wts_mean],optimset('TolX',1e-12,'TolFun',1e-10,'MaxFunEvals',1e10,'MaxIter',1e10));
        qwts = gauss_dist(0.6579*10^-3,xopt(1),wts);
    case 'exp'
        disp('Matching Exponential distribution, searching over mean of weight distribution')
        wts_mean=wts(floor(length(wts)/5));
        [xopt fval exitflag] = fminsearch(@(x) drift_err_exp(x,wts,opt_drift,fr),[wts_mean],optimset('TolX',1e-12,'TolFun',1e-10,'MaxFunEvals',1e10,'MaxIter',1e10));
        qwts = exp_dist(xopt(1),wts);
        xopt(1)=xopt(1)^2;
    case 'hyp'
        disp('Matching Hyperbolic distribution, searching over weight parameter used in weight distribution')
        wts_mean=wts(floor(length(wts)/5));
        [xopt fval exitflag] = fminsearch(@(x) drift_err_hyp(x,wts,opt_drift,fr),[wts_mean],optimset('TolX',1e-12,'TolFun',1e-10,'MaxFunEvals',1e10,'MaxIter',1e10));
        qwts = hyp_dist(xopt(1),wts);
        xopt(1)=xopt(1)^2;
    case 'logn'
        disp('Matching Lognormal distribution, searching over mu = sigma')
        wts_std=100;%log(opt_drift/fr);
        [xopt fval exitflag] = fminsearch(@(x) drift_err_logn(x,wts,opt_drift,fr),[wts_std],optimset('TolFun',1e-10,'MaxFunEvals',1e10,'MaxIter',1e10));
        std=1+xopt(1)*10^-3;
        m=log(0.6579*10^-3);
       qwts = logn_dist(m,std,wts);
end
if (exitflag ~= 1) 
    disp('Exit flag is not equal to 1 - non-optimal solution');
end;

%%%error functions
function err = drift_err_gauss(x,wts,opt_drift,fr);
    m=0.6579*10^-3;
    std=x(1);
    qwts=gauss_dist(m,std,wts);
    [drift diff]=calc_drift_diff(wts,qwts,fr);
    err = ((drift-opt_drift)/opt_drift)^2;

function err = drift_err_logn(x,wts,opt_drift,fr);
    std=1+x(1)*10^-3;
    m=log(0.6579*10^-3);%forces the lognormal median to be the same as the exp distibution
    qwts=logn_dist(m,std,wts);
    [drift diff]=calc_drift_diff(wts,qwts,fr);
    err = ((drift-opt_drift)/opt_drift)^2;

function err = drift_err_exp(x,wts,opt_drift,fr);
m=x(1);
qwts=exp_dist(m,wts);
[drift diff]=calc_drift_diff(wts,qwts,fr);
err = ((drift-opt_drift)/opt_drift)^2;

function err = drift_err_hyp(x,wts,opt_drift,fr);
m=x(1);
qwts=hyp_dist(m,wts);
[drift diff]=calc_drift_diff(wts,qwts,fr);
err = ((drift-opt_drift)/opt_drift)^2;


function [driftterm diffterm] = calc_drift_diff(wts,qwts,fr);
    driftterm = sum(wts.*qwts)*fr;
    diffterm = sum(wts.*wts.*qwts)*fr;

    
function qwts = gauss_dist(m,std,wts);
    qwts = exp(-(wts-m).^2/(2*std^2));
    qwts = qwts/sum(qwts);
    
function qwts = exp_dist(m,wts);
    m2=m^2;
    qwts = exp(-(wts/m2));
    qwts = qwts/sum(qwts);

function qwts = hyp_dist(m,wts);
    m=m*m;
    pwr=1; 
    qwts = 1./((1+(wts-wts(1))/m).^pwr);
    qwts = qwts/sum(qwts);
    
function qwts = logn_dist(m,std,wts);
    %eliminate zero wts
    zero_index= wts==0;
    wts(zero_index)=1;
    qwts = (1./wts).*exp(-(-m+log(wts)).^2/(2*std^2));
    qwts(zero_index)=0;
    qwts = qwts/sum(qwts);
