genEnsemble_main(p)

function varargout = genEnsemble_main(p)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Script for generating neuronal represenations of vector spaces
%%
%% Copyright (C) by Charles. H. Anderson and Chris Eliasmith (All Rights Reserved)
%% Dept. Anatomy and Neurobiology
%% Washington Univ. School of Medicine
%% St. Louis, MO
%% [e-mail]
%% [e-mail])

%% Modified Dec. 5, 2003 CHA

%% CHA Oct 13,2003 all reference to tmpfiles are removed. 
%% Test for output directory added.
%% Modified by CHA, Jan. 28, 2003 to make it compatible with the standalone
%% gen_ensemble gui interface.
%%  a) The tmp.mat file is now saved in the defined database directory, usually NeuronData 
%%  b) The SAVE DATA button now saves the data, but after creating a new dataset
%%  c) The precision dR is changed to dR/2 for the 2D case.

%% 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
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%  Initialize the parameters using the structure p


D = p.dimension;
N = p.numberNeurons;
Nname = p.population;

if length(p.outputDirectory)>=7;
    if (p.outputDirectory(1:6)=='<path>')
        DataBaseStr = lower(p.outputDirectory(7:end));
        mpath = lower(path);
        k = findstr(DataBaseStr,mpath);
        kend = k+length(DataBaseStr)-1;
        while((k>1)&(mpath(k)~=';'));
            k=k-1;
        end
        p.outputDirectory = mpath(k:kend);
    end
end
if(isdir(p.outputDirectory)==0)
    fprintf('Error: gen_ensemble output directory: "%s" does not exist\n',p.outputDirectory );
    beep;
    return;
end
Database = [p.outputDirectory filesep];
fileName = [Database,Nname,'_N',num2str(N),'D',num2str(D)];

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.

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. 

if (RandSeed==-1)
    RandSeed=abs(sum(rand('state'))*100-1500); % turn the 35 element vector into an integer
end

%%%%%%%% 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

%%% CHA Change 1/28/03
if(D==2)      
    dR = dR/2;  %% The 2D algorithm is not as precise as the other  cases.
end

CntlParms = [Radius,dR];
[Cmatrix Moments] = getDecVecParms(NeuronParms,ModelType,CntlParms,WeightFlag);


%%%%%%%% Compute the linear decoding weights
DecVec = getDecVec(EncVec,Cmatrix,Moments,Noise);

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

%%%%%%%% CHA clean up 10/13/03 

save(fileName,'D','N','epsilon','ModelType','tref','trc','noise','maxFR','Radius','dR','RandSeed'...
        ,'WeightFlag','SRange','TRange','NeuronParms','EncVec','DecVec','Cmatrix','Moments');

fprintf('Saved data in %s\n',fileName);