function [NUB,OUB,PUB] = starcover(v)
% starcover evaluates upper bounds on the minimal size of star coverings of
% the k-interaction model for k=2,...,v. 
% author: Guido Montufar
% date: December 2015
% references: See the publication ``Hierarchical Models as Marginals of Hierarchical Models''. 
% Preprint available at www.arxiv.org/abs/1508.03606
%
% v is the number of visible binary variables
% calling starcover(v) prints upper bounds on the minimal size of star
% coverings for the k-interaction model with v binary variables for
% k=2,...,v. These are upper bounds on the minimal number of hidden binary
% variables for which an RBM contains the k-interaction model in its
% closure. 
%
% NUB is the new upper bound for universal approximation 
% OUB is the old upper bound for universal approximation
% PUB is the parameter counting lower bound for universal approximation

indbds=1;                    % if true the function uses improved individual bounds for D(v,j)=C(v,v-j+1,v-j) where available. 
                             % otherwise the function uses only the exact results and general bounds from Lemma 7

for k=2:v
U(v,k,indbds)                % print new upper bound for approximation of the k-interaction model
end
NUB =U(v,v,indbds);          % new universal approximation upper bound
OUB =2^(v-1)-1;              % old universal approximation upper bound
PUB =ceil((2^(v))/(v+1) -1); % parameter counting lower bound 
end

%% parameter counting bound 
function [up]=UP(v,k)
up=0;
for j=1:k
up = up + nchoosek(v,j);
end
up = ceil((up + 1)/(v+1) -1);
end

%% star covering of layers 2--k
function [u]=U(v,k,indbds)
u=0;
for j=2:k
u = u + D(v,j,indbds);
end
end

%% star covering of layer j
function [d]=D(v,j,indbds)
if indbds==1
% improved individual bounds taken from ``Nurmela and Ostergard, New Coverings of t-Sets With (t+1)-Sets'' 
Mindbds = [0 0 0   0    0     0     0      0     ;
           0 0 0   0    0     0     0      0     ;
           0 0 0   0    0     0     0      0     ;
           0 0 1   0    0     0     0      0     ;
           0 0 4   1    0     0     0      0     ;
           0 0 6   5    1     0     0      0     ;
           0 0 12  9    6     1     0      0     ;
           0 0 14  20   12    7     1      0     ;
           0 0 25  30   30    16    8      1     ;
           0 0 30  51   50    45    20     9     ;
           0 0 47  66   100   84    63     25    ;
           0 0 57  113  132   176   126    84    ;
           0 0 78  157  245   264   297    185   ;
           0 0 91  230  371   508   471    471   ;
           0 0 124 295  580   825   979    789   ;
           0 0 140 405  808   1329  1722   1768  ;
           0 0 183 492  1213  2048  3040   3355  ;
           0 0 207 664  1547  3261  4690   6098  ;
           0 0 258 846  2175  4608  7949   10641 ;
           0 0 285 1083 2900  6765  12134  18590 ;
           0 0 352 1251 3979  9338  18894  29961 ;
           0 0 385 1573 4687  13244 27624  48855 ;
           0 0 466 1771 6169  17375 38962  75163 ;
           0 0 510 2237 7084  23276 56337  113227;
           0 0 600 2706 9321  29759 78012  168996;
           0 0 650 3306 11952 39080 105000 239815;
           0 0 763 3848 15210 50895 144080 334854;
           0 0 819 4550 18369 65286 193595 470340];

end

if j==1                       % in this case the exact value is known
  d=1;
elseif j==2                   % in this case the exact value is known
  d=v-1;
elseif j==3                   % in this case the exact value is known
  d=1;
  for t=(v-3):-1:1 
  d=ceil( ((v-t+1)/((v-2)-t+1)) * d );
  end
elseif j==v-3 && mod(v,12)~=7 % in this case the exact value is knonw
  d=ceil((v/4)*ceil(((v-1)/3)*ceil((v-2)/2)));
elseif j==v-2                 % in this case the exact value is known
  d=ceil((v/3)*ceil((v-1)/2));
elseif j==v-1                 % in this case the exact value is known
  d=ceil(v/2);
elseif j==v                   % in this case the exact value is known
  d=1;
else
  d=min(floor(((1+log(v-j+1))/(v-j+1))*nchoosek(v,j)),...% general bound 
        nchoosek(v-1,j-1));                              % simple bound 
  if indbds==1 && v<=28 && (v-j)<=8 && Mindbds(v,v-j)~=0 % (v-j)=t
  d = Mindbds(v,v-j);                                    % individual bounds
  end
end

end

