function [e,c,tau,k] = expFit(x,y,co,tauo,ko);
%[e,c,tau,k] = expFit(x,y,co,tauo,ko);
% Exponential fit to x,y data, can be compared to x,e
% with parameters such that e = c*exp(x/tau)+k

% set up simplex (3 parameters, so 4-vertex simplex)
widFac = 0.05;
if ko == 0
    ko = 0.0025; % dither k seed value (only one likely to be zero-valued)
end
x_o = [
    co            tauo                ko
    co*(1-widFac) tauo                ko
    co            tauo*(1-widFac)     ko
    co            tauo                ko*(1-widFac)
];

% define anonymous function call & handle
type = 2;
normFunction = @(x_o) expDiffs(x_o,y,x,type);

% perform fit via minimizing absolute value of sum of deviations
% and re-start minimization when done the first time all on L2 norm, then
% take a last pass on coefficient of determination
[x1,values,nEvals] = nmOpt(normFunction,x_o,1e-6);
[x2,values2,nEvals2] = nmOpt(normFunction,x1,1e-9);

% deconvolve output
c = x2(:,1);
tau = x2(:,2);
k = x2(:,3);
e = expCurve(mean(c),mean(tau),mean(k),x);

return

%%%%%%%%%%%%%%%PRIVATE FUNCTIONS%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = expCurve(c,tau,k,x);
y = c*exp(x/tau)+k;
return

%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
function val = expDiffs(x,y,t,type);
% calulates abs val diff or chi2 of fit vs. actual

% deconvolve x
c = x(1);
tau = x(2);
k = x(3);

if nargin == 3
    type = 2;
end

if type == 1
    val = sum(abs(y - expCurve(c,tau,k,t))); %L1 norm
elseif type == 8
    val = max(y - expCurve(c,tau,k,t));      %L8 norm
elseif type == 99 | type == 67
    val = 1-coeffDet(expCurve(c,tau,k,t),y); %coefficient of determination
else
    val = sum((y - expCurve(c,tau,k,t)).^2); %L2 norm ***
end

return
%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$


function [simplex,values,nEvals]=nmOpt(normFunction,simplex_o,vTol,mxEvals);
%[simplex,values,nEvals] = nmOpt(normFunction,simplex_o,valuesTol,maxEvals);
% Nelder-Mead simplex optimization scheme
%normFunction --> function handle to be minimized (use anonymous functions
%                 to tackle functions w/ parameters: e.g.:
%                                normFunction = @(x) foo(x,param1,param2);
%   simplex_o --> initial [d+1,d] simplex in search space
%   valuesTol --> exit tolerance on function values (default == 0.01)
%    maxEvals --> exit tolerance on function evaluations (default == 491)
%     simplex <== final [d+1,d] simplex, whose vertices give function 
%                 values all within valuesTol of one another (unless 
%                 maxEvals exceeded)
%      values <== [d+1,1] function values on simplex vertices
%      nEvals <== # of function evaluations taken

% (c) 23 Sep 2005, Curt A. L. Szuberla, 
% (c)              University of Alaska Fairbanks, all rights reserved

verNo = 1.0;
fnName = 'nmOpt';
%zcode(fnName);
%% version check
if nargin == 0
    disp(['(' fnName '.m) current version: ' num2str(verNo)]);
    return
end

% ascertain dimensions & initialize variables, set defaults
[diM,dim] = size(simplex_o);
simplex = simplex_o;
values = zeros(diM,1);
vTolDef = 0.01;   % function value exit tolerance
mxEvalsDef = 491; % max function evaluations before forced exit
EpS = eps*10^6;   % avoid divide-by-zero in exit tolerance checks
tempTol = mxEvalsDef;

% assign defaults if req'd
if nargin == 2
    vTol = vTolDef;
    mxEvals = mxEvalsDef;
elseif nargin == 3
    if isempty(vTol),vTol = vTolDef;end
    mxEvals = mxEvalsDef;
elseif nargin == 4
    if isempty(vTol),vTol = vTolDef;end
    if isempty(mxEvals),mxEvals = mxEvalsDef;end
end

% initialize values & sort
for k = 1:diM
    values(k) = normFunction(simplex(k,:));
end
nEvals = diM;
[values,idx] = sort(values);
simplex = simplex(idx,:);

% search pattern
while vTol < tempTol && nEvals < mxEvals
    % reflect worst point in simplex through opposite face
    [valueTest,simplex,values,nEvals] = ...
        metaReflect(-1,diM,dim,simplex,values,normFunction,nEvals);
    if valueTest >= values(dim) % refl pt worse than the next-worst vertex?
        valueTemp = values(diM); % set the worst vertex aside for a bit
        % contract in 1D along reflection line
        [valueTest,simplex,values,nEvals] = ...
            metaReflect(0.5,diM,dim,simplex,values,normFunction,nEvals);
        if valueTest >= valueTemp % contr. pt. worse than worst vertex?
            % contract entire simplex toward best vertex
            simplex(2:diM,:) = ... % no call to cols/rows here
                0.5*(simplex(2:diM,:) - ones(dim,1)*simplex(1,:));
            for k = 2:diM
                values(k) = normFunction(simplex(k,:));
            end
            nEvals = nEvals + dim;
        end
    elseif valueTest <= values(1) % was it better than the best vertex?
        % expansion in direction of reflection
        [valueTest,simplex,values,nEvals] = ...
            metaReflect(2,diM,dim,simplex,values,normFunction,nEvals);
    end
    % re-sort on exit
    [values,idx] = sort(values);
    simplex = simplex(idx,:);
    % calculate exit condition
    tempTol = abs(values(diM)-values(1))/(mean(abs(values))+EpS);
end

return

%% private functions
function [valueTest,simplex,values,nEvals] = ...
    metaReflect(scaleFactor,diM,dim,simplex,values,normFunction,nEvals);

% calculate meta-reflection point
scaleA = (1-scaleFactor)/dim;
simplexTest = sum(simplex)*scaleA - ...
    simplex(diM,:)*(scaleA-scaleFactor);

% attempt metaReflection
valueTest = normFunction(simplexTest);
nEvals = nEvals + 1;

% check result for improvement
if values(diM) > valueTest
    simplex(diM,:) = simplexTest;
    values(diM) = valueTest;
end

return