% [p,dp,d2p,b,pi] = penOWga(s,pot,tau,z,mu,rho,eta,theta) -
%                         Opper/Winther using convexity in ga Inner Loop penalty
% 
% pen(s) = -1/2 min_{b,pi} g(b,pi,s,z,mu,rho), where
%
%    g(b,pi,s,z,mu,rho) = z/pi - pi*s^2 + 2*b*s + 2*lZ/eta, with
%                  mun  = (mu-eta*b*rho)/(1-eta*pi*rho),
%                  rhon =           rho /(1-eta*pi*rho), and
%                  lZn  = \int exp(-(mun-t)^2/(2*rhon)) pot(t) dt.
%
% Here, theta are additional parameters for the potential function pot
% which is invoked by the call feval(pot, tau.*r, theta{:}).
%
%   See also PENFUNCTIONS.M.
% (c) by Hannes Nickisch, MPI for Biological Cybernetics, 2010 November 08

function [p,dp,d2p,b,pi,nfail] = penOWga(s,pot,tau,z,mu,rho,eta,varargin)

q = max([numel(s),numel(z),numel(mu),numel(rho)]);
pimax = 1./(eta*rho);
ub    = [ Inf(1,q); 0.999*pimax(:)'.*ones(1,q)];         % upper bound on (b;pi)
lb    = [-Inf(1,q); 0.001*pimax(:)'.*ones(1,q)];         % lower bound on (b;pi)

% find good initial values
b  = zeros(q,1); pi = pimax.*ones(q,1)/2;  % initial value
for j=1:5
  al    = 1-eta*pi.*rho;                                      % compute cavities
  mun   = (mu-eta*b.*rho)./al;
  rhon  =             rho ./al; 
  lZ    = feval(pot,tau.*mun,'EP',varargin{:},(tau.*tau).*rhon);
  be = tau.*lZ(:,2); nu = -tau.*tau.*lZ(:,3);
  X = solve_cubic( -(rho.*nu+2+eta^2*z(:).*rho), 1+2*rho.*nu, -rho.*nu );
  X = abs(X(:,3));
  pi = (1-X)./(eta*rho);
  b  = (mu+be.*rho-X.*s(:))./(eta*rho);
end
p = -critOWga([b;pi], s,pot,tau,z,mu,rho,eta,varargin{:})/2;

% iterate over individual
nfail = 0;
for j=1:q
  if numel(s)  ==q, sj   = s(j);   else sj   = s;   end
  if numel(z)  ==q, zj   = z(j);   else zj   = z;   end
  if numel(mu) ==q, muj  = mu(j);  else muj  = mu;  end
  if numel(rho)==q, rhoj = rho(j); else rhoj = rho; end
  if numel(eta)==q, etaj = eta(j); else etaj = eta; end
  if numel(tau)==q, tauj = tau(j); else tauj = tau; end
  try
    [bpi,npj] = wrap_lbfgsb([b(j); pi(j)], lb(:,j), ub(:,j), @critOWga, 20, ...
                                    sj, pot, tauj,zj,muj,rhoj,etaj,varargin{:});
    b(j) = bpi(1); pi(j) = bpi(2); p(j) = -npj/2;                % split outputs
  catch
    nfail = nfail+1;
  end
end

if nfail, fprintf('.'), end

if nargout>1, dp  = pi.*s(:)-b; end      % dir. deriv. of pen(s,z,mu,rho) wrt. s
if nargout>2, d2p = zeros(q,1); end                   % no second derivative yet