// kelley.sce

//----------------------------------
// DATA
//----------------------------------


// Convex function to be minimized
function [y,z]=oracle(x)
  y=x.^2;
  z=2.*x;
endfunction

// Interval on which we are working
x_min=-10.1;
x_max=12.3;

// Number of iterations of the algorithm
K=10;


//----------------------------------
// MACROS
//----------------------------------


// Returns the minimum and the approximate of the function
function [minimum,V]=kelley(f,x_min,x_max,K)
 
  // V contains the coordinates of the cutting lines
  V=zeros(K,2);
  V(:,2)=-%inf;
    
  // Initialization of the algorithm
  [y,z]=oracle(x_min);
  V(1,1)=z;
  V(1,2)=y-z*x_min;
  V(1,3)=x_min;
  //Matrix interpretation of the problem
  C=[1,0];
  
  A=zeros(K+2,2);
  A(1:3,:)=[0,-1;0,1;-1,V(1,1)];
  
  b=zeros(K+2,1);
  b(1:3,1)=[-x_min;x_max;-V(1,2)];
  
  
  // Main loop, search for a new cutting line
  for k=2:K
    
    //Call for an external linear solver
    [u_opt,c_opt]=linprog(C,A,b,[],[],lb=[-%inf;-%inf]);
    
    // tangent on f at the minimum of the approximate
    [y,z]=oracle(u_opt(2,1));
    V(k,1)=z;
    V(k,2)=y-z*u_opt(2,1);
    // abscissa of the minimum of the approximate
    V(k,3)=u_opt(2,1);
    //New constraints linked to the new cutting line
    A(k+2,:)=[-1,V(k,1)];
    b(k+2,1)=-V(k,2);
  end
  minimum=[u_opt(2,1),oracle(u_opt(2,1))];
  
endfunction


// Computation and graphs
[m,V]=kelley(oracle,x_min,x_max,K);
delta=(x_max-x_min)/100;
X=x_min:delta:x_max;

Y=zeros(K,101)
for i=1:101
  Y(:,i)=V(:,1)*X(i)+V(:,2);
end

xset("window",11);xbasc();
realtimeinit(1);
plot2d(X,[oracle(X)]',rect=[x_min,-25,x_max,150],style=1);
xtitle('The oracle function and the cutting lines given by the Kelley"+...
       " algorithm')
xpause(2,%t);

for i=1:size(V,'r') 
  if %f then realtime(i);end 
    
  plot2d(X,[Y(i,:)]',rect=[x_min,-25,x_max,150],style=2);
  x =V(i,3); y = oracle(x);
  plot2d(x,y,rect=[x_min,-25,x_max,150],style=-3);
  if %f then xset("wshow");
  else
    xclick(clearq=%t,getkey=%t)
  end
end