function [RVecDec, Errors, GammaInv] = getDecVecC(NeuronParms,TypeParms,CntlParms)
%% Compute the decoding vectors for genNeuronVecRep() neurons
%% Sept. 12 , 2000
%% Copyright (C) by Charles. H. Anderson (All Rights Reserved)
%% Dept. Anatomy and Neurobiology
%% Washington Univ. School of Medicine
%% St. Louis, MO
%% cha@shifter.wustl.edu

%% NeuronParms are generated by genNeuronVecRep.m
%% Type 1 Neurons: Modified Rectified Linear
%%  NeuronParms=[Rthres,Rbreak,Gain1,Gain2,Bias,EUV];
%% Type 2 Neurons: LIF
%%  NeuronsParms=[Rthres,epsilon,Gain,maxFR,EUV];

%% [N,D+3] = size(NeuronParms); N Neuron number, D dimension of space.

%% CntlParms = Radius,dr,Neuron Noise Level 
%% Default values are
%% Radius = 1.0;
%% dr     = 0.05*Radius; (Computation time scales as (1/dr)^2;
%% noise = 0.1; (of mean maximum firing rate)   
   
%% NeuronTypeParms;
%% First one tells the neuron type (1==rectified linear or 2=LIF)  

%% The computation is
%% gamma(n,m) = \int ... \int a_n(R) a_m(R) dR, integration over hypersphere
%% 
%% a_n(R) = G[EUV_n \dot R], where G[] is the neuron nonlinear response.
%% Choose x1 axis along U_n direction
%% and    x2 axis along U_m(1-U_m \dot U_n)
%% The integrand does not depend on the other axes so
%% gamma(n,m) = \int\int a_n(x1)a_m(x1,x2) V(1-sqrt(x1^2+x2^2)) dx1 dx2 
%% Where the volume is that of a hypersphere of dimension (D-2).

%% Vol unit D dimensional hypersphere V(D) = pi^(D/2)/gamma((D/2)+1);

%% To compute the decoding vectors we need V(n) = \int...\int R a_n(R) dR
%  V(n) = U(n) \int V2(1-sqrt(x^2) G[gain(n) x+bias(n)] dr
%  where V2() is the volume of a (D-1) hypersphere.

if(nargin>2)
   Radius  = CntlParms(1);
   dr      = CntlParms(2);
   noise   = CntlParms(3);
else
   Radius = 1.0;
   dr     = 0.05;
   noise  = 0.1;
end
Type = TypeParms(1);
if(Type==1)
   EUV     = NeuronParms(:,6:end); %% Encoding Unit Vectors \tilde \phi_i
elseif(Type==2);
   EUV     = NeuronParms(:,5:end);
end
[N,D] = size(EUV);
Rthres = NeuronParms(:,1);
if(D<1)
   error('NeuronParms vector space is less than 1');
end
Gamma = zeros(N,N);
Vn    = zeros(D,N);
noiseWeights = zeros(1,N); %Fraction for which a_n(R) > 0 (The noise is weighted by this factor) 
nW2 = zeros(1,N); % noiseWeights computed a second way. 
r = [-Radius:dr:Radius];
if(D==1)
   tic
   hsVolD = 2*Radius; % Vol of 1D sphere.
   R2est = sum(r.^2)*dr/hsVolD;
   an = zeros(N,length(r));
   Rvalues = EUV*r;
   an = genActivitiesB(NeuronParms,Rvalues,TypeParms);
	posValues = (an>0);
	noiseWeights=sum(posValues')/length(r);
   Vn = an*r'*(dr/hsVolD);
   Gamma= an*an'*(dr/hsVolD) + noise^2*diag(noiseWeights); % 
   t2 = toc;
else
	[Y X] = ndgrid(r,r);
	X = X(:);
	Y = Y(:);
	R2 = X.^2+Y.^2;
	index = find(R2<=Radius^2);
	X = X(index);
	Y = Y(index);
	R2 = R2(index);

%% Compute unit hypersphere volumes
   hsVolD  = pi^(D/2)/gamma(D/2+1)*Radius^D;            % Vol of unit hypersphere D 
   hsVolDm1= pi^((D-1)/2)/gamma((D-1)/2+1)*Radius^(D-1);    % Vol of unit hypersphere D-1
   hsVolDm2= pi^((D-2)/2)/gamma((D-2)/2+1)*Radius^(D-2);    % Vol of unit hypersphere D-2   
	Vol2 = hsVolDm1.*(1-r.^2/Radius^2).^((D-1)/2.0); % Weighting factors for V_i calculation
	if(D>2)                    
   	Vol = hsVolDm2.*(1-R2/Radius^2).^((D-2)/2.0); % Weighting factors for <a_i*a_j> cal.
   end
   index = find(r>0);
   % Numerical estimate of <r^2> (See 
	R2est = D*dr*sum(r.^2.*Vol2)/hsVolD;  % \int_0^1 r^2 (dV/dr) dr
	Nx = length(X);
	tic
	Xt    = zeros(Nx,1);
   Yt    = zeros(Nx,1);
   Rt    = zeros(Nx,D-1);
	am    = zeros(length(X),1);
	an    = zeros(length(X),1);
	for(n=1:N)
		index2 = find(r>=Rthres(n));
      an2 = genActivitiesB(NeuronParms(n,:),r(index2),TypeParms);
      noiseWeights(n) = sum(Vol2(index2)); % Fraction of time neuron is on
		VWeights(n) = sum(an2.*Vol2(index2).*r(index2)); %<xa_n(x)>   
		index=find(X>=Rthres(n));
		Xt = X(index);
      Yt = Y(index);
      an = genActivitiesB(NeuronParms(n,:),Xt',TypeParms);
		if(D==2)
			nW2(n) = length(index);
			Gamma(n,n) = an*an';
			for(m=n+1:N)
				Ux = EUV(n,:)*EUV(m,:)';
            Uy = sqrt(1-Ux^2); 
         	Rvalues = (Ux*Xt+Uy*Yt);
            am = genActivitiesB(NeuronParms(m,:),Rvalues',TypeParms);
				index = find(am>0);
         	if(length(index)>0)
            	Gamma(n,m)= am(index)*an(index)';
					Gamma(m,n) = Gamma(n,m);
				end
      	end
		else
			Vt = Vol(index);
			nW2(n) = sum(Vt);
			anp = an.*Vt';	          % Weight by D-2 hypersphere
			Gamma(n,n) = an*anp';		
			for(m=n+1:N) 
				Ux = EUV(n,:)*EUV(m,:)';
         	Uy = sqrt(1-Ux^2);
         	Rvalues = Ux*Xt+Uy*Yt;
            am = genActivitiesB(NeuronParms(m,:),Rvalues',TypeParms); 
				index = find(am>0);
				if(length(index)>0)
					Gamma(n,m)= am(index)*anp(index)';
					Gamma(m,n) = Gamma(n,m);
				end
			end
		end
	end
   t2 = toc;
% Comment out the calculated noiseWeights and use 1.0 instead.
%	noiseWeights = noiseWeights*(dr/hsVolD);
	noiseWeights = ones(1,N);
  	Vn = (VWeights'*ones(1,D)).*EUV*(dr/hsVolD); %dr because V's are 1D integrals  
	noiseMax = max(noiseWeights);
	nW2 = nW2*(dr^2/hsVolD);
	n2Max = max(nW2);
	aaMax = max(max(Gamma))*dr^2/hsVolD;
	Vmax = max(max(Vn));
%	fprintf('noiseMax=%8.3e, n2Max=%8.3e, aamax=%8.3e, vMax=%8.3f\n',noiseMax,n2Max,aaMax,Vmax);
   Gamma = (dr^2/hsVolD)*Gamma + noise^2.*diag(noiseWeights); %dr^2 because Gamma is a 2D integral
%	Gamma = (dr^2/hsVolD)*Gamma+ noise^2.*diag(nW2);
end
tic
%GammaInv = inv(Gamma); % Much faster and gives the same answer as pinv();
GammaInv = pinv(Gamma,1.0e-10);
% keyboard
RVecDec = GammaInv'*Vn;
t3= toc;
fprintf('Gamma time=%8.3f, InvGamma=%8.3f\n',t2,t3)
%% Residual Error = <r^2>/hsVolD - phi gamma phi;
R2analytic = D*pi^(D/2)/((D+2)*gamma(D/2+1))*Radius^(D+2)/hsVolD; %Analytic <r^2> = D*Vol(D)/(D+1)
PhiGammaPhi = sum(diag(RVecDec'*Gamma*RVecDec)); 
Error = R2est - PhiGammaPhi;
Error2 = R2analytic - PhiGammaPhi; %% PhiGammaPhi should ~ <r^r>
RVecDecSqr = RVecDec.*RVecDec;
if(D>1)
	RVecDecMag = sum(RVecDecSqr');
else
	RVecDecMag = RVecDecSqr';
end
Enoise  = noise^2*noiseWeights*RVecDecMag'; 
Edist = Error-Enoise;   % Using an numerical estimate for <r^2>
Edist2 = Error2-Enoise; % Using the analytic value for <r^2>
Errors = [Error2 Enoise Edist2];
fprintf('R2 = %7.3f, R2est=%7.3f, PhiR2=%8.3f\n',R2analytic,R2est,PhiGammaPhi);
fprintf('Total Error=%8.3e, Enoise=%8.3e, Edist=%8.3e, Edist2=%8.3e\n',Error,Enoise,Edist,Edist2);