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.
%% 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 (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, [weightFlag]
%% Default values are
%% Radius = 1.0;
%% dr = 0.05*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.
if(nargin>2)
Radius = CntlParms(1);
dr = CntlParms(2);
else
Radius = 1.0;
dr = 0.05;
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
