function [uSSlist, ySSlist, xSSlist, vallist]=doseResponse(m, Umin, Umax, xguess, userOptions)
%DOSERESPONSE computes the steady state response to the given input
%[uSSlist, ySSlist, xSSlist, vallist]=doseResponse(m, Umin, Umax,  xguess, ...)
%    computes the steady state response to the given input.
%
% m      - The kronecker bio model
% Umin   - The minimum (non zero) input concentration
% Umax   - The maximum input concnetration
%
% xguess - A guess for the steady state at Umin.  If the system is
%          multistable this parameter will chage which state is found
%
% 
% uSSlist - The concentration the steady states were evaluated at.
% ySSlist - The steady state values of the model outputs.
% xSSlist - The steady state values of all the model concentrations.
% vallist - The value of the residual norm(m.f(0,xss,uss))
%
% [...]=doseResponse(..., options) pass extra options
%   to the optimizer.  see optimset
%
% See also: optimset

% This work is licensed under the Creative Commons Attribution-Noncommercial-No Derivative 
% Works 3.0 United States License. To view a copy of this license, visit 
% http://creativecommons.org/licenses/by-nc-nd/3.0/us/ or send a letter to Creative Commons,
% 171 Second Street, Suite 300, San Francisco, California, 94105, USA.


% This code assumes this data is in columns
m.ic   = m.ic(:);
xguess = xguess(:);

% Set default optimization options
defOptions =  optimset('Hessian',    'on', ...
                       'GradObj',   'on',...
                       'DerivativeCheck', 'off',...
                       'Display',    'off', ...
                       'TolFun',     1e-20, ...
                       'TolX',       1e-20, ...
                       'MaxFunEvals', 8000, ...
                       'MaxIter',     2000, ...
                       'LargeScale', 'on', ...
                       'PrecondBandWidth', Inf);

% Use the user optimization options if provided
if(nargin > 4)
    options = optimset(defOptions, userOptions);
else
    options = defOptions;
end

% Find the mass conserved species and create constraints to preserve
% these quantites
N   = null(m.S', 'r')';
Aeq = N;
beq = full(N*m.ic);
   
% Solve for each value of the input
uSS     = Umin;     %Start at the begining
maxStep = 1;        %Maximum step is the size of the current concentration
uStep   = maxStep;  %Start with the largest step
i       = 0;       
u       = Umin;

%Q = 1e20*diag((typicalX./min(typicalX)).^(-2));
%Q=eye(m.nX);
xss = xguess;

fprintf(1, 'Step\tInput        \tError\n');
while u <= Umax
    i = i+1;

    while true
        xguess   = guessXss(xss, u, uSS);
        [v, val] = fminunc(@func4, xguess,options); 
        
        if(val>1e-6)
            if(i == 1)
                error('Failed to converge on first step, try another xGuess');
            end
            uStep  = uStep/2;
            uSS   = u*(1+uStep);
        else
            u            = uSS;
            xss          = v(1:m.nX);
            ySSlist(i,:) = m.c*xss;
            xSSlist(i,:) = xss;
            vallist(i)   = val;
            uSSlist(i)   = u;
            
            fprintf(1, '%04d\t%0.8e\t%0.8e\n', i, uSS, val);

            if(u >= Umax)
                % Then we are done!
                return;
            else
                % Step and try a bigger step this time
                uStep = min(uStep*1.8,maxStep);
                uSS   = min(u*(1+uStep), Umax);
                break;
            end
        end
    end
end

return
    
    function xguess =guessXss(x0, u0, u1)
        % Use a linerization of the model to guess the next steady state

        f    = m.f(0,x0,u0);
        dfdx = full(m.dfdx(0,x0,u0));
        dfdu = m.dfdu(0,x0,u0);
        du   = u1-u0;
        dX   = solveApprox(full(dfdx), -(f+dfdu*du), 1e-12);
        xguess = x0+dX;
    end

    
    function [F, G, H] = func4(x)
        % This objective function seeks to minimize the flux
        % and the mass constraint violation
        
        f=m.f(0, x, uSS);
        dfdx = m.dfdx(0,x,uSS);
        ftd2fdx2 = reshape(f'*2*m.a2,m.nX,m.nX);
        
        err = Aeq*x-beq;
        
        F = (f'*f + err'*err);
        G = 2*f'*dfdx + 2*err'*Aeq;
        H = 2*(ftd2fdx2+dfdx'*dfdx+Aeq'*Aeq);
    end

%     % The function being optimized over.
%     function [F, J] = func(v)
%         x     = v(1:m.nX);
%         lamda = v((m.nX+1):end);
%         
%         f=m.f(0, x, uSS);
%         dfdx = m.dfdx(0,x,uSS);
%         d2fdx2 = full(m.d2fdx2(0,x,uSS));
%         ftd2fdx2 = reshape(f'*d2fdx2',m.nX,m.nX);
%         F = [2*f'*dfdx+lamda'*M, (M*x-beq)']';
%         J = [2*ftd2fdx2+2*dfdx'*dfdx,              M';
%              M,                        zeros(numCons)];
%         J = full(J);
%     end

%     % The function being optimized over.
%     function [F, J] = func2(x)
%         F=full(m.f(0, x, uSS));
%         J=full(m.dfdx(0,x,uSS));
%     end

%     function [F, G, H] = func3(x)
%         f=m.f(0, x, uSS);
%         dfdx = m.dfdx(0,x,uSS);
%         ftd2fdx2 = reshape(f'*Q*m.a2,m.nX,m.nX);
%         
%         F = (f'*Q*f);
%         G = 2*f'*Q*dfdx;
%         H = 2*(ftd2fdx2+dfdx'*Q*dfdx);        
%     end
end