function [ySSlist, xSSlist, vallist]=doseResponse(m, uVals, xguess, userOptions)

%DOSERESPONSE computes the steady state response to the given input
%
% [ySSlist, xSSlist, vallist]=doseResponse(m, uVals, xguess, ...) computes
%    the steady state response to the given input using newtons method.
%
% m      - The kronecker bio model
% uVals  - A list of input concentrations (one for each dose)
% xguess - A guess for the steady state.  If the system is multistable this
%          parameter will chage which state is found
%
% 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('Jacobian',    'on', ...
                           'Display',    'off', ...
                           'TolFun',     1e-5, ...
                           'TolX',       1e-5, ...
                           'LargeScale', 'off', ...
                           'MaxFunEvals', 1000, ...
                           'MaxIter',     1000, ...
                           'DerivativeCheck', 'on', ...
                           'NonlEqnAlgorithm', 'gn');

    if(nargin > 4)
        options = optimset(defOptions, userOptions);
    else
        options = defOptions;
    end


    % Compute the reachable subspace
    Aeq = null(full(m.S'))';
    N   = null(Aeq); 
    x0  = m.ic;

    ySSlist = zeros(length(uVals), m.nY);
    xSSlist = zeros(length(uVals), m.nX);
    vallist = zeros(size(uVals));
    
    % Solve for each value of the input
    for i = 1:length(uVals)
        i
        uSS = uVals(i);
        dV  = N'*(xguess - x0);
        dX  = N*dV;
        F   = m.f(0, x0 + N*dV, uSS);
%         xss = x0 + N*dV;
        xss = x0;
        
        while norm(F) > 1e-20 && norm(dX) > 1e-20
           
           F   = m.f   (0, xss, uSS);
           J   = m.dfdx(0, xss, uSS);

           dX  = solveApprox(J, -F, 1e-9);
           xss = x0 + dX; 
           disp([norm(F), norm(dX)]);
        end
        
           disp('**************************');
        ySSlist(i,:) = m.c*xss;
        xSSlist(i,:) = xss;
        vallist(i)   = norm(F);

        if(i < length(uVals))
            dfdx   = m.dfdx(0,xss,uVals(i));
            dfdu   = m.dfdu(0,xss,uVals(i));
            du     = uVals(i+1)-uVals(i);
            v      = solveApprox(N'*dfdx*N, -N'*dfdu*du, 1e-9);
            xguess = x0 + N*v
        else
            xguess = x0 + N*v;
        end
    end
end

% function v = solveApprox(A, b, tol)
% 
%    [U,S,V] = svds(A);
% 
%    s  = spdiags(S);
% 
%    sn = s(s./s(1)>tol);
%    n  = length(sn);
% 
%    Sn = spdiags(sn,0,n,n);
%    Un = U(:, 1:n);
%    Vn = V(:, 1:n);
% 
%    v = Vn*Sn^(-1)*Un'*b;
% end