outdir = './RFModelFigs/';
outputIgor = false;

nlParams = [70.7 3.6667 -1.90 -1.6];

x = linspace(-1,1,50); 
softscale = 3;

nlfun = @(a,x)(a(1)*normcdf(a(2)*x+a(3),0,1)+a(4));
nldiff = @(a,x)(a(1)*normpdf(a(2)*x+a(3),0,1));

% Relu nonlinearity
relu = @(x)(max(0,x));
% Softplus nonlinearity
softplus = @(x)(log(1 + exp(x*softscale)));

y = nlfun(nlParams,x);
ydiff = [0, diff(y)]; % First derivative (Gaussian)
ydiff2 = [0,0, diff(diff(y))]; % Second derivative

% Linear
yL = softscale*x;
ydiffL = diff(yL(1:2))*ones(size(yL)); % First derivative
ydiffL2 = zeros(size(yL)); % Second derivative

% Relu
yR = relu(x);
ydiffR = [0,diff(yR)]; % First derivative
ydiffR2 = [0,0,diff(diff(yR))]; % Second derivative

% Softplus
yS = softplus(x);
ydiffS = 1./(1 + exp(-x*softscale)); % First derivative is logistic
ydiffS2 = (1 - 1./(1 + exp(-x*softscale))).*(1./(1 + exp(-x*softscale)));

couplingGain = 0.4;
xval = 0.5;
nSubunits = 50;
minVals = 1:nSubunits/2;

r = zeros(size(minVals));
rL = zeros(size(minVals));
rR = zeros(size(minVals));
rS = zeros(size(minVals));
for k = 1 : length(minVals)
    nMin = minVals(k);
    nMax = nSubunits - nMin;
    

    % Create converting and diverging coupled networks.
    xConverge = zeros(nSubunits,1);
    xDiverge = zeros(nSubunits,1);


    xConverge(1:nMin) = xval + couplingGain*xval*nMax/nMin;
    xConverge(nMin+1:nSubunits) = xval - couplingGain*xval;
    xDiverge(1:nMin) = xval - couplingGain*xval;
    xDiverge(nMin+1:nSubunits) = xval + couplingGain*xval*nMin/nMax;

    % Confirm that they sum to the same value.
%     disp(['convergent sum: ', num2str(sum(xConverge)), ' == divergent sum: ', num2str(sum(xDiverge))]);

    % Pipe through measured nonlinearity.
    yConverge = nlfun(nlParams, xConverge);
    yDiverge = nlfun(nlParams, xDiverge);
    
    % Pipe through Relu nonlinearity.
    yCRelu = relu(xConverge);
    yDRelu = relu(xDiverge);
    
    % Pipe through softplus nonlinearity.
    yCSoft = softplus(xConverge);
    yDSoft = softplus(xDiverge);

    rL(k) = sum(xDiverge) / sum(xConverge); % Linear
    r(k) = sum(yDiverge) / sum(yConverge); % Nonlinear
    rR(k) = sum(yDRelu) / sum(yCRelu); % Relu
    rS(k) = sum(yDSoft) / sum(yCSoft); % Softplus
end

%%
figure(100); clf;
xidx = (x >= 0.05 & x <= 0.55);

inflection = 0.2;
stepSize = 0.1;
x1 = find(abs(x - (inflection-stepSize)) == min(abs(x - (inflection-stepSize))), 1);
x2 = find(abs(x - (inflection)) == min(abs(x - (inflection))), 1);
x3 = find(abs(x - (inflection+stepSize)) == min(abs(x - (inflection+stepSize))), 1);

xTop = x([1 1 x3 x3 x2 x2 1]);
xBottom = x([1 1 x2 x2 x1 x1 1]);

