% 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 matched exponential, hyperbolic, Gaussian, and
%lognormal distributions, given a weight vector, a firing rate, and a drift term
%For the Gaussian and lognormal, the mean is set equal to the stdev in
%order to fit only one parameter

%Targets
dv=0.2*10^-3;   %V
th=20*10^-3;    %V
wts=dv:dv:th;   %Weight vector in V
drift=1;
fr = 1000;
fr_exp=fr;
fr_gauss=fr;
fr_hyp=fr;
fr_logn=fr;

if (match_drift_only==1)    %Match only the drift term, with all distributions having the same firing rate (and hence the same mean)
    
    %Exponential
    [wts_exp qwts_exp xopt_exp]=match_distribution_1_var(drift,wts,'exp',fr);
    disp(['Parameters for exponential: weight parameter = ',num2str(xopt_exp(1)),' V, firing rate = ',num2str(fr),' Hz']);
    disp(['Drift = ',num2str(sum(wts_exp.*qwts_exp)*fr),', Diffusion = ',num2str(sum(wts_exp.*wts_exp.*qwts_exp)*fr)]);

    %Hyperbolic
    [wts_hyp qwts_hyp xopt_hyp]=match_distribution_1_var(drift,wts,'hyp',fr);
    disp(['Parameters for hyperbolic: weight parameter = ',num2str(xopt_hyp(1)),' V, firing rate = ',num2str(fr),' Hz']);
    disp(['Drift = ',num2str(sum(wts_hyp.*qwts_hyp)*fr),', Diffusion = ',num2str(sum(wts_hyp.*wts_hyp.*qwts_hyp)*fr)]);

    %Gaussian
    [wts_gauss qwts_gauss xopt_gauss]=match_distribution_1_var(drift,wts,'gauss',fr);
    disp(['Parameters for Gaussian: mu = ',num2str(xopt_gauss(1)),' V, sigma = ',num2str(abs(xopt_gauss(1))),' firing rate = ',num2str(fr),' Hz']);
    disp(['Drift = ',num2str(sum(wts_gauss.*qwts_gauss)*fr),', Diffusion = ',num2str(sum(wts_gauss.*wts_gauss.*qwts_gauss)*fr)]);

    %Lognormal
    [wts_logn qwts_logn xopt_logn]=match_distribution_1_var(drift,wts,'logn',fr);
    disp(['Parameters for Lognormal: mu = ',num2str(exp(xopt_logn(1))),' V, sigma = ',num2str(sqrt(exp(xopt_logn(1)))),' firing rate = ',num2str(fr),' Hz']);
    disp(['Drift = ',num2str(sum(wts_logn.*qwts_logn)*fr),', Diffusion = ',num2str(sum(wts_logn.*wts_logn.*qwts_logn)*fr)]);

else        %Match drift and diffusion, allowing firing rates to change
    diff = 0.0018;    %V^2/s, diffusion corresponding to drift-only matched epx distribution
   
    %Exponential
    [wts_exp qwts_exp xopt_exp]=match_distribution_1_var_fr(drift,diff,wts,'exp',fr);
    fr = xopt_exp(2);
    fr_exp=fr;
    disp(['Parameters for exponential: weight parameter = ',num2str(xopt_exp(1)),' V, firing rate = ',num2str(fr),' Hz']);
    disp(['Drift = ',num2str(sum(wts_exp.*qwts_exp)*fr),', Diffusion = ',num2str(sum(wts_exp.*wts_exp.*qwts_exp)*fr)]);

    %Hyperbolic
    [wts_hyp qwts_hyp xopt_hyp]=match_distribution_1_var_fr(drift,diff,wts,'hyp',fr);
    fr = xopt_hyp(2);
    fr_hyp=fr;
    disp(['Parameters for hyperbolic: weight parameter = ',num2str(xopt_hyp(1)),' V, firing rate = ',num2str(fr),' Hz']);
    disp(['Drift = ',num2str(sum(wts_hyp.*qwts_hyp)*fr),', Diffusion = ',num2str(sum(wts_hyp.*wts_hyp.*qwts_hyp)*fr)]);

    %Gaussian
    [wts_gauss qwts_gauss xopt_gauss]=match_distribution_1_var_fr(drift,diff,wts,'gauss',fr);
    fr=xopt_gauss(2);
    fr_gauss=fr;
    disp(['Parameters for Gaussian: mu = ',num2str(xopt_gauss(1)),' V, sigma = ',num2str(abs(xopt_gauss(1))),' firing rate = ',num2str(fr),' Hz']);
    disp(['Drift = ',num2str(sum(wts_gauss.*qwts_gauss)*fr),', Diffusion = ',num2str(sum(wts_gauss.*wts_gauss.*qwts_gauss)*fr)]);

    %Lognormal
    [wts_logn qwts_logn xopt_logn]=match_distribution_1_var_fr(drift,diff,wts,'logn',fr);
    fr=xopt_logn(2);
    fr_logn=fr;
    disp(['Parameters for Lognormal: mu = ',num2str(exp(xopt_logn(1))),' V, sigma = ',num2str(sqrt(exp(xopt_logn(1)))),' firing rate = ',num2str(fr),' Hz']);
    disp(['Drift = ',num2str(sum(wts_logn.*qwts_logn)*fr),', Diffusion = ',num2str(sum(wts_logn.*wts_logn.*qwts_logn)*fr)]);

end
