function [lam,Gam,CI_eigvalue,CI_ratio,CI_eigvalue_log,CI_ratio_log] = ...
         pca_infer(X,Y,C,m,Normal_Dist,Cov_Mat,alpha)
%
% C is a d x g matrix.  Each column of C consists of coefficients
%   for a linear function of the eigenvalues.
%
% Set Normal_Dist to 1 to Assume Normality
% Set Normal_Dist to 0 to Not Assume Normality
% Default is Normal_Dist=1
%
% Set Cov_Mat to 1 to analyze the covariance matrix
% Set Cov_Mat to 0 to analyze the correlation matrix
% Default is Cov_Mat=1
%
% m = vector of eigenvalue multiplicities
% Default is m = ones(d,1)
%
[N,d]=size(Y);

if nargin ==3
  m=ones(d,1);
  Normal_Dist=1;
  Cov_Mat=1;
  alpha=0.05;
end

if nargin==4
  Normal_Dist=1;
  Cov_Mat=1;
  alpha=0.05;
    if isempty(m)
    m=ones(d,1);
  end
end

if nargin==5
  Cov_Mat=1;
  alpha=0.05;
  if isempty(m)
    m=ones(d,1);
  end
  if isempty(Normal_Dist)
    Normal_Dist=1;
  end
end

if nargin==6
  alpha=0.05;
  if isempty(m)
    m=ones(d,1);
  end
  if isempty(Normal_Dist)
    Normal_Dist=1;
  end
  if isempty(Cov_Mat)
    Cov_Mat=1;
  end
end

if nargin==7
  if isempty(m)
    m=ones(d,1);
  end
  if isempty(Normal_Dist)
    Normal_Dist=1;
  end
  if isempty(Cov_Mat)
    Cov_Mat=1;
  end
  if isempty(alpha)
    alpha=0.05;
  end
end
crit=norminv(1-alpha/2);
A=speye(N)-ppo(X);
r=rank(X);
n=N-r;
S=Y'*A*Y/n;
s=vec(S);
if Cov_Mat==1
  [Gam,Lam,Tmp]=svd(S);
else
  ds=diag(S);
  Dmhs=diag(1./sqrt(ds));
  R=Dmhs*S*Dmhs;
  [Gam,Lam,Tmp]=svd(R);
end
lam=diag(Lam);
Id=speye(d);
Ld=0;
for i=1:d
  ei=Id(:,i);
  Ld=Ld+kron(ei,ei)*ei';
end
Idd=commute(d,d);
Id2=speye(d^2);
twoNd=Id2+Idd;

k=length(m);

T=[];
for i=1:k;
  T=oplus(T,ones(m(i),1));
end;
Omega_n=twoNd*kron(S,S);
if Normal_Dist == 0
  U=A*Y;;
  c1=sum(diag(A).^2);
  c2=sum(sum(A.^4));
  c0=n*(n+2)*c2-3*c1^2;
  g1=n^2*c1/c0;
  g2=-n*(2*n*c2+(n-3)*c1^2)/((n-1)*c0);
  g3=-n^2*(c1^2-n*c2)/((n-1)*c0);
  Xi=0;
  for i=1:N
    ui=U(i,:)';
    ui2=kron(ui,ui);
    Xi=Xi+ui2*ui2';
  end
  Xi=Xi/n;
  Omega_n=g1*Xi+g2*s*s'+g3*Omega_n;
end


if Cov_Mat==0
  Wd=spdiag(vec(Id));
  L1=(eye(d^2)-.5*twoNd*kron(R,Id)*Wd)*kron(Dmhs,Dmhs);
  Omega_n=L1*Omega_n*L1';
end

D=T'*T;
Dmi=inv(D);
V_phi=Dmi*T'*Ld'*kron(Gam,Gam)'*Omega_n*kron(Gam,Gam)*Ld*T*Dmi;
V_lam=T*V_phi*T';
phi=Dmi*T'*lam;
lam=T*phi;

[d1 g]=size(C);
if abs(d1-d) > 0
  disp(['Error: Each column must have d components'])
end
CI_eigvalue=[];
CI_ratio=[];
CI_eigvalue_log=[];
CI_ratio_log=[];
h=ones(d,1);
for II=1:g
  c=C(:,II);
  psi=c'*lam;
  se_raw=sqrt(c'*V_lam*c)/sqrt(n);
  se_log=se_raw/abs(psi);
  CI_eigvalue=[CI_eigvalue; psi se_raw psi-se_raw*crit psi+se_raw*crit];
  CI_eigvalue_log=[CI_eigvalue_log
         psi se_log exp(log(psi)-se_log*crit) exp(log(psi)+se_log*crit)];

  Ph=h*inv(h'*lam)*lam';
  ImP=Id-Ph;
  psi=(c'*lam)/(h'*lam);
  se_raw=sqrt(c'*ImP'*V_lam*ImP*c/(n*(h'*lam)^2));
  se_log=se_raw/abs(psi);
  CI_ratio=[CI_ratio; psi se_raw psi-se_raw*crit psi+se_raw*crit];
  CI_ratio_log=[CI_ratio_log
     psi se_log exp(log(psi)-se_log*crit) exp(log(psi)+se_log*crit)];
end