subplot(221)
hold on
addLineToAxis(x,y, 'y',gca,'k','-','none');
addLineToAxis(x(xidx),y(xidx), 'highlight',gca,[0.8 0 0],'-','none');
addLineToAxis(xTop,y([x2 x3 x3 1 1 x2 x2]), 'top', gca, 'k','-','none');
addLineToAxis(xBottom, y([x1 x2 x2 1 1 x1 x1]), 'bottom', gca, 'k','-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[min(y) max(y)],'yzero',gca,'k','--','none');
hold off;
xlabel('input (x)'); ylabel('output [N(x)]');
if outputIgor
    makeAxisStruct(gca, 'VoltageVsRelease', 'basedir', outdir);
end

subplot(223)
hold on
addLineToAxis(x,ydiff, 'y',gca,'k','-','none');
addLineToAxis(x(xidx),ydiff(xidx), 'highlight',gca,[0.8 0 0],'-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[0 max(ydiff)],'yzero',gca,'k','--','none');
hold off;
xlabel('input (x)'); ylabel('dN(x)/dx');
if outputIgor
    makeAxisStruct(gca, 'VoltageVsDeltaRelease', 'basedir', outdir);
end

subplot(224)
hold on
addLineToAxis(x,ydiff2, 'y',gca,'k','-','none');
addLineToAxis(x(xidx),ydiff2(xidx), 'highlight',gca,[0.8 0 0],'-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[0 max(ydiff2)],'yzero',gca,'k','--','none');
hold off;
xlabel('input (x)'); ylabel('d2N(x)/dx2');
if outputIgor
    makeAxisStruct(gca, 'VoltageVsDeltaRelease2', 'basedir', outdir);
end

subplot(222)
hold on
% addLineToAxis(nSubunits./(minVals*2),r,'y',gca,'k','-','none');
addLineToAxis(minVals/nSubunits*100,rL,'linear',gca,[132 159 173]/255,'-','none');
addLineToAxis(minVals/nSubunits*100,r,'nonlinear',gca,[39 100 123]/255,'-','none');
addLineToAxis(minVals/nSubunits*100,rR,'relu',gca,[202 53 66]/255,'-','none');
addLineToAxis(minVals/nSubunits*100,rS,'soft',gca,[87 87 95]/255,'-','none');
hold off;
xlabel('% proximal subunits'); ylabel('Rdiverge : Rconverge');
if outputIgor
    makeAxisStruct(gca, 'SubunitDivergence', 'basedir', outdir);
end

% Linear input output
figure(101); clf;

subplot(221)
hold on
addLineToAxis(x,yL, 'y',gca,'k','-','none');
addLineToAxis(x(xidx),yL(xidx), 'highlight',gca,[0.8 0 0],'-','none');
addLineToAxis(xTop,yL([x2 x3 x3 1 1 x2 x2]), 'top', gca, 'k','-','none');
addLineToAxis(xBottom, yL([x1 x2 x2 1 1 x1 x1]), 'bottom', gca, 'k','-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[min(yL) max(yL)],'yzero',gca,'k','--','none');
hold off;
xlabel('input (x)'); ylabel('output [N(x)]');
if outputIgor
    makeAxisStruct(gca, 'VoltageVsReleaseL', 'basedir', outdir);
end

subplot(222)
hold on
addLineToAxis(x,ydiffL, 'y',gca,'k','-','none');
addLineToAxis(x(xidx),ydiffL(xidx), 'highlight',gca,[0.8 0 0],'-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[0 max(ydiffL)*2],'yzero',gca,'k','--','none');
hold off;
xlabel('input (x)'); ylabel('dN(x)/dx');
if outputIgor
    makeAxisStruct(gca, 'VoltageVsDeltaReleaseL', 'basedir', outdir);
end

subplot(223)
hold on
addLineToAxis(x,ydiffL2, 'y',gca,'k','-','none');
addLineToAxis(x(xidx),ydiffL2(xidx), 'highlight',gca,[0.8 0 0],'-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[0 max(ydiffL2)*2],'yzero',gca,'k','--','none');
hold off;
xlabel('input (x)'); ylabel('d2N(x)/dx2');
if outputIgor
    makeAxisStruct(gca, 'VoltageVsDeltaReleaseL2', 'basedir', outdir);
end

figure(102); clf;
% Relu
subplot(221)
hold on
addLineToAxis(x,yR, 'y',gca,'k','-','none');
addLineToAxis(x(xidx),yR(xidx), 'highlight',gca,[0.8 0 0],'-','none');
addLineToAxis(xTop,yR([x2 x3 x3 1 1 x2 x2]), 'top', gca, 'k','-','none');
addLineToAxis(xBottom, yR([x1 x2 x2 1 1 x1 x1]), 'bottom', gca, 'k','-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[min(yR) max(yR)],'yzero',gca,'k','--','none');
hold off;
xlabel('input (x)'); ylabel('output [N(x)]');
if outputIgor
    makeAxisStruct(gca, 'VoltageVsReleaseR', 'basedir', outdir);
end

subplot(222)
hold on
addLineToAxis(x,ydiffR, 'y',gca,'k','-','none');
addLineToAxis(x(xidx),ydiffR(xidx), 'highlight',gca,[0.8 0 0],'-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[0 max(ydiffR)*2],'yzero',gca,'k','--','none');
hold off;
xlabel('input (x)'); ylabel('dN(x)/dx');
if outputIgor
    makeAxisStruct(gca, 'VoltageVsDeltaReleaseR', 'basedir', outdir);
end

subplot(223)
hold on
addLineToAxis(x,ydiffR2, 'y',gca,'k','-','none');
addLineToAxis(x(xidx),ydiffR2(xidx), 'highlight',gca,[0.8 0 0],'-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[0 max(ydiffR2)*2],'yzero',gca,'k','--','none');
hold off;
xlabel('input (x)'); ylabel('d2N(x)/dx2');
if outputIgor
    makeAxisStruct(gca, 'VoltageVsDeltaReleaseR2', 'basedir', outdir);
end

figure(103); clf;

subplot(221)
hold on
addLineToAxis(x,yS, 'y',gca,'k','-','none');
addLineToAxis(x(xidx),yS(xidx), 'highlight',gca,[0.8 0 0],'-','none');
addLineToAxis(xTop,yS([x2 x3 x3 1 1 x2 x2]), 'top', gca, 'k','-','none');
addLineToAxis(xBottom, yS([x1 x2 x2 1 1 x1 x1]), 'bottom', gca, 'k','-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[min(0,min(yS)) max(yS)],'yzero',gca,'k','--','none');
hold off;
xlabel('input (x)'); ylabel('output [N(x)]');
if outputIgor
    makeAxisStruct(gca, 'VoltageVsReleaseS', 'basedir', outdir);
end

subplot(222)
hold on
addLineToAxis(x,ydiffS, 'y',gca,'k','-','none');
addLineToAxis(x(xidx),ydiffS(xidx), 'highlight',gca,[0.8 0 0],'-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[0 max(ydiffS)],'yzero',gca,'k','--','none');
hold off;
xlabel('input (x)'); ylabel('dN(x)/dx');
if outputIgor
    makeAxisStruct(gca, 'VoltageVsDeltaReleaseS', 'basedir', outdir);
end

subplot(223)
hold on
addLineToAxis(x,ydiffS2, 'y',gca,'k','-','none');
addLineToAxis(x(xidx),ydiffS2(xidx), 'highlight',gca,[0.8 0 0],'-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[0 max(ydiffS2)*2],'yzero',gca,'k','--','none');
hold off;
xlabel('input (x)'); ylabel('d2N(x)/dx2');
if outputIgor
    makeAxisStruct(gca, 'VoltageVsDeltaReleaseS2', 'basedir', outdir);
end

subplot(224)
hold on
% addLineToAxis(nSubunits./(minVals*2),r,'y',gca,'k','-','none');
addLineToAxis(minVals/nSubunits*100,rL,'linear',gca,[132 159 173]/255,'-','none');
addLineToAxis(minVals/nSubunits*100,r,'nonlinear',gca,[39 100 123]/255,'-','none');
addLineToAxis(minVals/nSubunits*100,rR,'relu',gca,[202 53 66]/255,'-','none');
addLineToAxis(minVals/nSubunits*100,rS,'soft',gca,[87 87 95]/255,'-','none');
hold off;
xlabel('% proximal subunits'); ylabel('Rdiverge : Rconverge');
if outputIgor
    makeAxisStruct(gca, 'SubunitDivergence2', 'basedir', outdir);
end

% Illustrate measured nonlinearity for converging/diverging conditions.
figure(105); clf;
% xidx = (x >= 0.05 & x <= 0.55);

inflection = 0.5;
stepUp = 0.5;
stepDown = 0.2;
x1 = find(abs(x - (inflection-stepDown)) == min(abs(x - (inflection-stepDown))), 1);
x2 = find(abs(x - (inflection)) == min(abs(x - (inflection))), 1);
x3 = find(abs(x - (inflection+stepUp)) == min(abs(x - (inflection+stepUp))), 1);
xTop = x([1 1 x3 x3 x2 x2 1]);
xBottom = x([1 1 x2 x2 x1 x1 1]);

subplot(221)
hold on
addLineToAxis(x,y, 'y',gca,'k','-','none');
addLineToAxis(x(xidx),y(xidx), 'highlight',gca,[0.8 0 0],'-','none');
addLineToAxis(xTop,y([x2 x3 x3 1 1 x2 x2]), 'top', gca, 'k','-','none');
addLineToAxis(xBottom, y([x1 x2 x2 1 1 x1 x1]), 'bottom', gca, 'k','-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[min(y) max(y)],'yzero',gca,'k','--','none');
hold off;
xlabel('input (x)'); ylabel('output [N(x)]');
if outputIgor
    makeAxisStruct(gca, 'VoltageConverge', 'basedir', outdir);
end

inflection = 0.5;
stepUp = 0.1;
stepDown = 0.2;
x1 = find(abs(x - (inflection-stepDown)) == min(abs(x - (inflection-stepDown))), 1);
x2 = find(abs(x - (inflection)) == min(abs(x - (inflection))), 1);
x3 = find(abs(x - (inflection+stepUp)) == min(abs(x - (inflection+stepUp))), 1);

xTop = x([1 1 x3 x3 x2 x2 1]);
xBottom = x([1 1 x2 x2 x1 x1 1]);

subplot(222)
hold on
addLineToAxis(x,y, 'y',gca,'k','-','none');
addLineToAxis(x(xidx),y(xidx), 'highlight',gca,[0.8 0 0],'-','none');
addLineToAxis(xTop,y([x2 x3 x3 1 1 x2 x2]), 'top', gca, 'k','-','none');
addLineToAxis(xBottom, y([x1 x2 x2 1 1 x1 x1]), 'bottom', gca, 'k','-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[min(y) max(y)],'yzero',gca,'k','--','none');
hold off;
xlabel('input (x)'); ylabel('output [N(x)]');
if outputIgor
    makeAxisStruct(gca, 'VoltageDiverge', 'basedir', outdir);
end

% Make a pdf of the input-output functions for illustration purposes.
figure(106); clf;

subplot(221)
hold on
addLineToAxis(x,y, 'y',gca,'k','-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[min(y) max(y)],'yzero',gca,'k','--','none');
hold off; axis tight; axis square;
xlabel('input (x)'); ylabel('output [N(x)]');

subplot(222)
hold on
addLineToAxis(x,yL, 'y',gca,'k','-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[min(yL) max(yL)],'yzero',gca,'k','--','none');
hold off; axis tight; axis square;
xlabel('input (x)'); ylabel('output [N(x)]');

subplot(223)
hold on
addLineToAxis(x,yR, 'y',gca,'k','-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[min(yR) max(yR)],'yzero',gca,'k','--','none');
hold off; axis tight; axis square;
xlabel('input (x)'); ylabel('output [N(x)]');

subplot(224)
hold on
addLineToAxis(x,yS, 'y',gca,'k','-','none');
addLineToAxis(x([1 end]),zeros(1,2),'xzero',gca,'k','--','none');
addLineToAxis(zeros(1,2),[min(yS) max(yS)],'yzero',gca,'k','--','none');
hold off; axis tight; axis square;
xlabel('input (x)'); ylabel('output [N(x)]');

print('InputOutputs','-dpdf')