function [best_params] = boltzmannFit(contrasts, onset_fi)
% use a gradient descent approach to fit the fi curve owith a Boltzmann
% sigmoidal function
%
% Arguments:
% constrasts: the x-values (must be sorted ascendingly)
% onset_fi: the measured neuronal responses
%
% Returns the parameters leading to the best fit.

% we will make some educated guesses for the start parameters; minimum:
% minimum fi response; saturation: maximum fi response; inflection: zero
% contrast; slope: response range/contrast range
start_params = [min(onset_fi), max(onset_fi), 0, ...
    (max(onset_fi) - min(onset_fi))/(max(contrasts) - min(contrasts))];

best_params = gradientDescent(start_params, contrasts, onset_fi);


%% subfunctions that could be separate m-files, but are only used in this context

    function bp = gradientDescent(p, x, y)
    % performs the gradient descent. Input arguments: 
    % p: the list of starting parameters for the search
    % x: the x-values
    % y: the measured y-values
    %
    % returns the parameters leading to the best fit
        
       delta = 0.001;
       gradient = getGradient(p, x, y, delta);
       eps = [0.001 0.001, 0.00001, 0.00001]; % key point; use smaller changes for the more delicate parameters
       while norm(gradient) > 0.5
           p = p - eps .* gradient;
           gradient = getGradient(p, x, y, delta);
       end  
       bp = p;
    end

    function g = getGradient(p, x, y, delta)
    % calculate the gradient in the error surface for a small step
    % along all parameters
    % p: the parameter list
    % x: the x-values
    % y: the measured y-values
    % delta: the step used or the gradient
    %
    % returns the gradient
        g = zeros(size(p));
        
        for i = 1:length(p)
            np = p;
            np(i) = np(i) + delta;
            g(i) = (costFunction(np, x, y) - costFunction(p, x, y)) / delta;
        end
    end

    function cost = costFunction(p, x, y)
        % cost function that relates the model to the data. Applies a
        % Boltzman model and calculates the error in a mean square error
        % sense.
        % p: the set of parameters used for the model
        % x: the x-values
        % y: the measured y-values
        %
        % returns the costs related to that parameterization
        
        y_est = boltzmannModel(p, x);
        cost = msqError(y, y_est);
    end

    function e = msqError(y, y_est)
    % calculates the deviation between y and y_est in the mean square
    % error sense
    % y and y_est must have the same size
    %
    % returns the error
        e = mean((y - y_est).^2);
    end

    

end