function [] = JacobianKronecker(sizes)
%JACOBIANKRONECKER computes the Jacobian of a Kronecker product model for
%some random parameters, taking as argument sizes = [a,b]. Here a is the dimension of
%the visible exponential family plus one and b is the dimension of the
%hidden exponential family plus one.

%AUTHOR: Guido Montufar, montufar@mis.mpg.de
%DATE: Nov 10, 2015
%PURPOSE: Compute the dimension of marginals of Kronecker products of
%exponential families
%REFERENCES:See paper ``Dimension of Marginals of Kronecker Product Models;
%Geometry of hidden-visible products of exponential families''

% if no arguments are provided choose some prespecified values
if ~exist('sizes','var')
a=4+1; % dimension of the visible model plus one
b=2+1; % dimension of the hidden model plus one

% else set [a,b] as provided
else
a=sizes(1);
b=sizes(2);
end
clf;

% some random parameters
Theta = randn(b,a);

% define the sufficient statistics matrix of the visible exponential family
A = [ones(1,2^(a-1));binlist(a-1)']; % (a-1)-bit independence model

% we compare two types of models
for model=[1 2]

% define the first type of sufficient statistics matrix of the hidden exponential family
if model==1
B = eye(b); B(1,:) = ones(1,b); % (b-1)-dimensional simplex

% define the second type of sufficient statistics matrix of the hidden exponential family
else
B = [ones(1,2^(b-1));binlist(b-1)']; % (b-1)-bit independence model
end

% compute the conditional distribution of y given x for the parameters Theta
p = exp(B'*Theta*A);
p = p./repmat(sum(p,1),size(B,2),1); % normalize

% compute the expectation of B for each p(.|x)
EB = B*p;

% compute the Jacobian matrix
for x=1:size(A,2)
J(:,x) = kron(A(:,x),EB(:,x));
end

% compute the rank of the Jacobian
R = rank(J);

% visualize Jacobian matrix and its rank
if 1==1
figure(1)
subplot(1,2,model)
imagesc(J)
if model==1
title(sprintf(strcat('Mixture Model\n','rank(J)=',num2str(R),', expected=',num2str(min(size(A,2), a*b)),', parameters=',num2str(a*b) )'))
elseif model==2
title(sprintf(strcat('Hadamard Model\n','rank(J)=',num2str(R),', expected=',num2str(min(size(A,2), a*b)),', parameters=',num2str(a*b) )'))
end
end

% visualize expectation parameters of the conditional distributions
if b==3 % hidden model of dimension two
figure(2)
subplot(1,2,model)
scatter(EB(2,:),EB(3,:))
hold on
if a>=2
for i=[0:4:(2^(a-1)-4)]
ed = i+[1 2 4 3 1]; plot(EB(2,ed),EB(3,ed),'r')
end
end
axis equal
xlim([0, 1]);ylim([0, 1.1])
axis off
hold on
for x=1:2^(a-1)
text(EB(2,x),EB(3,x)+.02,num2str(x))
end
if model==1
plot(B(2,[1 2 3 1]),B(3,[1 2 3 1]),'k')
title(sprintf(strcat('Mixture Model\n','rank(J)=',num2str(R),', expected=',num2str(min(size(A,2), a*b)),', parameters=',num2str(a*b) )'))
elseif model==2
plot(B(2,[1 2 4 3 1]),B(3,[1 2 4 3 1]),'k')
title(sprintf(strcat('Hadamard Model\n','rank(J)=',num2str(R),', expected=',num2str(min(size(A,2), a*b)),', parameters=',num2str(a*b) )'))
end

% visualize e-geodesics
if 1==2
figure(2)
subplot(1,2,model)
clear e
for ix=1:2:size(A,2)
for iy=1:size(B,2)
e(iy,:) = p(iy,ix) * (p(iy,ix)/p(iy,ix+1)).^[-100:.1:100];
end
e = e./repmat(sum(e,1),size(B,2),1);
EB = B * e;
plot(EB(2,:),EB(3,:))
hold on
end
axis equal
xlim([0, 1]);ylim([0, 1])
end

hold off
end

end
end

function [v] = binlist(V,q)
%BINLIST gives a full list of binary vectors of length V
% binlist(V) outputs a matrix with rows equal to all binary vectors of
% length V.
% binlist(V,q) outputs a matrix with rows equal to all q-ary vectors of length V in
% entrywise one-hot representation.
if ~exist('q','var')
v = zeros(2^V,V);
str= dec2bin(0:2^V-1);
for j=1:V;
v(:,j) = str2num(str(:,j));
end
else
v=kron(ones(q^V,V) , 1:q)== kron(combn(1:q , V),ones(1,q));
end
end