if nargin < 3, A = [];
                     end, end, end, end, end, end, end, end

% Set up constant strings
medium =  'medium-scale: active-set';
large = 'large-scale';

if nargout > 4
   computeLambda = 1;
else 
   computeLambda = 0;
end

% Options setup
largescale = isequal(optimget(options,'LargeScale',defaultopt,'fast'),'on');
diagnostics = isequal(optimget(options,'Diagnostics',defaultopt,'fast'),'on');
switch optimget(options,'Display',defaultopt,'fast')
case {'off', 'none'}
   verbosity = 0;
case 'iter'
   verbosity = 2;
case 'final'
   verbosity = 1;
case 'testing'
   verbosity = Inf;
otherwise
   verbosity = 1;
end
mtxmpy = optimget(options,'HessMult',defaultopt,'fast');
% check if name clash
if isequal(mtxmpy,'hmult')
   warnstr = sprintf('%s\n%s\n%s\n', ...
            'Potential function name clash with a Toolbox helper function:',...
            'Use a name besides ''hmult'' for your HessMult function to',...
            'avoid errors or unexpected results.');
   warning(warnstr)
end

% Set the constraints up: defaults and check size
[nineqcstr,numberOfVariablesineq]=size(A);
[neqcstr,numberOfVariableseq]=size(Aeq);
if isa(H,'double') & ( isempty(mtxmpy) )
   lengthH = length(H);
else % HessMult in effect, so H can be anything
   lengthH = 0;
end

numberOfVariables = ...
    max([length(f),lengthH,numberOfVariablesineq,numberOfVariableseq]); % In case A or Aeq is empty
ncstr = nineqcstr + neqcstr;

if isempty(f), f=zeros(numberOfVariables,1); end
if isempty(A), A=zeros(0,numberOfVariables); end
if isempty(B), B=zeros(0,1); end
if isempty(Aeq), Aeq=zeros(0,numberOfVariables); end
if isempty(Beq), Beq=zeros(0,1); end

% Expect vectors
f=f(:);
B=B(:);
Beq=Beq(:);

if ~isequal(length(B),nineqcstr)
    error('The number of rows in A must be the same as the length of b.')
elseif ~isequal(length(Beq),neqcstr)
    error('The number of rows in Aeq must be the same as the length of beq.')
elseif ~isequal(length(f),numberOfVariablesineq) & ~isempty(A)
    error('The number of columns in A must be the same as the length of f.')
elseif ~isequal(length(f),numberOfVariableseq) & ~isempty(Aeq)
    error('The number of columns in Aeq must be the same as the length of f.')
end

[X0,lb,ub,msg] = checkbounds(X0,lb,ub,numberOfVariables);
if ~isempty(msg)
   exitflag = -1;
   output = []; X=X0; fval = []; lambda = [];
   if verbosity > 0
      disp(msg)
   end
   return
end

