genNeuronVecRep(N,D,TypeParms,TRange,SRange,maxFR0,RandSeed)

function [Neurons] = genNeuronVecRep(N,D,TypeParms,TRange,SRange,maxFR0,RandSeed)
%% Generates a population of N neurons that represent a D dimensional vector space.

%% Type 1 Neurons: Rectified Linear
%%  Neurons=[Rthres,Gain,maxFR,Sat,zeros(N,1),EUV];
%%  The zeros(N,1) is utilize to maintain compatibility with the LIF Neuron data set.
%% Type 2 Neurons: LIF
%%  Neurons=[Rthres,Gain,maxFR,tRC,Jbias,EUV];
%% 
%%  EUV are random unit encoding vectors in the space. 

%% Rthres is the threshold value along the preferred direction
%% Sat sets the degree of saturation at R=Radius (TypeParms(2)). (Sat must be < 1.0) 
%% maxFR0 = 1/(refractory period) for Type 2
%% maxRF0 = firing rate at r=Radius for Type 1 neurons
%% 
%% Trange [Tmin,Tmax] sets the range for Rthres = [TRange(2)-TRange(1)]*rand()+TRange(1).
%% 
%% SRange[Smin, Smax]  sets the range for the degree of saturation at R=Radius
%%
%% If RandSeed is >0 it resets the random number generator seed.
%% 

%% Feb. 25, 2001 Consolidation

%% August 14, 2000
%% Modified Dec. 1, 2000
%% Modified Jan. 17, 2001 (General cleaning up of neuron subroutines on this date)
%% Copyright (C) by Charles. H. Anderson (All Rights Reserved)
%% Dept. Anatomy and Neurobiology
%% Washington Univ. School of Medicine
%% St. Louis, MO
%% [e-mail]

Type = TypeParms(1);
if((Type<1)|(Type>2))
   error('Type must be between 1 and 2\n');
end
Radius = TypeParms(2);
if(length(TRange)<2)
   error('TRange must be a 2D vector');
end
if(length(SRange)<2)
   error('SRange must be a 2D vector');
end
if(TRange(1)>TRange(2))
   error('TRange(1) must be less than TRange(2)');
end
EUV = genuniD(N,D,RandSeed);  % Use -1 for true random numbers
Rthres = (TRange(2)-TRange(1))*rand(N,1)+TRange(1);
if(D==1)
    iOn = find(EUV>0);
    Rthres(iOn) = sort(Rthres(iOn));
    Rthres(iOn) = flipud(Rthres(iOn));
    iOff = find(EUV<0);
    Rthres(iOff) = sort(Rthres(iOff));
    Rthres(iOff) = flipud(Rthres(iOff));
else    
    Rthres = sort(Rthres); % Sort them highest threshold
    Rthres = flipud(Rthres);    % to smallest threshold
end
%%%%%
Sat = (SRange(2)-SRange(1))*rand(N,1)+SRange(1);
%%%%%%%
if(Type==1)
   maxFR = Sat*maxFR0;
   Gain =  maxFR./(Radius-Rthres);
   Jbias = -Gain.*Rthres;
   Neurons=[Rthres,Gain,maxFR,Sat,Jbias,EUV]; % 
elseif(Type==2)
   epsilon=TypeParms(3);
   maxFR = maxFR0*ones(N,1); 
   tRC   = epsilon./maxFR; % Note that tRC and maxFR are constants!!!
   Gain= 1./((exp((1-Sat)./(Sat.*epsilon))-1).*(Radius-Rthres));   
   Jbias = 1.0-Gain.*Rthres;
   Neurons=[Rthres,Gain,maxFR,tRC,Jbias,EUV];
end