1. lft_llt


//Main function: Conjugate of an univariate function
function [Conj]=pl_lft_llt(X,Y,S,fusionopt)
	if (argn(2)>3 & fusionopt~=1) 
 		fusioncall=_pl_fusion;//scilab syntax for fast execution
 	else
 		fusioncall=_pl_fusionsci;//programming syntax for linear-time execution
 	end
 	if (Y == %inf.*Y | Y == %inf.*(-1.*Y))
 		Conj = %inf.*(-1.*Y);
 	else
 		[XX,YY]=pl_bb(X,Y); // Compute the convex hull (input sensitive algorithm).	
 		h=size(XX,1);
		C=(YY(2:h)-YY(1:h-1))./(XX(2:h)-XX(1:h-1)); // Compute the slopes associated with the primal points
 		[H]=fusioncall(C,S);// Merge sequences C and S
 		// j=H(i,1) is the first indice of the point X(j) such that the line
 		// with slope S(i) supports the epigraph of the function at X(j).
 		Conj= S .* XX(H(:,1)) -YY(H(:,1));
 	end
 endfunction
//OLD VERSION


2. opt
function y=_AttouchWetsKernel(x),
  global mu p
  y=0.5^(p-1)/(mu*p)*abs(x)^p;
endfunction

function y=_dAttouchWetsKernel(x),
  global mu p
  y=0.5^(p-1)/(mu)*sign(x)*abs(x)^(p-1);
endfunction


function [y,argmin] = opt_me(f1,df1,domf1,x,lambda,output)
  //Moreau envelope M_lambda(f)(x)=inf_y [f(y)+||x-y||^2/lambda]
  if argn(2)<6 then output=0;end;
  //cannot use dynamic scoping to pass parameters since the function is called by
  //an external (fortran) function. So use global variables
  global mu p
  mu = lambda;
  p = 2;
  f2=_zero;df2=_zero;domf2=[-%inf,%inf];
  //this gives f1 \infconv ||.||^p/(mu * p) so for p=2 get the Moreau envelope of f1
  [y,argmin] = opt_ka(f1,df1,domf1,f2,df2,domf2,_AttouchWetsKernel,_dAttouchWetsKernel,[-%inf,%inf],x,0.5,output);  
  clearglobal mu p
endfunction


function [_f,_g,ind] = _costf_ic(y,ind,x),
  //input:
  //ind=1: not output; ind=2: only compute f; ind=3: only compute g; ind=4: compute f and g
  //output: ind<0 cannot compute f at x; ind=0 interrupt optimization
  if ind==2 | ind==4 then 
    _f=f1(y)+f2(x-y);
  end
  if ind==3 | ind==4 then 
    _g=df1(y)-df2(x-y);
  end

  //workaround bug in optim with method 'nd': when gradient is exactly 0 errors out with optim: Wrong input parameters.
  //workaround is to replace 0 with a small value
  if _g==0 then _g=1E-12;end;
  //what do we return for a point outside the domain? The choice below while not the most
  //efficient works on our examples.
  if _f==%inf then _f=number_properties('huge');_g=%inf;end;

  //printf("y=%f, f=%f, g=%f\n",y,_f,_g);
endfunction




function x0 = _findFeasiblePoint(xtry)
  //heuristic search to find a feasible x0
  //look for x0 on an enlarged x grid
  //WARNING: dynamic scoping used to pass: x,lambda
  y=xtry;z = _costf_ka(y,2,x(1),lambda);
  if isinf(z) | isnan(z) then
    n=max(30,min(length(x),100));//n is between 30 and 100
    xmin=min(min(x),-10);xmax=max(max(x),10);
    Y=linspace(xmin,xmax,n);
    for y=Y
      z = _costf_ka(y,2,x(1),lambda);
      if ~isinf(z) then break;end;
    end
  end
  x0=y;
  if isinf(z) then x0=%nan;end;//search failed
endfunction
