function [Cmatrix, Moments] = getDecVecParms(NeuronParms,Type,CntlParms,WeightFlag)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Compute Cmatrix and Moments for genNeuronVecRep() neurons
%%
%% Cmatrix(n,m) = \int ... \int a_n(R) a_m(R) Weight(R) dR, integration over hypersphere
%%
%% The Moments are the first 5 moments of <a_n(R)R^l> if D=1,
%% otherwise they are the 0th, 1st and 2 second moments as defined below.
%%
%% NeuronParms are the parameters that define your neuron responses.
%% You will also have to define the function genActivities (which you've
%% already written) that accepts those parameters and returns the neuron
%% responses. genActivities needs to accept NeuronParms, the input values
%% to be encoded (i.e. x), and a Type flag that determines the kind of
%% neuron response (i.e. rectified or LIF).  
%%
%% Suggested structure for NeuronParms: [Rthres,TRC,Tref,Gain,Jbias,EncVec]
%% Whatever structure you choose, make sure it is used properly below
%%
%% [N,D+3] = size(NeuronParms); N Neuron number, D dimension of space.
%%
%% CntlParms = Radius,dr,Neuron Noise Level, [weightFlag] 
%% Default values are
%% Radius = 1.0;
%% dr     = 0.02*Radius; (Computation time scales as (1/dr)^2;
%% No noise is used in these calculations !!!
%%  
%% The computation is
%% Cmatrix(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
%% Cmatrix(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.
%%
%% Jan. 2, 2001 Inclusion of weighting function sqrt(1-R^2) 
%%          This is accomplished by changing the D to D+1 hyperphere factors
%%          \int dx1--dxD = V^{D} = 2 \int sqrt(1-\sum_{i=1}^{D-1} x_i^2) dx1--dx(D-1) 
%% So \int sqrt(1-R^2) dx1...dxD = V^{D+1}/2 !!! 
%% The factor of 2 is lost in computing <f(R)sqrt(1-R^2)> = \int f(R) sqrt(1-R^2)dx.../\int sqrt(1-R^2)dx..
%%
%% Jan. 8, 2001, Moment Computations 
%% M0 = <a_i(\X)>                        ; 0th moment
%% M1 = <\X \cdot \phi_i a_i(\X)>        ; 1st moment
%% M2 = <(\X \cdot \phi_i)^2 a_i(\X)>    ; 2nd moment 1st one
%% M3 = <\X^2 a_i(\X)>                   ; 2nd moment 2nd one 
%% M4 = <(\X \cdot \phi_i)^3 a_i(\X)>    ; 3rd moment 1st one
%% M5 = <(\X \cdot \phi_i) \X^2 a_i(\X)> ; 3rd moment 2nd one
%%
%% Sept. 9, 2001, Comments on weighting.
%% The averages are weighted by the factor [Radius^2-r^2]^WeightFlag.
%% The default is WeightFlag=0, which produces uniform weighting at all radii.
%% A positive WeightFlag emphasizes the errors near the origin.
%% The program will generate an error if the weightflag is < 0.
%%
%% Oct. 29, 2002 Released with NESim
%% Feb. 25, 2001 Consolidation of many older routines.
%% Modified Jan. 28, 2001 (Changed Gamma to Cmatrix to avoid Matlab's `Gamma' BUG)
%% Jan. 8, 2001
%% Derived from Jan. 2, 2001  version of getDecVecW.m 
%%
%% Copyright 2002 C. H. Anderson (cha@shifter.wustl.edu) and
%% C. Eliasmith (eliasmith@uwaterloo.ca)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(nargin>2)
   Radius  = CntlParms(1);
   dr      = CntlParms(2);
else
   Radius = 1.0;
   dr     = 0.02;
end
if(nargin<4)
    WeightFlag=0; % Default version is not to use the weighting factor.
else
    if(WeightFlag<0)
        fprintf('WeightFlag=%5.2f\n',WeightFlag);
        error('getDecVecParms: WeightFlag must be positive');
    end
end

EUV = NeuronParms(:,6:end); %% Encoding Unit Vectors \tilde \phi_i
[N,D] = size(EUV);
Rthres = NeuronParms(:,1);
if(D<1)
   error('NeuronParms vector space is less than 1');
end
Cmatrix    = zeros(N,N);
Moments    = zeros(N,6);  
r = [-Radius:dr:Radius];
r2 = r.^2;
r3 = r.^3;
if(D==1)
    r4 = r.^4;
    r5 = r.^5;
    hsVolD = 2*Radius; % Vol of 1D sphere.
 %   an = zeros(N,length(r));
    Rvalues = EUV*r;
    an = genActivities(NeuronParms,Rvalues,Type);
    if(WeightFlag) 
        WeightFactor = sqrt(Radius^2-r2).^WeightFlag;
        hsVolD = dr*sum(WeightFactor);
        Moments(:,1) = an*WeightFactor'*(dr/hsVolD);
        Moments(:,2) = an*(WeightFactor.*r)'*(dr/hsVolD);
        Moments(:,3) = an*(WeightFactor.*r2)'*(dr/hsVolD);
        Moments(:,4) = an*(WeightFactor.*r3)'*(dr/hsVolD);
        Moments(:,5) = an*(WeightFactor.*r4)'*(dr/hsVolD);
        Moments(:,6) = an*(WeightFactor.*r5)'*(dr/hsVolD);
        Cmatrix= (an.*(ones(N,1)*WeightFactor))*an'*(dr/hsVolD); 
    else
        Moments(:,1) = transpose(sum(an')*(dr/hsVolD));
        Moments(:,2) = an*r'*(dr/hsVolD);
        Moments(:,3) = an*r2'*(dr/hsVolD);
        Moments(:,4) = an*r3'*(dr/hsVolD);
        Moments(:,5) = an*r4'*(dr/hsVolD);
        Moments(:,6) = an*r5'*(dr/hsVolD); 
        Cmatrix= an*an'*(dr/hsVolD);  
    end
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 Weighting factors for computations

                      % Modification on Jan. 28, 2001     
    d = D+WeightFlag; % Increase D to impose the weighting factor sqrt(Radius^2-r^2)^WeightFlag 
    KD  = pi^(d/2)/gamma(d/2+1);          % Vol of unit D   hypersphere 
    KD_1= pi^((d-1)/2)/gamma((d-1)/2+1);  % Vol of unit D-1 hypersphere
    KD_2= pi^((d-2)/2)/gamma((d-2)/2+1);  % Vol of unit D-2 hypersphere
    hsVolD = KD*Radius^d;
    MD =  sqrt(pi)*(gamma((d)/2)/gamma((d+1)/2)-gamma((d+2)/2)/gamma((d+3)/2));
    F0 = (dr/hsVolD)*KD_1*(Radius^2-r2).^((d-1)/2);    % Weights for <(\phi \cdot \X)^n a_i(\X)>_\X
    F1 = (dr/hsVolD)*KD_2*MD*(Radius^2-r2).^((d+1)/2); % Weights for <(\phi \cdot \X)^2 x_2^2 a_i(x1)>
   	Fg = (dr^2/hsVolD)*KD_2*(Radius^2-R2).^((d-2)/2);  % Weighting factors for <a_i*a_j> cal.
	Nx = length(X);
	Xt    = zeros(Nx,1);
    Yt    = zeros(Nx,1);
    if(D==2)
        Rvalues = ones(N,1)*r;
        an2 = genActivities(NeuronParms,Rvalues,Type);
        Moments(:,1) = an2*F0';       % <a_n(\X)>
		Moments(:,2) = an2*(F0.*r)';  % < x_1 a_n(x_1)>   
        Moments(:,3) = an2*(F0.*r2)'; % <(x_1)^2 a_n(x_1)>
        Moments(:,4) = an2*F1';       % <(x_2)^2 a_n(x_1)>
%            Moments(:,5) = an2*(F0.*r3)'; % < (x_1)^3 a_n(x_1)>
%            Moments(:,6) = an2*(F1.*r)';  % < x_1(x_2)^2 a_n(x_1)>
        R = [X,Y];
        Rvalues = EUV*R';
        an = genActivities(NeuronParms,Rvalues,Type);
        an2 = an.*(ones(N,1)*Fg');
        Cmatrix = an*an2';
    else 
        am    = zeros(length(X),1);
	    an    = zeros(length(X),1);
        F01 = F0.*r;
        F02 = F0.*r2;
	    for(n=1:N)
		    index2 = find(r>Rthres(n)); % Only compute at pointes where the neurons are on
            an2 = genActivities(NeuronParms(n,:),r(index2),Type);
            Moments(n,1) = sum(an2.*F0(index2));                      % <a_n(\X)>
		    Moments(n,2) = sum(an2.*F01(index2));  % < x_1 a_n(x_1)>   
            Moments(n,3) = sum(an2.*F02(index2)); % <(x_1)^2 a_n(x_1)>
            Moments(n,4) = sum(an2.*F1(index2));             % < (x_2)^2 a_n(x_1)>
%            Moments(n,5) = sum(an2.*F0(index2).*r3(index2)); % < (x_1)^3 a_n(x_1)>
%            Moments(n,6) = sum(an2.*F1(index2).*r(index2));  % < x_1(x_2)^2 a_n(x_1)>
		    index=find(X>Rthres(n));
		    Xt = X(index);
            Yt = Y(index);
            an = genActivities(NeuronParms(n,:),Xt',Type);
            anp = an.*Fg(index)';
	        Cmatrix(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 = genActivities(NeuronParms(m,:),Rvalues',Type); 
			    index = find(am>0);
			    if(length(index)>0)
				    Cmatrix(n,m)= am(index)*anp(index)';
				    Cmatrix(m,n) = Cmatrix(n,m);
			    end
		    end
	    end
    end
end
