function s = subset(o, idx)

% Return a new dprior object restricted to the selected parameter indices.
%
% INPUTS
% - o    [dprior]            Prior object.
% - idx  [logical|integer]   Parameter selection (logical mask or index vector).
%
% OUTPUTS
% - s    [dprior]            Prior object for the selected parameters.

% Copyright © 2026 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.

if islogical(idx)
    idx = find(idx);
end

s = dprior();

% Slice parameter vectors.
if ~isempty(o.p1), s.p1 = o.p1(idx); end
if ~isempty(o.p2), s.p2 = o.p2(idx); end
if ~isempty(o.p3), s.p3 = o.p3(idx); end
if ~isempty(o.p4), s.p4 = o.p4(idx); end
if ~isempty(o.p5), s.p5 = o.p5(idx); end
if ~isempty(o.p6), s.p6 = o.p6(idx); end
if ~isempty(o.p7), s.p7 = o.p7(idx); end
if ~isempty(o.p11), s.p11 = o.p11(idx); end
if ~isempty(o.lb), s.lb = o.lb(idx); end
if ~isempty(o.ub), s.ub = o.ub(idx); end
if ~isempty(o.name), s.name = o.name(idx); end

% Remap distribution index vectors.  Original indices refer to positions
% in the full parameter vector; we need positions in the subset (1..n).
names = {'idbeta', 'idgamma', 'idgaussian', 'idinvgamma1', 'iduniform', 'idinvgamma2', 'idweibull'};
for k = 1:length(names)
    [member, pos] = ismember(o.(names{k}), idx);
    s.(names{k}) = pos(member);
end

return % --*-- Unit tests --*--

%@test:1
% Setup: build a 14-parameter prior (2 of each distribution family)
% incrementally with add().
prior_trunc = 1e-10;
Prior = dprior();
for rep = 1:2
    Prior.add(sprintf('p%d_beta',      rep), 'shape', 'beta',       'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_gamma',     rep), 'shape', 'gamma',      'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_normal',    rep), 'shape', 'normal',     'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_invgamma1', rep), 'shape', 'inv_gamma',  'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_uniform',   rep), 'shape', 'uniform',    'lower_bound', 0, 'upper_bound', 1, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_invgamma2', rep), 'shape', 'inv_gamma2', 'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_weibull',   rep), 'shape', 'weibull',    'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
end

try
    idx = (1:7)';
    Sub = Prior.subset(idx);
    t(1) = true;
catch
    t(1) = false;
end

if t(1)
    t(2) = length(Sub)==7;
    t(3) = length(Sub.idbeta)==1;
    t(4) = length(Sub.idgamma)==1;
    t(5) = length(Sub.idgaussian)==1;
    t(6) = length(Sub.idinvgamma1)==1;
    t(7) = length(Sub.iduniform)==1;
    t(8) = length(Sub.idinvgamma2)==1;
    t(9) = length(Sub.idweibull)==1;
    t(10) = Sub.idbeta==1;
    t(11) = Sub.idgamma==2;
    t(12) = Sub.idgaussian==3;
    t(13) = Sub.idinvgamma1==4;
    t(14) = Sub.iduniform==5;
    t(15) = Sub.idinvgamma2==6;
    t(16) = Sub.idweibull==7;
    t(17) = isequal(Sub.p6, Prior.p6(idx));
    t(18) = isequal(Sub.p7, Prior.p7(idx));
    t(19) = isequal(Sub.lb, Prior.lb(idx));
    t(20) = isequal(Sub.ub, Prior.ub(idx));
end
T = all(t);
%@eof:1

%@test:2
% Test that density on a subset matches density on a fresh prior built with
% only those parameters.
prior_trunc = 1e-10;
Prior = dprior();
for rep = 1:2
    Prior.add(sprintf('p%d_beta',      rep), 'shape', 'beta',       'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_gamma',     rep), 'shape', 'gamma',      'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_normal',    rep), 'shape', 'normal',     'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_invgamma1', rep), 'shape', 'inv_gamma',  'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_uniform',   rep), 'shape', 'uniform',    'lower_bound', 0, 'upper_bound', 1, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_invgamma2', rep), 'shape', 'inv_gamma2', 'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_weibull',   rep), 'shape', 'weibull',    'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
end

try
    x = Prior.draw();
    % Pick beta + gamma + normal parameters (indices 1, 2, 3, 8, 9, 10).
    idx = [1, 2, 3, 8, 9, 10];
    Sub = Prior.subset(idx);
    lpd_sub = Sub.density(x(idx));
    % Reference: rebuild those six parameters in order (beta, gamma, normal,
    % beta, gamma, normal) so the id* mapping matches Sub's.
    Ref = dprior();
    for rep = 1:2
        Ref.add(sprintf('p%d_beta',   rep), 'shape', 'beta',   'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
        Ref.add(sprintf('p%d_gamma',  rep), 'shape', 'gamma',  'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
        Ref.add(sprintf('p%d_normal', rep), 'shape', 'normal', 'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    end
    lpd_ref = Ref.density(x(idx));
    t(1) = true;
catch
    t(1) = false;
end

if t(1)
    t(2) = abs(lpd_sub - lpd_ref) < 1e-12;
end
T = all(t);
%@eof:2

%@test:3
% Test subset with logical indexing. Uses add() per parameter.
prior_trunc = 1e-10;
Prior = dprior();
Prior.add('x1', 'shape', 'normal', 'mean', 0.5, 'std', 0.1, 'prior_trunc', prior_trunc);
Prior.add('x2', 'shape', 'normal', 'mean', 0.3, 'std', 0.1, 'prior_trunc', prior_trunc);
Prior.add('x3', 'shape', 'beta',   'mean', 0.4, 'std', 0.1, 'lower_bound', 0, 'upper_bound', 1, 'prior_trunc', prior_trunc);
Prior.add('x4', 'shape', 'beta',   'mean', 0.6, 'std', 0.1, 'lower_bound', 0, 'upper_bound', 1, 'prior_trunc', prior_trunc);
Prior.add('x5', 'shape', 'gamma',  'mean', 0.8, 'std', 0.2, 'prior_trunc', prior_trunc);

try
    % Logical mask: select the two normal parameters (1, 2).
    mask = [true; true; false; false; false];
    Sub = Prior.subset(mask);
    t(1) = true;
catch
    t(1) = false;
end

if t(1)
    t(2) = length(Sub)==2;
    t(3) = isequal(Sub.idgaussian, [1; 2]);
    t(4) = isempty(Sub.idbeta);
    t(5) = isempty(Sub.idgamma);
    t(6) = isequal(Sub.p6, Prior.p6(mask));
end
T = all(t);
%@eof:3
