function result = ligap(time_vector,data_vector,labels)
%LIGAP  The LIGAP method is widely applicable to quantify differential
%       time-course dynamics of many types of datasets and generalizes
%       to any number of conditions
%   
%   usage: result = ligap(time_vector, data_vector, labels)
%
%   where:
%       result          is the result vector, probabilities has the
%                       posterior probabilities for the models specified
%                       by models variable (see below), profiles has the
%                       means and standard deviations of the fitted GPs
%                       per condition (evaluation time points are given)
%       time_vector     is a cell vector of length N (N conditions)
%                       containing the measurement time points
%       data_vector     is a cell vector of length N (N conditions) 
%                       containing the corresponding measurements
%                       (log2 transformed for microarrays) for the
%                       time points given in time_vector, order should be
%                       the same
%       labels          is a cell vector of length N holding the labels of
%                       the experimental conditions
%
%   description of the model notation:
%       {cond1cond2cond3}: all conditions are similar
%       {cond1cond2}{cond3}: cond3 is different cond1 and cond2 (which are
%                             similar)
%       {cond1cond3}{cond2}: cond2 is different cond1 and cond3 (which are
%                             similar)
%       {cond1}{cond2cond3}: cond1 is different cond2 and cond3 (which are
%                             similar)
%       {cond1}{cond2}{cond3}: all conditions are different
%
%   examples:
%       two conditions, three time points (1,2,3) and three replicates:
%           result = ligap({[1,1,1,2,2,2,3,3,3]',[1,1,1,2,2,2,3,3,3]'},
%                           {rand(9,1),rand(9,1)},{'cond1','cond2'})
%
%       three conditions, three time points (1,2,3), and varying number of
%       replicates
%           result = ligap({[1,1,1,2,2,2,3,3,3]',[1,1,2,2,3,3]',
%                           [1,1,1,2,3,3,3]'},{rand(9,1),rand(6,1),
%                            rand(7,1)},{'cond1','cond2','cond3'})
%
%   The GPML Matlab code is written by Carl Edward Rasmussen. The partition
%   code is written by Matt Fig.
%
%   Author: Tarmo ?ij? <tarmo.aijo@aalto.fi>
%
%   Updated: 03/05/2013


if nargin < 3
    error('not enough parameters')
end

if length(time_vector) ~= length(data_vector)
    error('cell vectors time_vector and data_vector should have the same length')
elseif length(time_vector) < 2
    error('there should be at least two conditions')
end

scaling_factor = max(cellfun(@max,data_vector));
data_vector = cellfun(@(x) x/scaling_factor,data_vector,'UniformOutput',false);

result = ligap_calc(time_vector,data_vector,labels,scaling_factor);

end

function result = ligap_calc(timeVector,profiles,labels,scaling_factor)

% scale the time vector
time_scaling_factor = max(cellfun(@max,timeVector));
timeVector_orig = timeVector;
timeVector = cellfun(@(x) x/time_scaling_factor,timeVector,'UniformOutput',false);

% define the covariance function
covFunc = {'covSum',{'covNNone','covNoise'}};

% number of the different initial values
initValIter = 10;

% find the different models
combinations = partitions(length(profiles));
numCombinations = length(combinations);

% construct the labels
labelsFull = cell(numCombinations,1);
for combinationIdx=1:numCombinations
    str = [];
    for foo=1:length(combinations{combinationIdx})
        str = [str '{' [labels{combinations{combinationIdx}{foo}}] '}'];
    end
    labelsFull{combinationIdx} = str;
end

% go through the different models
margLLHs = zeros(numCombinations,1);
for combinationIdx=1:length(combinations)
    tmpVector = zeros(length(combinations{combinationIdx}),1);
    for foo=1:length(combinations{combinationIdx})
        tmp = cat(1,profiles{combinations{combinationIdx}{foo}});
        Y = tmp;
        timeVectorLong = cat(1,timeVector{combinations{combinationIdx}{foo}});
        tmpVector(foo) = calculateMargLLH(timeVectorLong,Y,covFunc,initValIter);
    end
    % calculate the product of the marginal likelihoods
    margLLHs(combinationIdx) = prod(tmpVector);