caller = 'quadprog';
% Check out H and make sure it isn't empty or all zeros
if isa(H,'double') & isempty(mtxmpy)
   if  norm(H,'inf')==0 | isempty(H)
      H=[]; 
      % Really a lp problem
      caller = 'linprog';
      warning('Hessian is empty or all zero; calling LINPROG');
      [X,fval,exitflag,output,lambda]=linprog(f,A,B,Aeq,Beq,lb,ub,X0,options);
      return
   else
      % Make sure it is symmetric
      if norm(H-H',inf) > eps
         if verbosity > -1
            warning('Your Hessian is not symmetric.  Resetting H=(H+H'')/2')
         end
         H = (H+H')*0.5;
      end
   end
end

% Use large-scale algorithm or not?
% Determine which algorithm and make sure problem matches.

%    If any inequalities, 
%    or both equalities and bounds, 
%    or more equalities than variables,
%    or no equalities and no bounds and no inequalities
%    or asked for active set (~largescale) then call qpsub
if ( (nineqcstr > 0) | ...
      ( neqcstr > 0 & (sum(~isinf(ub))>0 | sum(~isinf(lb)) > 0)) | ...
      (neqcstr > numberOfVariables) | ...
      (neqcstr==0 & nineqcstr==0 & ... 
          all(eq(ub, inf)) & all(eq(lb, -inf))) | ...  % unconstrained
      ~largescale)
   % (has linear inequalites  OR both equalities and bounds) OR 
   % ~largescale, then call active-set code
   output.algorithm = medium;
%    if largescale  & ...
%          (  issparse(H)  | issparse(A) | issparse(Aeq) )% asked for sparse
%       warnstr = sprintf('%s\n%s\n', ...
%          'This problem formulation not yet available for sparse matrices.',...
%          'Converting to full matrices and switching to medium-scale method.');
%       warning(warnstr);
%    elseif largescale % and didn't ask for sparse
%          warning(['Large-scale method does not currently solve this problem formulation,',...
%                sprintf('\n'), 'switching to medium-scale method.',sprintf('\n')])
%       
%    end
   if ~isa(H,'double') | ( ~isempty(mtxmpy) ) 
      error('H must be specified explicitly for medium-scale algorithm: cannot use HessMult option.');
   end
   H = full(H); A = full(A); Aeq = full(Aeq);
else % call sqpmin when just bounds or just equalities
   output.algorithm = large;
   if isempty(mtxmpy) 
     H = sparse(H);
   end
   A = sparse(A); Aeq = sparse(Aeq);
end

if diagnostics 
   % Do diagnostics on information so far
   gradflag = []; hessflag = []; line_search=[];
   constflag = 0; gradconstflag = 0; non_eq=0;non_ineq=0;
   lin_eq=size(Aeq,1); lin_ineq=size(A,1); XOUT=ones(numberOfVariables,1);
   funfcn{1} = [];ff=[]; GRAD=[];HESS=[];
   confcn{1}=[];c=[];ceq=[];cGRAD=[];ceqGRAD=[];
   msg = diagnose('quadprog',output,gradflag,hessflag,constflag,gradconstflag,...
      line_search,options,defaultopt,XOUT,non_eq,...
      non_ineq,lin_eq,lin_ineq,lb,ub,funfcn,confcn,ff,GRAD,HESS,c,ceq,cGRAD,ceqGRAD);
end

% if any inequalities, or both equalities and bounds, or more equalities than bounds,
%    or asked for active set (~largescale) then call qpsub
if isequal(output.algorithm, medium)
   if isempty(X0), 
      X0=zeros(numberOfVariables,1); 
   end
   [X,lambdaqp,exitflag,output]= ...
      qpsub(H,f,[Aeq;A],[Beq;B],lb,ub,X0,neqcstr,...
      verbosity,caller,ncstr,numberOfVariables,options,defaultopt); 
   output.algorithm = medium; % have to reset since call to qpsub obliterates
   
elseif isequal(output.algorithm,large)  % largescale: call sqpmin when just bounds or just equalities
   [X,fval,output,exitflag,lambda]=...
      sqpmin(f,H,X0,Aeq,Beq,lb,ub,verbosity,options,defaultopt,computeLambda,varargin{:});
   
   if exitflag == -2  % Problem not handled by sqpmin at this time
      if largescale  & ( issparse(H) | issparse(A) | issparse(Aeq) )% asked for sparse
         warnstr = sprintf('%s\n%s\n', ...
         'This problem formulation not yet available for sparse matrices.',...
         'Converting to full matrices and switching to medium-scale method.');
         warning(warnstr);
      elseif largescale
         warning(['Large-scale method does not currently solve this problem formulation,',...
               sprintf('\n'), 'switching to medium-scale method.',sprintf('\n')])
      end
      
      if isempty(X0), 
         X0=zeros(numberOfVariables,1); 
      end
      output.algorithm = medium;
      if ~isa(H,'double') | ( ~isempty(mtxmpy)  ) 
        error('H must be specified explicitly for medium-scale algorithm: cannot use HessMult option.');
      end
      H = full(H); A = full(A); Aeq = full(Aeq);
      
      [X,lambdaqp,exitflag,output]= ...
         qpsub(H,f,[Aeq;A],[Beq;B],lb,ub,X0,neqcstr,...
         verbosity,caller,ncstr,numberOfVariables,options,defaultopt);
      output.algorithm = medium; % have to reset since call to qpsub obliterates
   end
end


if isequal(output.algorithm , medium)
   fval = 0.5*X'*(H*X)+f'*X; 
   llb = length(lb); 
   lub = length(ub);
   lambda.lower = zeros(llb,1);
   lambda.upper = zeros(lub,1);
   arglb = ~isinf(lb); lenarglb = nnz(arglb);
   argub = ~isinf(ub); lenargub = nnz(argub);
   lambda.eqlin = lambdaqp(1:neqcstr,1);
   lambda.ineqlin = lambdaqp(neqcstr+1:neqcstr+nineqcstr,1);
   lambda.lower(arglb) = lambdaqp(neqcstr+nineqcstr+1:neqcstr+nineqcstr+lenarglb);
   lambda.upper(argub) = lambdaqp(neqcstr+nineqcstr+lenarglb+1: ...
                                  neqcstr+nineqcstr+lenarglb+lenargub);
   
   output.firstorderopt=[];
   output.cgiterations =[];
   
   if verbosity > 0
      if ( exitflag ==1 )
         disp('Optimization terminated successfully.');   
      end
      if ( exitflag == 2)
         % do some sort of check here to see how unreliable
         disp('Optimization completed.'); 
      end
      if (exitflag ==0)
         disp('Maximum number of iterations exceeded;')
         disp('   increase options.MaxIter')
      end
      
   end
end


function errstring = consist(model, type, inputs, outputs)

errstring = '';

% If type string is not empty
if ~isempty(type)
  % First check that model has type field
  if ~isfield(model, 'type')
    errstring = 'Data structure does not contain type field';
    return
  end
  % Check that model has the correct type
  s = model.type;
  if ~strcmp(s, type)
    errstring = ['Model type ''', s, ''' does not match expected type ''',...
	type, ''''];
    return
  end
end

% If inputs are present, check that they have correct dimension
if nargin > 2
  if ~isfield(model, 'nin')
    errstring = 'Data structure does not contain nin field';
    return
  end

  data_nin = size(inputs, 2);
  if model.nin ~= data_nin
    errstring = ['Dimension of inputs ', num2str(data_nin), ...
	' does not match number of model inputs ', num2str(model.nin)];
    return
  end
end

% If outputs are present, check that they have correct dimension
if nargin > 3
  if ~isfield(model, 'nout')
    errstring = 'Data structure does not conatin nout field';
    return
  end
  data_nout = size(outputs, 2);
  if model.nout ~= data_nout
    errstring = ['Dimension of outputs ', num2str(data_nout), ...
	' does not match number of model outputs ', num2str(model.nout)];
    return
  end

% Also check that number of data points in inputs and outputs is the same
  num_in = size(inputs, 1);
  num_out = size(outputs, 1);
  if num_in ~= num_out
    errstring = ['Number of input patterns ', num2str(num_in), ...
	' does not match number of output patterns ', num2str(num_out)];
    return
  end
end

function net = svm(nin, kernel, kernelpar, C, use2norm, qpsolver, qpsize)

%

% 

if nargin < 7,
  qpsize = 50;
end
if nargin < 6,
  qpsolver = '';
end
if nargin < 5,
  use2norm = 0;
end
if nargin < 4,
  C = 1;
end
if nargin < 3,
  kernelpar = [];
end

net.type = 'svm';
net.nin = nin;
net.nout = 1;
net.kernel = kernel;
net.kernelpar = kernelpar;
net.c = C;
net.use2norm = use2norm;

net.nbexamples = 0;
net.alpha = [];
net.svcoeff = [];
net.sv = [];
net.svind = [];
net.bias = [];
net.normalw = [];

net.qpsolver = qpsolver;
net.qpsize = qpsize;
net.alphatol = 1e-2;
net.kkttol = 5e-2;
net.chunksize = 500;
%     'chunksize' = Large matrix operations (for example when evaluating
%       the kernel functions) are split up into submatrices with maximum
%       size [NET.chunksize, NET.chunksize]. Default value: 500
net.recompute = Inf;
%     'recompute' = During training, the SVM outputs are updated
%       iteratively. After NET.recompute iterations the SVM outputs are
%       built again from scratch. Lower this when high precision is required.


function [net, CVErr, paramSeq] = svmcv(net, X, Y, range, step, nfold, Xv, Yv, dodisplay)

% 

% Check arguments for consistency
errstring = consist(net, 'svm', X, Y);
if ~isempty(errstring);
  error(errstring);
end
if nargin<9,
  dodisplay = 1;
end
if nargin<8,
  Xv = [];
end
if nargin<7,
  Yv = [];
end
if nargin<6,
  nfold = 10;
end
if (~isempty(Xv)) & (~isempty(Yv)),
  errstring = consist(net, 'svm', Xv, Yv);
  if ~isempty(errstring);
    error(errstring);
  end
  if (nfold~=1),
    error('Input parameters XV and YV may only be used with NFOLD==1');
  end
end
if nargin<5,
  step = 0;
end

range = range(:)';
N = size(X, 1);
if N<nfold,
  error('At least NFOLD (default 10) training examples must be given');
end

if (length(range)>2) | isempty(step),
  % If range parameter has more than only min/max entries: Use this as
  % the sequence of parameters to test
  paramSeq = range;
else
  paramSeq = [];
  switch net.kernel
    case 'rbf'
      if step==0,
        step = sqrt(2);
      end
    % Multiplicative update, step size < 1 : start with max value
    if abs(step)<1,
      param = max(range);
      while (param>=min(range)),
        paramSeq = [paramSeq param];
        param = param*abs(step);
      end
    else
      % Multiplicative update, step size > 1 : start with min value
      param = min(range);