% rndseq - Sequence of values sampled from a random variable model
%
% rndseq ::
% 	  model xdom:[[1,XD]]   ~ random variable model, data size is xdom,
% 	  sdom:[[1,SD]],        ~ size of sample,
% -> seq([[xdom sdom]]) ~ size of rndseq is product of xdom and sdom
%
% rndseq ::
% 	  model xdom:[[1,XD]]   ~ random variable model, data size is xdom,
% 	  sdom:[[1,SD]],        ~ size of sample,
% 	  rndstate              ~ initial state of generators
% -> rndseq([[xdom sdom]]) ~ size of rndseq is product of xdom and sdom
%
% If an initial rndstate is supplied, rndseq is purely functional
% and referentially transparent. Otherwise, the initial construction
% uses the current state of the random generators. After this, the
% entire sequence is fully determined.
% 
% EXAMPLE
%
%  rndseq(gaussian,[2,200]) :: seq [[2,200]]
% 	rndseq(dirichlet(3,0.5),6) :: seq [[3,6]]
function d=rndseq(model,sdom,k)
	if nargin<2, sdom=1; end
	if isa(model,'struct')
		if model.nparams>0,
			error('Model has unbound parameters');
		end
		gen=sampler(model.sample,sdom);
	elseif iscell(model)
		gen=sampler(model{1},sdom);
	end

	if nargin>2, 
		if ~iscell(k), s={k,k}; else s=k; end 
	else
		s=getrndstate;
	end
	d=unfoldseq(gen,s);
end

% this creates a rand-state-managed version of g applied to args,
% ie,
% sampler :: 
%    ((A1,A2,....) -> random B) ~'an action which generates a value'
% -> (rndstata -> B, rndstate)  ~'a deterministic random state transformer'.
function f=sampler(g,varargin)
	f=@ufn;
	function [x,s]=ufn(s)
		setrndstate(s); x=g(varargin{:});
		s=getrndstate;
	end
end

