1;

## Fit a monotone non-decreasing curve to the data points (Xi, Yi).
##
## The fitted curve is (Xi, Hi).
##
## Assume the Xi are increasing, and therefore the constraint is that the
## Hi are nondecreasing.
##
## The loss function is:
##        sum_i (Yi - Hi)^2 - \lambda * sum_i log(H(i+1) - Hi)



## Calculate the objective function and gradient.
function [V,G] = fun(Y, X, H, lambda)

  ## The number of data points.
  n = length(H);

  ## Test for a nonfeasible point.
  if (any(H(2:n) <= H(1:n-1)))
    V = Inf;
    G = [];
    return;
  endif

  ## The loss function value.
  V = sumsq(Y-H) - lambda*sum(log(H(2:n)-H(1:n-1)));

  ## The loss function gradient.
  G = -2*(Y-H);
  G(2:n) = G(2:n) - lambda./(H(2:n) - H(1:n-1));
  G(1:n-1) = G(1:n-1) + lambda./(H(2:n) - H(1:n-1));

endfunction


##
## Simulate data according to the model.
##

## Sample size.
n = 50;

## Abscissas.
X = sort(5*rand(n,1));

## True ordinates.
YT = sort(log(5*rand(n,1)));

## Observed ordinates.
Y = YT + 0.1*randn(size(YT));



##
## Optimization.
##

## A starting value.
beta = cov(Y, X) / var(X);
alpha = mean(Y) - beta*mean(X);
if (beta < 0)
  beta = 0.1;
endif
H = alpha + beta*X;

## The penalty function value.
for lambda = [1e-1, 1e-2, 1e-3, 1e-4]

  ## Initialize the gradient and search direction.
  [u,g] = fun(Y, X, H, lambda);
  h = g;

  ## Iterations loop for each method.
  for it=1:1000
    
    ## Reset the search direction each time we cycle through one set of
    ## conjugate directions.
    if (rem(it-1, n) == 0)
      [u,h] = fun(Y, X, H, lambda);
      g = h;
    endif
    
    ##
    ## Line search in the direction of h.
    ##
    
    ## Get an initial bracket.
    L = [-1,0,1];
    F = [fun(Y, X, H-L(1)*h, lambda), fun(Y, X, H, lambda), ...
	 fun(Y, X, H-L(3)*h, lambda)];
    while (F(2) > min(F(1), F(3)))
      L = L*2;
      F(1) = fun(Y, X, H-L(1)*h, lambda);
      F(3) = fun(Y, X, H-L(3)*h, lambda);
    endwhile
    
    ## Bisection.
    while (L(3)-L(1) > 1e-8)
      
      if (L(3)-L(2) > L(2)-L(1))
	L1 = (L(2)+L(3)) / 2;
	F1 = fun(Y, X, H-L1*h, lambda);
	if (F1 > F(2))
 	  L = [L(1), L(2), L1];
 	  F = [F(1), F(2), F1];
	else
 	  L = [L(2), L1, L(3)];
 	  F = [F(2), F1, F(3)];
	endif
      else
	L1 = (L(1)+L(2)) / 2;
	F1 = fun(Y, X, H-L1*h, lambda);
	if (F1 > F(2))
 	  L = [L1, L(2), L(3)];
 	  F = [F1, F(2), F(3)];
	else
 	  L = [L(1), L1, L(2)];
 	  F = [F(1), F1, F(2)];
	endif
      endif
      
    endwhile
    
    ## Move the point.
    H = H - L(2)*h;
    
    ##
    ## Update the search direction.
    ##
    
    ## Get the new gradient.
    [u, g1] = fun(Y, X, H, lambda);
    
    ## Test for convergence.
    if (norm(g1) < 1e-4)
      break;
    endif
    
    ## Update the Fletcher-Powell search direction.
    gamma = sumsq(g1) / sumsq(g);
    h = g1 + gamma*h;
    g = g1;
    
#    gset key below;
#    A = [X, Y, YT, H];
#    gplot A(:,[1,2]) title "Observed" with points, ...
#	A(:,[1,3]) title "True" with points, ...
#	A(:,[1,4]) title "Fitted" with lines; 
    
  disp(norm(g1));
  endfor 
endfor

if (norm(g1) > 1e-4)
  disp("Warning: did not converge\n");
endif