function [a, rSquared] = fitDifferenceOfGaussianConstraint(x, y, gaussianSign)
%%
% DoG = @(a) a(1) * exp(-((x - a(2)) / a(3)).^2) + ...
%            a(4) * exp(-((x - a(2)) / a(5)).^2);
% lsq = @(a) sum(F(a) - y).^2;

DoG = @(a) differenceOfGaussians1D(a, x);
% DoG = @(a) a(1) * exp(-((x - a(2)) / a(3)).^2) + ...
%            a(4) * exp(-((x - a(2)) / a(5)).^2);
lsq = @(a) sum((DoG(a) - y).^2);

[minY, argMinY] =  min(y);
[maxY, argMaxY] = max(y);
xStep = mean(diff(x));
positiveAmpBounds = [0, 2 * range(y)];
negativeAmpBounds = [- 2 * range(y), 0];
positionBounds = x([1, end]);

if gaussianSign > 0
    centerAmpBounds = positiveAmpBounds;
    centerWidthBounds = [xStep, range(x)];
    surroundAmpBounds = negativeAmpBounds;
    surroundWidthBounds = [2 * xStep, range(x)];
    a0 = [maxY, argMaxY, xStep, minY, range(x) / 2];
else
    centerAmpBounds = negativeAmpBounds;
    centerWidthBounds = [2 * xStep, range(x)];
    surroundAmpBounds = positiveAmpBounds;
    surroundWidthBounds = [2 * xStep, range(x)];
    a0 = [minY, argMinY, xStep, maxY, range(x) / 2];
end

% center should be narrower than surround: a(3) < a(5) <=> a(3) - a(5) < 0
A = [0 0 3 0 -5]; 
b = 0;
Aeq = []; 
beq = [];
lb = [centerAmpBounds(1), positionBounds(1), centerWidthBounds(1), ...
      surroundAmpBounds(1), surroundWidthBounds(1)];
ub = [centerAmpBounds(2), positionBounds(2), centerWidthBounds(2), ...
      surroundAmpBounds(2), surroundWidthBounds(2)];

[a, fitResidue] = fmincon(lsq, a0, A, b, Aeq, beq, lb, ub);

rSquared = 1 - fitResidue / sum((y - mean(y)).^2);

end