function [simdata,pfailure] = Simulate(model,ndata);
% [simdata,pfailure] = Simulate(model,ndata)
% Generates (ndata) Monte Carlo random numbers drawn from a model with the stated parameters
% Returns this along with the proportion of failures.

APsucceed = ( rand(1,ndata) < model.pstim );
noise = randn(1,ndata) * model.noise;
siterelease = ( rand(model.n,ndata) < model.p ) ;
if (model.QVtype==1)
    qsize = ones(model.n,ndata)*model.q;
    qvariance = randn(model.n,ndata) * model.qsigma;
    simdata = noise + APsucceed .* sum ( siterelease .* ( qsize + qvariance ) , 1 ) + model.offset;
elseif model.QVtype==2 % "flat"
    qvariance = randn(1,ndata) * model.qsigma;
    simdata = noise + APsucceed .* sum ( siterelease .* model.q  , 1 ) + APsucceed .* (sum( siterelease,1)>0) .* qvariance + model.offset;
else
    errordlg('error in SImulate');
end
pfailure = 1-sum(APsucceed.*(sum(siterelease,1)~=0))/ndata;

% May 19th 2008: Finally, allow for a glitch rate, if present:
if isfield(model,'glitchrate') && model.glitchrate>0
    % Say: each EPSP has a [model.glitchrate] chance of being replaced with a
    % new value picked randomly from the interval [u1,u2]
    jreplace = find(rand(size(simdata))<model.glitchrate);
    simdata(jreplace) = model.glitchmin + (model.glitchmax-model.glitchmin)*rand(size(jreplace)) + model.offset;
end


