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('LargeScale', 'on', ...
                       'GradObj',    'on', ...
                       'Hessian',    'on', ...
                       'Display',    'on', ...
                       'TolFun',     1e-15, ...
                       'TolX',       1e-15, ...
                       'MaxFunEvals', 1e5, ...
                       'MaxIter',     1e4, ...
                       'DerivativeCheck', 'off');

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

% Compute the reachable subspace
M   = massConservation(m)';
b   = M*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);
    objFcn = constructObj(m, uSS);
    [xss G] = fmincon(objFcn, xguess, [], [], M, b, [], [], [], options);
    ySSlist(i,:) = m.c*xss;
    xSSlist(i,:) = xss;
    vallist(i)   = sqrt(G);
    
    semilogx(uVals(1:i), ySSlist(1:i,:), '-')
    set(gca, 'XLim', [uVals(1) uVals(end)])
    drawnow
    
    % either update at the most recent solution
    xguess       = xss;
    
end
end

function fcn = constructObj(m, u)
fcn = @obj;

    function [f J H] = obj(x)
        f = m.f(0, x, u);
        if nargout > 1
            dfdx = m.dfdx(0,x,u);
            J = 2*f'*dfdx;
            if nargout > 2
                d2fdx2 = 2*m.a2;
                H = 2*dfdx'*dfdx + reshape(2*f'*d2fdx2, m.nX, m.nX);
            end
        end
        f = f'*f;
    end
end

% function v = solveApprox(A, b, tol, nX)
% 
%    [U,S,V] = svds(A, nX);
% 
%    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