end
% uniform prior
priors = ones(length(margLLHs),1); 
priors = priors/sum(priors);
% posterior probabilities
probabilities = (priors.*margLLHs)/sum(priors.*margLLHs);

curves = cell(length(profiles),1);
for profileIdx=1:length(profiles)
    % points where the GPs are evaluated
    timeVectorTest = linspace(min(timeVector_orig{profileIdx}),max(timeVector_orig{profileIdx}),100)';
    [foo,bar] = calculateMargLLHDraw(timeVector{profileIdx},profiles{profileIdx},covFunc,initValIter,scaling_factor,timeVectorTest/max(timeVectorTest));
    curves{profileIdx}.mean = foo;
    curves{profileIdx}.std = bar;
    curves{profileIdx}.timepoints = timeVectorTest;
end

result.probabilities = probabilities;
result.models = labelsFull;
result.profiles = curves;

end

function z = calculateMargLLH(timeVector,profile,covFunc,initValIter)
tmpfX = Inf;
tmpX = zeros(3,1);
for idx = 1:initValIter
    % fit the model
    status = true;
    while status
        % sample a set of the initial values
        X = [0; -3; -3] + ([1; 1; 1]-[0; -3; -3]).*rand(3,1);
        status = false;
        try
            [logHyper,fX] = minimize(X, 'gpr', -1000, covFunc, timeVector, profile);
        catch all
            status = true;
        end
    end
    if tmpfX > fX(end)
        tmpX = logHyper;
        tmpfX = fX(end);
    end
end
z = exp(-tmpfX);
end

function [mu,S2] = calculateMargLLHDraw(timeVector,profile,covFunc,initValIter,scaling_factor,timeVectorTest)
tmpfX = Inf;
tmpX = zeros(3,1);
for idx = 1:initValIter
    % fit a model
    status = true;
    while status
        % sample a set of initial values
        X = [0; -3; -3] + ([1; 1; 1]-[0; -3; -3]).*rand(3,1);
        status = false;
        try
            [logHyper,fX] = minimize(X, 'gpr', -1000, covFunc, timeVector, profile);
        catch all
            status = true;
        end
    end
    if tmpfX > fX(end)
        tmpX = logHyper;
        tmpfX = fX(end);
    end
end
[mu,S2] = gpr(tmpX, covFunc, timeVector, profile, timeVectorTest);
mu = mu*scaling_factor;
S2 = (S2 - exp(2*tmpX(3)))*scaling_factor^2;
end

function [out1, out2] = gpr(logtheta, covfunc, x, y, xstar)
% gpr - Gaussian process regression, with a named covariance function. Two
% modes are possible: training and prediction: if no test data are given, the
% function returns minus the log likelihood and its partial derivatives with
% respect to the hyperparameters; this mode is used to fit the hyperparameters.
% If test data are given, then (marginal) Gaussian predictions are computed,
% whose mean and variance are returned. Note that in cases where the covariance
% function has noise contributions, the variance returned in S2 is for noisy
% test targets; if you want the variance of the noise-free latent function, you
% must substract the noise variance.
%
% usage: [nlml dnlml] = gpr(logtheta, covfunc, x, y)
%    or: [mu S2]  = gpr(logtheta, covfunc, x, y, xstar)
%
% where:
%
%   logtheta is a (column) vector of log hyperparameters
%   covfunc  is the covariance function
%   x        is a n by D matrix of training inputs
%   y        is a (column) vector (of size n) of targets
%   xstar    is a nn by D matrix of test inputs
%   nlml     is the returned value of the negative log marginal likelihood
%   dnlml    is a (column) vector of partial derivatives of the negative
%                 log marginal likelihood wrt each log hyperparameter
%   mu       is a (column) vector (of size nn) of prediced means
%   S2       is a (column) vector (of size nn) of predicted variances
%
% For more help on covariance functions, see "help covFunctions".
%
% (C) copyright 2006 by Carl Edward Rasmussen (2006-03-20).

