function fitOutputStruct = fitGaussian1(x, yArray, gaussianSign, varargin)
% Fit one dimensional Gaussian to trial data or to many cells (ROIs).
switch nargin
    case {1, 2, 3}
        dataStructure = 'cells';
    case 4
        dataStructure = varargin{1};
end

switch dataStructure
    case 'cells'
    case 'trials'
        originalArray = yArray;
        yArray = mean(yArray, 1);
end

if ~exist('gaussianSign', 'var') || isempty(gaussianSign) || gaussianSign == 0
    convexityProxy = approximateConvexity(yArray);
else
    convexityProxy = gaussianSign * ones(1, size(yArray, 1));
end

fitOutputStruct(size(yArray, 1)).fit = nan;
fitOutputStruct(size(yArray, 1)).gof = nan;
fitOutputStruct(size(yArray, 1)).info = nan;


for iRF = 1: size(yArray, 1)
%     x = (1 : numel(yArray(iRF, :)));
    % Idea: recover the smoothing filter through deconvolution to
    % deconvolve the fitted gaussian.
%     y = smoothdata(yArray(iRF, :), 'gaussian');
    y = (yArray(iRF, :));
    options = fitoptions('gauss1');
    options.Lower = [-10 1 1];
    options.Upper = [10 numel(y) - 1 numel(y)];
    if sign(convexityProxy) == 1
        options.Lower(1) = 0;
        [~, maxInd] = max(y);
        options.StartPoint = [max(y) maxInd 1];
    else
        options.Upper(1) = 0;
        [~, minInd] = min(y);
        options.StartPoint = [min(y) minInd 1];
    end
    
    switch dataStructure
        case 'cells'
            [curve, gof, output] = fit(x', y','gauss1', options);
        case 'trials'
            x = repelem(x, size(originalArray, 1));
            [curve, gof, output] = fit(x', ColumnVector(originalArray),'gauss1', options);
    end

    fitOutputStruct(iRF).fit = curve;
    fitOutputStruct(iRF).gof = gof;
    fitOutputStruct(iRF).info = output;
end

end



function convexityProxy = approximateConvexity(yArray)
% Approximate the sign of the curvature, by summing lowest to highest
% values.
firstDecile = @(x) quantile(x, 0.1, 2);
lastDecile = @(x) quantile(x, 0.9, 2);

% concavity = gradient(gradient(yArray));
lowValInds = yArray <= firstDecile(yArray);
highValInds = yArray >= lastDecile(yArray);

lowSum = zeros(1, size(yArray, 1));
highSum = zeros(1, size(yArray, 1));

% Carefull with smoothing level since bars and noise have different spatial
% sampling rates.
for iRF = 1: size(yArray, 1)
%     x = (1 : numel(yArray(iRF, :)));
%     y = smoothdata(yArray(iRF, :), 'gaussian', 3);
    y = (yArray(iRF, :));
    lowSum(iRF) = sum(y(lowValInds(iRF, :)));
    highSum(iRF) = sum(y(highValInds(iRF, :)));
end
convexityProxy = lowSum + highSum;
end