% softplus_function: Fits a softplus function to the input data and
% extracts the softplus growth rate.

% This function fits a softplus function to the given time series data and
% calculates the softplus growth rate.

function [softplus_growth] = softplus_function(years,values)


%start_asymp = max(values);
%start_inflect = min(years) + (max(years) - min(years));
%start_growth = 1/max(values) - min(values);
max_asymp = max(values);
% Define softplus equation
%softplus_equation =  fittype('(L ./ k) * log(1 + exp(k * (x - thalf)))', ...
 %   'dependent', 'y', 'independent', 'x', ...
  %  'coefficients', {'L', 'k', 'thalf'});

  ft = fittype('(L / k)*log(1 + exp(k * (x - thalf)))', 'coeff', ...
      {'L', 'k', 'thalf'});
  options = fitoptions(ft);
  options.Lower = [max_asymp 0 1900];
  options.Upper = [max_asymp*2 2 2100 ];

    try 
    % Fit the softplus equation to the input data
      % softpl =  fit(years, values, ...
       %         '(L / k)(log(1 + exp(k * (x - thalf))))', ...
        %        "StartPoint",[start_asymp start_growth start_inflect]) ;
        %softpl =  fit(years, values, ...
         %       '(L / k) (log[1 + exp(k * {x - thalf})]');
        
        softpl=fit(years, values, ft, options);
        % Extract coefficients from softpl fit
        softplus_coefficients = coeffvalues(softpl);
        % Extract growth rate
        softplus_growth = softplus_coefficients(2);
        
    catch

        % If softplus fitting fails, define softplus fit as zero
        softplus_growth = 0;
    end
end