function varargout = gui_genEnsemble(p)
%% Script for generating neuronal represenations of vector spaces
%% Dec. 13, 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
%% Much revised version Jan. 10, 2000
%% Feb. 25, 2001 Consolidation
%% Modified by John Harwell to interface with gen_Ensemble.m
%% Cleaned up version by CHA, Sept. 10, 2001
%% Modified by CHA, Sept. 26, 2002
%% GEN ENSEMBLE data is saved in tmp.mat
%% tmp.mat is changed to desired filename when SAVE is hit
%% If tmp.mat does not exist when hitting SAVE, then the data is generated
%% and then saved with the desired filename.
%% CLOSE removes tmp.mat.
%%%%%%%% Initialize the parameters using the structure p
Database = [p.outputDirectory filesep];
tmpfile = [Database,'tmp.mat'];
Nname = p.population;
D = p.dimension;
N = p.numberNeurons;
ModelType = p.modelType;
RandSeed = p.randomSeed;
Radius = p.radius;
dR = p.dRadius;
TRange = [p.threshRangeMin, p.threshRangeMax];
SRange = [p.satRangeMin, p.satRangeMax];
noise = p.noise;
WeightFlag = p.weight;
tref = p.tauRefractory;
trc = p.tauRC;
epsilon = trc / tref;
maxFR = 1.0 / tref;
onsOnly = p.onsOnly; %!*************** ons only ***********************
if (D > 1)
onsOnly = 0;
end
%%%%%%%% Print out the parameter values.
if p.saveData == 1
status = exist(tmpfile);
if(status)
if(ModelType==1)
outfile = [Database,Nname,'_N',num2str(N),'D',num2str(D),'L.mat'];
else
outfile = [Database,Nname,'_N',num2str(N),'D',num2str(D),'.mat'];
end
status2 = dos(['mv ',tmpfile,' ',outfile]);
if(status2==0)
fprintf('Saved file %s\n',outfile) ;
return;
else
fprintf('Could not save output file %s.\n',outfile);
fprintf('Please make sure the directory %s exists.\n',Database);
return;
end
end
end % if p.saveData
fprintf('************\n');
fprintf('Name %s: N=%d, D=%d\n',Nname,N,D)
fprintf('ModelType=%d, WeightFlag=%d, Noise=%5.2f, RandSeed=%d \n',ModelType,WeightFlag,noise,RandSeed);
fprintf('Radius = %5.2f, TRange=[%5.2f,%5.2f], SRange=[%5.2f,%5.2f]\n',Radius,TRange,SRange);
fprintf('tref = %5.2f(ms), trc = %5.2f(ms), maxFR %5.1f at r=R\n',1000/maxFR,epsilon*maxFR/1000,maxFR*SRange(2));
%%%%%%%% Print out the errors.
%%%%%%%% Create the neuron parameters %%%%%%%%%%%%%%
Noise = noise*SRange(2)*maxFR; %% Rescale the noise level
TypeParms=[ModelType,Radius,epsilon]; %% Used in genNeuronVecRep()
if onsOnly==1
NeuronParms = genNeuronVecRepOn(N,D,TypeParms, TRange, SRange, maxFR, RandSeed);
else
NeuronParms = genNeuronVecRep(N,D,TypeParms,TRange,SRange,maxFR,RandSeed);
end
%%%%%%% NeuronParms(:,1:5) = [Rthreshold,Gain,maxFR,tRC,Jbias] for each neuron.
EncVec = NeuronParms(:,6:end); %% Pulls out the encoding vectors from NeuronParms.
%%%%%%%% Compute Cmatrix = <an(R)am(R)> and various moments
CntlParms = [Radius,dR];
[Cmatrix Moments] = getDecVecParms(NeuronParms,ModelType,CntlParms,WeightFlag);
%%%%%%%% Compute the linear decoding weights
DecVec = getDecVec(EncVec,Cmatrix,Moments,Noise);
%%%%%%%% Set up the activities for plots along the preferred directions.
r = [-Radius:0.05*Radius:Radius];
Nr = length(r);
if(D==1)
Rvalues = EncVec*r;
a = genActivities(NeuronParms,Rvalues,ModelType);
else
a = genActivities(NeuronParms,ones(N,1)*r,ModelType);
end
if p.plotPopulation == 1
fh1 = figure(1);clf;
plot(r,a);
title('Neuron firing rates along prefered direction');
xlabel('R');
ylabel('Firing Rate');
set(fh1, 'Name', 'Population');
else
figure(1);clf;
end
%%%%%%%% Compute the linearity, or precision of the representation.
%%% This code shows how to
%%% 1) Set up sample vectors in the various space.
%%% 2) Encode these vectors.
%%% 3) The compute the decoded values.
if(D==1)
r = r'; %% Uniform samples along the axis.
Nsamples = length(r);
elseif(D==2)
[X,Y] = ndgrid(r,r);
R2 = X.^2+Y.^2;
index = find(R2<=Radius^2); %% Uniform sampling over the circle.
x=X(index);
y=Y(index);
r=[x,y];
Nsamples = length(index);
else
Nsamples = 1000;
r = zeros(Nsamples,D);
cnt=0; %% Random samples in the D hypersphere.
while(cnt ~= Nsamples)
r0 = 2*Radius*rand(Nsamples-cnt,D)-Radius*ones(Nsamples-cnt,D);
r2 = sum((r0').^2);
index = find(r2<Radius.^2);
r(cnt+1:cnt+length(index),:) = r0(index,:);
cnt = cnt+length(index);
end
end
%%%%%%%% RValues are the projection of the vector r along the prefered directions.
RValues = EncVec*r';
%%%%%%%% Generate the activites using the RValues.
a = genActivities(NeuronParms,RValues,ModelType);
%%%%%%%% Decode the encoded vectors
Rest = a'*DecVec;
%%%%%%%% Compute the errors
DeltaR = r-Rest;
MSE = sum(sum((DeltaR').^2))/Nsamples;
DistRMSError = sqrt(MSE);
NoiseRMSError = Noise*sqrt(sum(sum(DecVec.^2))); % noise^2*sum_n |DecVec(n)|^2
%%%%%%%% Generate the plots
if p.plotLinearity == 1
fh2 = figure(2);clf;
if(D==1)
plot(r,Rest);
hold on;
plot(r,r,'r');
ErrScale = 0.25*Radius/DistRMSError;
plot(r,ErrScale*(r-Rest),'g');
plotText = ['Linearity Plot, rms_DistError=', num2str(DistRMSError)];
title(plotText);
legendText = [num2str(ErrScale,'%6.1f'),'*(r-r_{est})'];
legend('r','r_{est}',legendText,2);
xlabel('r');
ylabel('r_est');
elseif(D==2)
quiver(x,y,DeltaR(:,1),DeltaR(:,2));
plotText = ['2D Error Vectors, rmsError=', num2str(DistRMSError)];
title(plotText);
else
plot(sqrt(sum((DeltaR').^2))/Radius);
plotText = ['Absolute error at selected points, rmsError=', num2str(DistRMSError)];
title(plotText);
ylabel('RMS error');
xlabel('Sample Number');
end
set(fh2, 'Name', 'Linearity');
else
figure(2);clf;
end
if p.plotDecodingVectors == 1
fh3 = figure(3);clf;
if(D==1)
plot(sqrt(DecVec.^2))
else
plot(sum(DecVec'.*DecVec'));
end
xlabel('Neuron Number');
title('DecVec Magnitude');
set(fh3, 'Name', 'Decoding Vectors');
else
figure(3);clf;
end % if p.plotDecodingVectors
%%%%%%%% Save the data.
if p.saveData == 1
if(ModelType==1)
outfile = [Database,Nname,'_N',num2str(N),'D',num2str(D),'L.mat'];
else
outfile = [Database,Nname,'_N',num2str(N),'D',num2str(D),'.mat'];
end
else
outfile = tmpfile;
end % if p.saveData
save(outfile,'D','N','epsilon','ModelType','tref','trc','noise','maxFR','Radius','dR','RandSeed'...
,'WeightFlag','SRange','TRange','NeuronParms','EncVec','DecVec','Cmatrix','Moments');
%%%%%%%% Print out the errors.
totalRMSError = sqrt(DistRMSError^2+NoiseRMSError^2);
fprintf('Distortion RMSError = %7.4f, Noise RMSError =%7.4f, Total = %7.4f \n', DistRMSError, NoiseRMSError,totalRMSError);
if p.saveData == 1
fprintf('Saved file %s\n',outfile);
else
fprintf('DATA NOT SAVED\n');
end
