% Opper/Winther using convexity in ga Inner Loop criterion
%
% f(b,pi,s,z) = z/pi - pi*s^2 + 2*b*s + 2*lZn/eta, where
%               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.
%
% (c) by Hannes Nickisch, MPI for Biological Cybernetics, 2010 November 08

function [f,dbpi] = critOWga(bpi, s,pot,tau,z,mu,rho,eta,varargin)

b = bpi(1:end/2); pi = bpi(end/2+1:end);                          % split inputs
s = s(:); tau = tau(:); z = z(:); mu = mu(:); rho = rho(:);

al    = 1-eta*pi.*rho;                                        % compute cavities
mu_n  = (mu-eta*b.*rho)./al;
rho_n =            rho ./al;
lZn   = ( log(tau.*tau)+log(rho_n) + mu_n.*(mu./rho-eta*b) )/2;
lZ    = feval(pot,tau.*mu_n,'EP',varargin{:},(tau.*tau).*rho_n);

f     = lZ(:,1) + lZn;
f     = z./pi - pi.*s.*s + 2*b.*s + 2*f/eta;                      % add the rest
if nargout>1
  dmu_n  =  tau.*lZ(:,2);
  drho_n = (dmu_n.*dmu_n + (tau.*tau).*lZ(:,3))/2;
  db     = -eta*(rho_n.*dmu_n + mu_n);
  dpi    =  eta* rho_n.*(mu_n.*(dmu_n+(mu./rho-eta*b)/2)+drho_n.*rho_n +.5);
  db     = 2*s + 2*db/eta;                                        % add the rest
  dpi    = -z./(pi.*pi) - s.*s + 2*dpi/eta;
  dbpi   = [db;dpi];                                              % join outputs
end

% Below, there is the unrolled code for computing lZ, dlZ/db and dlZ/dpi
%   dmun_db   = -eta*rho./(1-eta*pi.*rho);
%   dmun_dpi  =  eta*rho.*(mu-eta*b.*rho)./(1-eta*pi.*rho).^2;
%   drhon_dpi =  eta*rho.^2./(1-eta*pi.*rho).^2;
%   lZ = feval(pot,mun,'EP',varargin{:},rhon);
%   d2lZ_dmun2 = lZ(:,3); 
%   dlZ_dmun   = lZ(:,2);
%   dlZ_drhon  = (d2lZ_dmun2.^2+dlZ_dmun)/2;
%   lZ = lZ(:,1);
%   dlZ_db  = dlZ_dmun.*dmun_db;
%   dlZ_dpi = dlZ_dmun.*dmun_dpi + dlZ_drhon.*drhon_dpi;