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