function [x, finalStats] = simple_mlp_optimCGP(initX, uFunc, uPar, ...
				 stepsize, maxIter, fgoal, gnsgoal)
% A "conjugate gradient parabolic thingy" method for training an MLP.
% There is no line search, only a parabolic fit. If there is no descent,
% try half the length. Very heuristic, but seems to work quite well.
%
% uFunc = "user function", in our case MLP evaluator, but could be anything
% uPar  = "user parameters", in our case training data, but could be
%         anything

nFC = 0;
x = initX; 
step = stepsize;

% Initially set everything to the negative gradient.
[f, grad2] = uFunc(x, uPar);
r = -grad2; rp = r; d = r;

for iter=1:maxIter
  statstr = ' ';

  %fprintf(1, 'f %.17f, gradns %.17f\n',f,grad2'*grad2);
  % Stop if any goal is reached
  if f <= fgoal, break, end;
  if (grad2' * grad2) <= gnsgoal, break, end;
  
  % Evaluate at a point further off in direction d.
  x2 = x + step * d;
  f2 = uFunc(x2, uPar);       % (Gradient not needed here.)

  % Find parameters for parabola satisfying f(x1), grad(x1), f(x2).
  b = -rp' * d;
  a = (f2 -b*step - f) / step^2;
  
  % Decide step size to reach bottom of the parabola.
  stepB = -b / max(2*a,1e-8);
  
  % Try step into the supposed directional minimum.
  x2 = x + stepB * d;
  [f2, grad2] = uFunc(x2, uPar);

  % Reality calling; it is possibly off. Then backstep.
  while f2>f && stepB^2 > 1e-7
    statstr = 'r';
    stepB = stepB / 2;
    x2 = x + stepB * d;
    [f2, grad2] = uFunc(x2, uPar);
    nFC = nFC + 1;
  end

  % If we had to backstep, try a smaller fit interval next time.
  if strcmp(statstr,'r'), step = max(0.6*step, 1e-6); end

  r = -grad2;

  % Update beta and search direction as in the Polak-Ribiere scheme.
  beta = sum(r .* (r - rp));
  beta = max(beta/(rp' * rp), 0);
  d = r + beta*d;

  % Update the current location and previous gradient information.
  x = x2; f = f2; rp = r;

%  disp(sprintf('%d/%d%s: f=%.17f  gns=%.17f  sz=%.17f', ...
%	       iter, maxIter, statstr, f, grad' * grad, step));

  step = step * 1.01;  % Always increase the step size a bit. This
                       % way it makes a heuristic, greedy wobble..
end
%disp(sprintf('%d/%d: f=%.17f  gns=%.17f\n', iter, maxIter, f, grad2' * grad2));

finalStats.outerIters = iter;
