function OUTP = mysourcedb(THIS,range,varargin)
% mysourcedb  [Not a public function] Create model-specific source database.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (C) 2007-2012 Jaromir Benes.

% Parse options.
opt = passvalopt('model.sourcedb',varargin{:});

ndraw = opt.ndraw;
ncol = opt.ncol;

if strcmp(opt.dtrends,'auto')
    opt.dtrends = ~opt.deviation;
end

if ~isnumeric(range)
    error('Incorrect type of input argument(s).');
end

%**************************************************************************

nalt = size(THIS.Assign,3);

if (ncol > 1 && nalt > 1) ...
        || (ndraw > 1 && nalt > 1)
    utils.error('model', ...
        ['The options `''ncol=''` or `''ndraw=''` can be used only with ', ...
        'single parameterisation models.']);
end

nf = sum(imag(THIS.solutionid{2}) > 2);
maxlag = -(min(imag(THIS.solutionid{2}(nf+1:end))) - 1);
nper = length(range);
xrange = range(1)-maxlag : range(end);
nxper = length(xrange);

n1 = sum(THIS.nametype == 1);
n2 = sum(THIS.nametype == 2);
n3 = sum(THIS.nametype == 3);

if ncol > 1 && nalt > 1
    utils.error('model', ...
        ['Input argument NCOL can be used only with ', ...
        'single-parameterisation models.']);
end

n = n1 + n2 + n3;
nloop = max([nalt,ncol,ndraw]);
OUTP = struct();

if opt.deviation
    X = zeros(n,nxper,nalt);
else
    tvec = double(round(xrange - xrange(1)));
    X = zeros(n,nxper,nalt);
    X(1:n1+n2,:,:) = mytrendarray(THIS,1:n1+n2,tvec,false,Inf);
end

if opt.dtrends
    D = mydtrendsrequest(THIS,'range',xrange);
    X(1:n1,:,:) = X(1:n1,:,:) + D;
end

X(THIS.log(1:n),:,:) = exp(X(THIS.log(1:n),:,:));

if nloop > 1
    X = X(:,:,ones(1,nloop));
end

tmp = tseries();
for i = find(THIS.nametype <= 3)
    OUTP.(THIS.name{i}) = replace(tmp,permute(X(i,:,:),[2,3,1]), ...
        xrange(1),THIS.namelabel{i});
end

% Generate random residuals if requested.
if opt.randomshocks
    OUTP = xxrandomshocks(THIS,OUTP,nper,maxlag,nloop,ndraw);
elseif ~isempty(opt.residuals)
    % SOURCEDB not implemented yet for models with cross-correlations.
    ne = sum(THIS.nametype == 3);
    if any(any(THIS.stdcorr(1,ne+1:end,:) ~= 0))
        utils.error('model', ...
            ['Source databases with ''residuals'' option ', ...
            'not implemented yet for models with cross-correlated shocks.']);
    end
    OUTP = xxaddrand(THIS,OUTP,opt.residuals,nxper,nloop);
end

% Add parameters.
for i = find(THIS.nametype == 4)
    OUTP.(THIS.name{i}) = permute(THIS.Assign(1,i,:),[1,3,2]);
end

end

% Subfunctions.

%**************************************************************************
function D = xxaddrand(THIS,D,FCN,NPER,NLOOP)
nalt = size(THIS.Assign,3);
ne = sum(THIS.nametype == 3);
elist = THIS.name(THIS.nametype == 3);
stdvec = THIS.stdcorr(1,1:ne,:);
for ie = 1 : ne
    x = zeros(1,NPER,NLOOP);
    for iloop = 1 : NLOOP
        if iloop <= nalt
            std = stdvec(1,ie,iloop);
        end
        x(1,:,iloop) = FCN(std,[1,NPER]);
    end
    D.(elist{ie}) = replace(D.(elist{ie}),permute(x,[2,3,1]));
end
end
% do_addrand().

%**************************************************************************
function D = xxrandomshocks(THIS,D,NPER,MAXLAG,NLOOP,NDRAW)
ne = sum(THIS.nametype == 3);
nalt = size(THIS.Assign,3);
elist = THIS.name(THIS.nametype == 3);
for iloop = 1 : NDRAW
    if iloop <= nalt
        Omg = covfun.stdcorr2cov(THIS.stdcorr(1,:,iloop),ne);
        F = covfun.factorise(Omg);
    end
    e = F*randn(ne,NPER);
    for i = 1 : ne
        name = elist{i};
        if NDRAW == 1
            D.(name).data(MAXLAG+1:end,1:NLOOP) = e(i*ones(1,NLOOP),:).';
        else
            D.(name).data(MAXLAG+1:end,iloop) = e(i,:).';
        end
    end
end
end
% xxrandomshocks().