if ischar(covfunc), covfunc = cellstr(covfunc); end % convert to cell if needed
[n, D] = size(x);
if eval(feval(covfunc{:})) ~= size(logtheta, 1)
  error('Error: Number of parameters do not agree with covariance function')
end

K = feval(covfunc{:}, logtheta, x);    % compute training set covariance matrix

L = chol(K)';                        % cholesky factorization of the covariance
alpha = solve_chol(L',y);

if nargin == 4 % if no test cases, compute the negative log marginal likelihood

  out1 = 0.5*y'*alpha + sum(log(diag(L))) + 0.5*n*log(2*pi);

  if nargout == 2               % ... and if requested, its partial derivatives
    out2 = zeros(size(logtheta));       % set the size of the derivative vector
    W = L'\(L\eye(n))-alpha*alpha';                % precompute for convenience
    for i = 1:length(out2)
      out2(i) = sum(sum(W.*feval(covfunc{:}, logtheta, x, i)))/2;
    end
  end

else                    % ... otherwise compute (marginal) test predictions ...

  [Kss, Kstar] = feval(covfunc{:}, logtheta, x, xstar);     %  test covariances

  out1 = Kstar' * alpha;                                      % predicted means

  if nargout == 2
    v = L\Kstar;
    out2 = Kss - sum(v.*v)';
  end  

end
end

function [A, B] = covNNone(loghyper, x, z)

% Neural network covariance function with a single parameter for the distance
% measure. The covariance function is parameterized as:
%
% k(x^p,x^q) = sf2 * asin(x^p'*P*x^q / sqrt[(1+x^p'*P*x^p)*(1+x^q'*P*x^q)])
%
% where the x^p and x^q vectors on the right hand side have an added extra bias
% entry with unit value. P is ell^-2 times the unit matrix and sf2 controls the
% signal variance. The hyperparameters are:
%
% loghyper = [ log(ell)
%              log(sqrt(sf2) ]
%
% For more help on design of covariance functions, try "help covFunctions".
%
% (C) Copyright 2006 by Carl Edward Rasmussen (2006-03-24)

if nargin == 0, A = '2'; return; end              % report number of parameters

persistent Q K;                 
[n D] = size(x);
ell = exp(loghyper(1)); em2 = ell^(-2);
sf2 = exp(2*loghyper(2));
x = x/ell;

if nargin == 2                                             % compute covariance
  Q = x*x';
  K = (em2+Q)./(sqrt(1+em2+diag(Q))*sqrt(1+em2+diag(Q)'));
  A = sf2*asin(K);                 
elseif nargout == 2                              % compute test set covariances
  z = z/ell; 
  A = sf2*asin((em2+sum(z.*z,2))./(1+em2+sum(z.*z,2)));
  B = sf2*asin((em2+x*z')./sqrt((1+em2+sum(x.*x,2))*(1+em2+sum(z.*z,2)')));
else                                                % compute derivative matrix
  % check for correct dimension of the previously calculated kernel matrix
  if any(size(Q)~=n)  
    Q = x*x';
  end
  % check for correct dimension of the previously calculated kernel matrix
  if any(size(K)~=n)  
    K = (em2+Q)./(sqrt(1+em2+diag(Q))*sqrt(1+em2+diag(Q)'));
  end
  if z == 1                                                   % first parameter
    v = (em2+sum(x.*x,2))./(1+em2+diag(Q));
    A = -2*sf2*((em2+Q)./(sqrt(1+em2+diag(Q))*sqrt(1+em2+diag(Q)'))- ...
                            K.*(repmat(v,1,n)+repmat(v',n,1))/2)./sqrt(1-K.^2);
    clear Q;
  else                                                       % second parameter
    A = 2*sf2*asin(K);
    clear K;
  end
end
end

function C = sq_dist(a, b, Q);
% sq_dist - a function to compute a matrix of all pairwise squared distances
% between two sets of vectors, stored in the columns of the two matrices, a
% (of size D by n) and b (of size D by m). If only a single argument is given
% or the second matrix is empty, the missing matrix is taken to be identical
% to the first.
%
% Special functionality: If an optional third matrix argument Q is given, it
% must be of size n by m, and in this case a vector of the traces of the
% product of Q' and the coordinatewise squared distances is returned.
%
% NOTE: The program code is written in the C language for efficiency and is
% contained in the file sq_dist.c, and should be compiled using matlabs mex
% facility. However, this file also contains a (less efficient) matlab
% implementation, supplied only as a help to people unfamiliar with mex. If
% the C code has been properly compiled and is avaiable, it automatically
% takes precendence over the matlab code in this file.
%
% Usage: C = sq_dist(a, b)
%    or: C = sq_dist(a)  or equiv.: C = sq_dist(a, [])
%    or: c = sq_dist(a, b, Q)
% where the b matrix may be empty.
%
% where a is of size D by n, b is of size D by m (or empty), C and Q are of
% size n by m and c is of size D by 1.
%
% Copyright (c) 2003, 2004, 2005 and 2006 Carl Edward Rasmussen.
% 2006-03-09.

if nargin < 1 | nargin > 3 | nargout > 1
  error('Wrong number of arguments.');
end

if nargin == 1 | isempty(b)                   % input arguments are taken to be
  b = a;                                   % identical if b is missing or empty
end 

[D, n] = size(a); 
[d, m] = size(b);
if d ~= D
  error('Error: column lengths must agree.');
end

if nargin < 3
  C = zeros(n,m);
  for d = 1:D
    C = C + (repmat(b(d,:), n, 1) - repmat(a(d,:)', 1, m)).^2;
  end
  % C = repmat(sum(a.*a)',1,m)+repmat(sum(b.*b),n,1)-2*a'*b could be used to 
  % replace the 3 lines above; it would be faster, but numerically less stable.
else
  if [n m] == size(Q)
    C = zeros(D,1);
    for d = 1:D
      C(d) = sum(sum((repmat(b(d,:), n, 1) - repmat(a(d,:)', 1, m)).^2.*Q));
    end
  else
    error('Third argument has wrong size.');
  end
end
end

function [A, B] = covSum(covfunc, logtheta, x, z)
% covSum - compose a covariance function as the sum of other covariance
% functions. This function doesn't actually compute very much on its own, it
% merely does some bookkeeping, and calls other covariance functions to do the
% actual work.
%
% For more help on design of covariance functions, try "help covFunctions".
%
% (C) Copyright 2006 by Carl Edward Rasmussen, 2006-03-20.

for i = 1:length(covfunc)                   % iterate over covariance functions
  f = covfunc(i);
  if iscell(f{:}), f = f{:}; end          % dereference cell array if necessary
  j(i) = cellstr(feval(f{:}));
end

if nargin == 1,                                   % report number of parameters
  A = char(j(1)); for i=2:length(covfunc), A = [A, '+', char(j(i))]; end
  return
end

[n, D] = size(x);

v = [];              % v vector indicates to which covariance parameters belong
for i = 1:length(covfunc), v = [v repmat(i, 1, eval(char(j(i))))]; end

switch nargin
case 3                                              % compute covariance matrix
  A = zeros(n, n);                       % allocate space for covariance matrix
  for i = 1:length(covfunc)                  % iteration over summand functions
    f = covfunc(i);
    if iscell(f{:}), f = f{:}; end        % dereference cell array if necessary
    A = A + feval(f{:}, logtheta(v==i), x);            % accumulate covariances
  end

case 4                      % compute derivative matrix or test set covariances
  if nargout == 2                                % compute test set cavariances
    A = zeros(size(z,1),1); B = zeros(size(x,1),size(z,1));    % allocate space
    for i = 1:length(covfunc)
      f = covfunc(i);
      if iscell(f{:}), f = f{:}; end      % dereference cell array if necessary
      [AA BB] = feval(f{:}, logtheta(v==i), x, z);   % compute test covariances
      A = A + AA; B = B + BB;                                  % and accumulate
    end
  else                                            % compute derivative matrices
    i = v(z);                                       % which covariance function
    j = sum(v(1:z)==i);                    % which parameter in that covariance
    f = covfunc(i);
    if iscell(f{:}), f = f{:}; end        % dereference cell array if necessary
    A = feval(f{:}, logtheta(v==i), x, j);                 % compute derivative
  end

end
end

function [A, B] = covNoise(logtheta, x, z)
% Independent covariance function, ie "white noise", with specified variance.
% The covariance function is specified as:
%
% k(x^p,x^q) = s2 * \delta(p,q)
%
% where s2 is the noise variance and \delta(p,q) is a Kronecker delta function
% which is 1 iff p=q and zero otherwise. The hyperparameter is
%
% logtheta = [ log(sqrt(s2)) ]
%
% For more help on design of covariance functions, try "help covFunctions".
%
% (C) Copyright 2006 by Carl Edward Rasmussen, 2006-03-24.

if nargin == 0, A = '1'; return; end              % report number of parameters

s2 = exp(2*logtheta);                                          % noise variance

if nargin == 2                                      % compute covariance matrix
  A = s2*eye(size(x,1));
elseif nargout == 2                              % compute test set covariances
  A = s2;
  B = 0;                               % zeros cross covariance by independence
else                                                % compute derivative matrix
  A = 2*s2*eye(size(x,1));
end
end

function [X, fX, i] = minimize(X, f, length, varargin)
% Minimize a differentiable multivariate function. 
%
% Usage: [X, fX, i] = minimize(X, f, length, P1, P2, P3, ... )
%
% where the starting point is given by "X" (D by 1), and the function named in
% the string "f", must return a function value and a vector of partial
% derivatives of f wrt X, the "length" gives the length of the run: if it is
% positive, it gives the maximum number of line searches, if negative its
% absolute gives the maximum allowed number of function evaluations. You can
% (optionally) give "length" a second component, which will indicate the
% reduction in function value to be expected in the first line-search (defaults
% to 1.0). The parameters P1, P2, P3, ... are passed on to the function f.
%
% The function returns when either its length is up, or if no further progress
% can be made (ie, we are at a (local) minimum, or so close that due to
% numerical problems, we cannot get any closer). NOTE: If the function
% terminates within a few iterations, it could be an indication that the
% function values and derivatives are not consistent (ie, there may be a bug in
% the implementation of your "f" function). The function returns the found
% solution "X", a vector of function values "fX" indicating the progress made
% and "i" the number of iterations (line searches or function evaluations,
% depending on the sign of "length") used.
%
% The Polack-Ribiere flavour of conjugate gradients is used to compute search
% directions, and a line search using quadratic and cubic polynomial
% approximations and the Wolfe-Powell stopping criteria is used together with
% the slope ratio method for guessing initial step sizes. Additionally a bunch
% of checks are made to make sure that exploration is taking place and that
% extrapolation will not be unboundedly large.
%
% See also: checkgrad 
%
% Copyright (C) 2001 - 2006 by Carl Edward Rasmussen (2006-09-08).

INT = 0.1;    % don't reevaluate within 0.1 of the limit of the current bracket
EXT = 3.0;                  % extrapolate maximum 3 times the current step-size
MAX = 20;                         % max 20 function evaluations per line search
RATIO = 10;                                       % maximum allowed slope ratio
SIG = 0.1; RHO = SIG/2; % SIG and RHO are the constants controlling the Wolfe-
% Powell conditions. SIG is the maximum allowed absolute ratio between
% previous and new slopes (derivatives in the search direction), thus setting
% SIG to low (positive) values forces higher precision in the line-searches.
% RHO is the minimum allowed fraction of the expected (from the slope at the
% initial point in the linesearch). Constants must satisfy 0 < RHO < SIG < 1.
% Tuning of SIG (depending on the nature of the function to be optimized) may
% speed up the minimization; it is probably not worth playing much with RHO.

% The code falls naturally into 3 parts, after the initial line search is
% started in the direction of steepest descent. 1) we first enter a while loop
% which uses point 1 (p1) and (p2) to compute an extrapolation (p3), until we
% have extrapolated far enough (Wolfe-Powell conditions). 2) if necessary, we
% enter the second loop which takes p2, p3 and p4 chooses the subinterval
% containing a (local) minimum, and interpolates it, unil an acceptable point
% is found (Wolfe-Powell conditions). Note, that points are always maintained
% in order p0 <= p1 <= p2 < p3 < p4. 3) compute a new search direction using
% conjugate gradients (Polack-Ribiere flavour), or revert to steepest if there
% was a problem in the previous line-search. Return the best value so far, if
% two consecutive line-searches fail, or whenever we run out of function
% evaluations or line-searches. During extrapolation, the "f" function may fail
% either with an error or returning Nan or Inf, and minimize should handle this
% gracefully.

if max(size(length)) == 2, red=length(2); length=length(1); else red=1; end
if length>0, S='Linesearch'; else S='Function evaluation'; end 

i = 0;                                            % zero the run length counter
ls_failed = 0;                             % no previous line search has failed
[f0 df0] = feval(f, X, varargin{:});          % get function value and gradient
fX = f0;
i = i + (length<0);                                            % count epochs?!
s = -df0; d0 = -s'*s;           % initial search direction (steepest) and slope
x3 = red/(1-d0);                                  % initial step is red/(|s|+1)

while i < abs(length)                                      % while not finished
  i = i + (length>0);                                      % count iterations?!

  X0 = X; F0 = f0; dF0 = df0;                   % make a copy of current values
  if length>0, M = MAX; else M = min(MAX, -length-i); end

  while 1                             % keep extrapolating as long as necessary
    x2 = 0; f2 = f0; d2 = d0; f3 = f0; df3 = df0;
    success = 0;
    while ~success && M > 0
      try
        M = M - 1; i = i + (length<0);                         % count epochs?!
        [f3 df3] = feval(f, X+x3*s, varargin{:});
        if isnan(f3) || isinf(f3) || any(isnan(df3)+isinf(df3)), error(''), end
        success = 1;
      catch                                % catch any error which occured in f
        x3 = (x2+x3)/2;                                  % bisect and try again
      end
    end
    if f3 < F0, X0 = X+x3*s; F0 = f3; dF0 = df3; end         % keep best values
    d3 = df3'*s;                                                    % new slope
    if d3 > SIG*d0 || f3 > f0+x3*RHO*d0 || M == 0  % are we done extrapolating?
      break
    end
    x1 = x2; f1 = f2; d1 = d2;                        % move point 2 to point 1
    x2 = x3; f2 = f3; d2 = d3;                        % move point 3 to point 2
    A = 6*(f1-f2)+3*(d2+d1)*(x2-x1);                 % make cubic extrapolation
    B = 3*(f2-f1)-(2*d1+d2)*(x2-x1);
    x3 = x1-d1*(x2-x1)^2/(B+sqrt(B*B-A*d1*(x2-x1))); % num. error possible, ok!
    if ~isreal(x3) || isnan(x3) || isinf(x3) || x3 < 0 % num prob | wrong sign?
      x3 = x2*EXT;                                 % extrapolate maximum amount
    elseif x3 > x2*EXT                  % new point beyond extrapolation limit?
      x3 = x2*EXT;                                 % extrapolate maximum amount
    elseif x3 < x2+INT*(x2-x1)         % new point too close to previous point?
      x3 = x2+INT*(x2-x1);
    end
  end                                                       % end extrapolation

  while (abs(d3) > -SIG*d0 || f3 > f0+x3*RHO*d0) && M > 0  % keep interpolating
    if d3 > 0 || f3 > f0+x3*RHO*d0                         % choose subinterval
      x4 = x3; f4 = f3; d4 = d3;                      % move point 3 to point 4
    else
      x2 = x3; f2 = f3; d2 = d3;                      % move point 3 to point 2
    end
    if f4 > f0           
      x3 = x2-(0.5*d2*(x4-x2)^2)/(f4-f2-d2*(x4-x2));  % quadratic interpolation
    else
      A = 6*(f2-f4)/(x4-x2)+3*(d4+d2);                    % cubic interpolation
      B = 3*(f4-f2)-(2*d2+d4)*(x4-x2);
      x3 = x2+(sqrt(B*B-A*d2*(x4-x2)^2)-B)/A;        % num. error possible, ok!
    end
    if isnan(x3) || isinf(x3)
      x3 = (x2+x4)/2;               % if we had a numerical problem then bisect
    end
    x3 = max(min(x3, x4-INT*(x4-x2)),x2+INT*(x4-x2));  % don't accept too close
    [f3 df3] = feval(f, X+x3*s, varargin{:});
    if f3 < F0, X0 = X+x3*s; F0 = f3; dF0 = df3; end         % keep best values
    M = M - 1; i = i + (length<0);                             % count epochs?!
    d3 = df3'*s;                                                    % new slope
  end                                                       % end interpolation

  if abs(d3) < -SIG*d0 && f3 < f0+x3*RHO*d0          % if line search succeeded
    X = X+x3*s; f0 = f3; fX = [fX' f0]';                     % update variables
    %fprintf('%s %6i;  Value %4.6e\r', S, i, f0);
    s = (df3'*df3-df0'*df3)/(df0'*df0)*s - df3;   % Polack-Ribiere CG direction
    df0 = df3;                                               % swap derivatives
    d3 = d0; d0 = df0'*s;
    if d0 > 0                                      % new slope must be negative
      s = -df0; d0 = -s'*s;                  % otherwise use steepest direction
    end
    x3 = x3 * min(RATIO, d3/(d0-realmin));          % slope ratio but max RATIO
    ls_failed = 0;                              % this line search did not fail
  else
    X = X0; f0 = F0; df0 = dF0;                     % restore best point so far
    if ls_failed || i > abs(length)         % line search failed twice in a row
      break;                             % or we ran out of time, so we give up
    end
    s = -df0; d0 = -s'*s;                                        % try steepest
    x3 = 1/(1-d0);                     
    ls_failed = 1;                                    % this line search failed
  end
end
%fprintf('\n');
end

function x = solve_chol(A, B)
% solve_chol - solve linear equations from the Cholesky factorization.
% Solve A*X = B for X, where A is square, symmetric, positive definite. The
% input to the function is R the Cholesky decomposition of A and the matrix B.
% Example: X = solve_chol(chol(A),B);
%
% NOTE: The program code is written in the C language for efficiency and is
% contained in the file solve_chol.c, and should be compiled using matlabs mex
% facility. However, this file also contains a (less efficient) matlab
% implementation, supplied only as a help to people unfamiliar with mex. If
% the C code has been properly compiled and is avaiable, it automatically
% takes precendence over the matlab code in this file.
%
% Copyright (c) 2004, 2005, 2006 by Carl Edward Rasmussen. 2006-02-08.

if nargin ~= 2 | nargout > 1
  error('Wrong number of arguments.');
end

if size(A,1) ~= size(A,2) | size(A,1) ~= size(B,1)
  error('Wrong sizes of matrix arguments.');
end

x = A\(A'\B);
end

function C = partitions(M,K)
%PARTITIONS Find all partitions of a set.
% C = PARTITIONS(N), for scalar N, returns a cell array, wherein each cell 
% has a various number of other cells representing the partitions of the  
% set of numbers {1,2,3,...N}.  The length of C is the Bell number for N.
% C = PARTITIONS(N), for vector N, returns the partitions of the vector
% elements treated as members of a set.
% C = PARTITIONS(N), for cell N, returns the partitions of the cell
% elements treated as members of a set.
% PARTITIONS(N,K) returns only the partitions of which have K members.  
%
% K must be less than or equal to N for scalar N, or length(N).
%
%
% Examples:
%
%     C = partitions(4); % Find the partitions of set {1 2 3 4}
%     home  % Makes for nice display for small enough N.  Else use clc.
%     partdisp(C) % File accompanying PARTITIONS
%
%     C = partitions({'peanut','butter','yummy'});
%     home
%     partdisp(C)
%
%     C = partitions([5 7 8 50],3); % partitions of set {5 7 8 50}
%     home
%     partdisp(C,3)
%
%     C = partitions(['a','b','c','d']); % for longer chars, use {}, not []
%     home
%     partdisp(C)
%
% Class support for inputs N,K:
%      float: double, single, char (N only)
% 
% See also, nchoosek, perms, npermutek (On the FEX)
%
% Author: Matt Fig,   
% Contact: popkenai@yahoo.com
% Date: 5/17/2009

Kflg = false;  % No K passed in.
Cflg = false;  % A set not passed in.

if length(M)>1 
    N = length(M); % User passed {'here','there','everywhere'} for example.
    Cflg = true;
else
    if iscell(M)  % User passed {2} for example.
        error('Set arguments must have more than one element.')
    end
    
    N = M;
end

if nargin>1
    Kflg = true;  % Used in while loop below, K passed in.
    S = stirling(N,K);  % The number of partitions, Stirling number.
    C = cell(S,min(1,K)); % Main Cell.
    cnt = 1;  % Start the while loop counter.
else
    K = 0;  % Since user doesn't want only certain partitions.  
    S = ceil(sum((1:2*N).^N./cumprod(1:2*N))/exp(1)); % Bell number.
    C = cell(S,1); % Main Cell.
    
    if Cflg
        C{1} = {M};
    else
        C{1} = {1:N}; % First one is easy.
    end
    
    cnt = 2; % Start the while loop counter.
end

if K~=ceil(K) || numel(K)~=1  
    error('Second argument must be an integer of length 1.  See help.')
end

if N<0 || K>N || K<0
   error('Arguments must be greater than zero, and K cannot exceed N.') 
end

if N==0
    C = {{}}; % Easy case.
    return
end

if K==1
    if Cflg
        C{1} = {M};  % Easy case.
    else
        C{1} = {1:N};  % Easy case.
    end
    return
end

if Cflg
    NV = M;
else
    NV = 1:M; % Vector: base partition, RGF indexes into this guy.
end

NU = 1; % Number of unique indices in current partition.
stp1 = N; % Controls assigning of indices.
RGF = ones(1,N); % Holds the indexes.
BLD = cell(1,N); % Smaller cell array will be used in creating larger.

while cnt<=S
    idx1 = N; % Index into RGF.
    stp2 = RGF(idx1); % Works with stp1.

    while stp1(stp2)==1
        RGF(idx1) = 1;  % Assign value to RGF.
        idx1 = idx1 - 1; % Need to increment idx1 for translation below.
        stp2 = RGF(idx1); % And set this guy for stp1 assign below.
    end

    NU = NU + idx1 - N;  % Get provisional number of unique vals.
    stp1(1) = stp1(1) + N - idx1;

    if stp2==NU % Increment the number of unique elements.
        NU = NU +1;
        stp1(NU) = 0;
    end

    RGF(idx1) = stp2 + 1;  % Increment this position.
    stp1(stp2) = stp1(stp2) - 1;  % Translate indices of these two. 
    stp1(stp2+1) = stp1(stp2+1) + 1; % Re-assign stopper.

    if NU==(~Kflg * NU + Kflg * K)
    % We could use:  C{cnt} = accumarray(RGF',NV',[],@(x) {x});   (SLOW!!)
    % or the next lines to C{cnt} = TMP;
        TMP = BLD(1:NU); % Build subcell of correct size.
        TMP{1} = NV(RGF==1); % These first two are always here.... no loop.
        TMP{2} = NV(RGF==2);

        for ii = 3:NU % Build the rest of cell array, if any.
            TMP{ii} = NV(RGF==ii);
        end

        C{cnt} = TMP;  % Assign cell at jj. 
        cnt = cnt + 1;
    end
end
end

function S = stirling(N,K)
% Calculate the Stirling number of the second kind. Subfunc to partitions.

for ii = K:-1:0
    S(ii+1) = (-1)^ii * prod(1:K) / (prod(1:(K-ii)) *...
              (prod(1:ii)))*(K - ii)^N;  %#ok
end

S = 1/prod(1:K) * sum(S);
end
