plotEnsemble(p,DecVec,NeuronParms,EncVec,Noise)

function plotEnsemble(p,DecVec,NeuronParms,EncVec,Noise)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Script for plotting neuronal representations of vector spaces
%%
%% Modified by CDE, Oct. 27, 2002
%%          Renamed, integrated with Java GUI, released with NESim
%% 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.
%% Cleaned up version by CHA, Sept. 10, 2001
%% Modified by John Harwell to interface with gen_Ensemble.m
%% Feb. 25, 2001 Consolidation 
%% Much revised version Jan. 10, 2000 
%% Dec. 13, 2000
%%
%% Copyright 2002 C. H. Anderson ([e-mail]) and
%% C. Eliasmith ([e-mail])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

D=p.dimension;
Radius=p.radius;
ModelType=p.modelType;
N=p.numberNeurons;

r = [-Radius:0.05*Radius:Radius];

if(D==1)
    Rvalues = EncVec*r; 
    a = genActivities(NeuronParms,Rvalues,ModelType);
else
    a = genActivities(NeuronParms,ones(N,1)*r,ModelType);
end

% Population plots    
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
    if ~isempty(findobj('Name','Population'))
        close(1);
    end
end
%Linearity plots
if p.plotLinearity == 1
    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 = EncVec*r';
    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 
    
    fh2 = figure(2);clf;
	if(D==1)
        [AX,h1,h2]=plotyy(r,Rest,r,(r-Rest),'plot');
        set(get(AX(1),'YLabel'),'String','r_{est}');
        set(get(AX(2),'YLabel'),'String','Error');
        %set(h1,'Color','r');
        %set(h2,'Color','g');
        hold on;
        h3=plot(r,r,'r');
        plotText = ['Linearity Plot,  rms_{DistError}=', num2str(DistRMSError)];        
        title(plotText);
        legend([h3 h1 h2],'r','r_{est}','r-r_{est}',2);
        xlabel('r');
	elseif(D==2)
        quiver(x,y,DeltaR(:,1),DeltaR(:,2));
        plotText = ['2D Error Vectors,  rmsError=', num2str(DistRMSError)];
        title(plotText);
        keyboard;
	else
        plot(sqrt(sum((DeltaR').^2)));
        plotText = ['Absolute error at selected points,  rmsError=', num2str(DistRMSError)];
        title(plotText);
        ylabel('Error');
        xlabel('Sample Number');
	end
    set(fh2, 'Name', 'Linearity');
    
    %display the errors
    totalRMSError = sqrt(DistRMSError^2+NoiseRMSError^2);
    fprintf('Distortion RMSError = %7.4f, Noise RMSError =%7.4f, Total = %7.4f \n', DistRMSError, NoiseRMSError,totalRMSError);
else
    if ~isempty(findobj('Name','Linearity'))
        close(2);
    end
end

% Decoding vector plots
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
    if ~isempty(findobj('Name','Decoding Vectors'))
        close(3);
    end
end  % if p.plotDecodingVectors