function [res,upper,lower] = mykl(J, fn)
  if nargin < 2, fn = 'kl'; end
  J = J / sum(sum(J));
  % picture of J:
  %  A \ B    true    false
  %  true     1,1     1,2
  %  false    2,1     2,2
  pA = [ J(1,1)+J(1,2) ; J(2,1)+J(2,2) ]
  pB = [ J(1,1)+J(2,1) ; J(1,2)+J(2,2) ];
  pAgivenB = [ J(1,1)/pB(1)     J(1,2)/pB(2) ;
               J(2,1)/pB(1)     J(2,2)/pB(2) ];
  pBgivenA = [ J(1,1)/pA(1)     J(1,2)/pA(1) ;
               J(2,1)/pA(2)     J(2,2)/pA(2) ]
  
  numRep = 100;
  res = zeros(numRep*numRep*numRep,7);
  upper = []; lower = [];
  
  div = numRep+1;
  %for rep=1:numRep,
  rep = 0;
  for rep1=1:numRep,
  if mod(rep1,10) == 0, fprintf(1, 'rep = %d/%d\n', rep1, numRep); end;
  for rep2=1:numRep,
  for rep3=1:numRep,
    rep = rep + 1;

    qA = [ rep1/div  ;  1-rep1/div ];
    qBgivenA = [ rep2/div      1-rep2/div ;
                 rep3/div      1-rep3/div ];
    
    %qA = noisify(pA,noiseRate); qA = qA / sum(qA);
    %qBgivenA = noisify(pBgivenA, 0.001);
    %qBgivenA = qBgivenA ./ [sum(qBgivenA,2) sum(qBgivenA,2)];
    
    qAB = [ qA(1)*qBgivenA(1,1)     qA(1)*qBgivenA(1,2) ;
            qA(2)*qBgivenA(2,1)     qA(2)*qBgivenA(2,2) ];
    qAgivenB = qAB ./ [sum(qAB) ; sum(qAB)];

    klBgivenA = ...
        pA(1) * dist(fn, qBgivenA(1,:), pBgivenA(1,:)) + ...
        pA(2) * dist(fn, qBgivenA(2,:), pBgivenA(2,:));
    
    klAgivenB = ...
        pB(1) * dist(fn, qAgivenB(:,1), pAgivenB(:,1)) + ...
        pB(2) * dist(fn, qAgivenB(:,2), pAgivenB(:,2));
    
    res(rep,:) = [ dist(fn, qA, pA)   klBgivenA   dist(fn, J,qAB)   klAgivenB   rep1/div   rep2/div   rep3/div ];

    if res(rep,1) > 0.14,
      x = res(rep,1);
      y = res(rep,3);
      if y > 0.7077 * x + 0.0338,
        upper = [upper qA(1)];
      else
        lower = [lower qA(1)];
        end
    end
    
    %if abs(res(rep,1) - 0.4) < 0.01 && abs(res(rep,3) < 0.249) < 0.01,
      %disp 'lower'
      %qA
      %qBgivenA
      %lower = [lower qA(1)];
    %end
    %if abs(res(rep,1) - 0.4) < 0.001 && abs(res(rep,3) < 0.369) < 0.001,
      %disp 'upper'
      %qA
      %qBgivenA
      %upper = [upper qA(1)];
    %end

      
    
%    fprintf(1, 'kl(qA||pA) = %g\n', 0);
%    fprintf(1, 'kl(qBgivenA||pBgivenA;pA) = %g\n', 0);
%    fprintf(1, 'kl(q||p) = %g\n', 0);
%    fprintf(1, 'kl(qAgivenB||pAgivenB;pB) = %g\n', 0);
  end ; end ; end

  minX = min(res(:,1)); maxX = max(res(:,1));
  minY = min(res(:,2)); maxY = max(res(:,2));
  
  spacing = 20;
  xrange = linspace(minX, maxX, spacing);
  yrange = linspace(minY, maxY, spacing);
  vals3= zeros(spacing,spacing);
  vals4= zeros(spacing,spacing);
  xs   = zeros(spacing,spacing);
  ys   = zeros(spacing,spacing);
  xsigma = (maxX - minX) / spacing / 10;
  ysigma = (maxY - minY) / spacing / 10;
  col = 4;
  for ii=1:spacing,
    for jj=1:spacing,
      sumWeight = 0.;
      for n=1:(numRep*numRep*numRep),
        dist = (1 / xsigma) * (res(n,1) - xrange(ii)) ^ 2 + (1 / ysigma) * (res(n,2) - yrange(jj)) ^ 2;
        weight = exp(- dist);
        sumWeight = sumWeight + weight;
        vals3(ii,jj) = vals3(ii,jj) + weight * res(n,3);
        vals4(ii,jj) = vals4(ii,jj) + weight * res(n,4);
      end
      vals3(ii,jj) = vals3(ii,jj) / sumWeight;
      vals4(ii,jj) = vals4(ii,jj) / sumWeight;
      xs(ii,jj) = xrange(ii);
      ys(ii,jj) = yrange(jj);
    end
  end
  figure;
  subplot(321)
  mesh(xs,ys,vals3);
  xlabel([fn  ' A']);
  ylabel([fn  ' B|A']);
  title([fn  ' joint']);

  subplot(323)
  plot(res(:,1), res(:,3), 'b.');
  xlabel([fn  ' A']);
  ylabel([fn  ' joint']);

  subplot(325)
  plot(res(:,2), res(:,3), 'b.');
  xlabel([fn  ' B|A']);
  ylabel([fn  ' joint']);

  subplot(322)
  mesh(xs,ys,vals4);
  xlabel([fn  ' A']);
  ylabel([fn  ' B|A']);
  title([fn  ' A|B']);

  subplot(324)
  plot(res(:,1), res(:,4), 'b.');
  xlabel([fn  ' A']);
  ylabel([fn  ' A|B']);

  subplot(326)
  plot(res(:,2), res(:,4), 'b.');
  xlabel([fn  ' B|A']);
  ylabel([fn  ' A|B']);
  
  suptitle(sprintf('p(a,b) = [%g %g ; %g %g]', J));
  
function z = noisify(p, noiseRate)
  z = exp(log(p) + noiseRate * randn(size(p)));
  
function kl = kl(q,p)  % kl of q FROM p == kl(p||q)
  if abs(sum(sum(q)) - 1) > 1e-5,
    fail
  elseif abs(sum(sum(p)) - 1) > 1e-5,
    fail
  end
  kl = 0.;
  for ii=1:size(p,1),
    for jj=1:size(p,2),
      kl = kl + p(ii,jj) * log(p(ii,jj) / q(ii,jj));
    end
  end
  
  
function d = l1(q,p)
  d = sum(sum(abs(q-p)));

function d = variational(q,p)
  d = max(max(abs(q-p)));

function d = dist(fn, q, p)
  if strcmp(fn, 'kl'),
      d = kl(q,p);
  elseif strcmp(fn, 'l1'),
      d = l1(q,p);
  elseif strcmp(fn, 'var'),
      d = variational(q,p);
  else
    fail
  end
  