function [lpd, dlpd, d2lpd, info] = density(o, x)

% Evaluate the logged prior density at x.
%
% INPUTS
% - o       [dprior]
% - x       [double]   m×1 vector, point where the prior density is evaluated.
%
% OUTPUTS
% - lpd     [double]   scalar, value of the logged prior density at x.
% - dlpd    [double]   m×1 vector, first order derivatives.
% - d2lpd   [double]   m×1 vector, second order derivatives.
% - info    [double]   vector of indices of parameters for which the density is zero or infinite (only returned if nargout==4).
%
% REMARKS
% Second order derivatives holder, d2lpd, has the same rank and shape than dlpd because the priors are
% independent (we would have to use a matrix if non orthogonal priors were allowed in Dynare).
%
% EXAMPLE
%
% >> Prior = dprior();
% >> Prior.add('alpha', 'shape', 'beta', 'mean', 0.5, 'std', 0.1);
% >> lpd = Prior.density(x)

% Copyright © 2023-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/>.

lpd = 0.0;
if nargout>1
    dlpd = zeros(1, length(x));
    if nargout>2
        d2lpd = dlpd;
        if nargout>3
            m = length(x);
            info = NaN(m, 1);
            ninfo = 0;
        end
    end
end
if o.has('uniform')
    oob_uniform = (x(o.iduniform)-o.p3(o.iduniform)<0) | (x(o.iduniform)-o.p4(o.iduniform)>0);
    if any(oob_uniform)
        lpd = -Inf ;
        if nargout==4
            idx = o.iduniform(oob_uniform);
            info(ninfo+1:ninfo+length(idx)) = idx;
            ninfo = ninfo + length(idx);
        else
            return
        end
    else
        lpd = lpd - sum(log(o.p4(o.iduniform)-o.p3(o.iduniform))) ;
        if nargout>1
            dlpd(o.iduniform) = zeros(length(o.iduniform), 1);
            if nargout>2
                d2lpd(o.iduniform) = zeros(length(o.iduniform), 1);
            end
        end
    end
end
% Dispatch table for continuous distributions: {index_vector, @lpdf_fn, {args}}.
% Anonymous function handles propagate nargout to the underlying lpdf* function.
families = {};
if o.has('gaussian')
    id = o.idgaussian;
    families{end+1} = {id, @lpdfnorm, {x(id), o.p6(id), o.p7(id)}};
end
if o.has('gamma')
    id = o.idgamma;
    families{end+1} = {id, @lpdfgam, {x(id)-o.p3(id), o.p6(id), o.p7(id)}};
end
if o.has('beta')
    id = o.idbeta;
    families{end+1} = {id, @lpdfgbeta, {x(id), o.p6(id), o.p7(id), o.p3(id), o.p4(id)}};
end
if o.has('invgamma1')
    id = o.idinvgamma1;
    families{end+1} = {id, @lpdfig1, {x(id)-o.p3(id), o.p6(id), o.p7(id)}};
end
if o.has('invgamma2')
    id = o.idinvgamma2;
    families{end+1} = {id, @lpdfig2, {x(id)-o.p3(id), o.p6(id), o.p7(id)}};
end
if o.has('weibull')
    id = o.idweibull;
    families{end+1} = {id, @lpdfgweibull, {x(id), o.p6(id), o.p7(id), o.p3(id)}};
end
for k = 1:length(families)
    id = families{k}{1};
    fn = families{k}{2};
    args = families{k}{3};
    switch nargout
      case 1
        lpd = lpd + sum(fn(args{:}));
        if isinf(lpd), return, end
      case 2
        [tmp, dlpd(id)] = fn(args{:});
        lpd = lpd + sum(tmp);
        if isinf(lpd), return, end
      case 3
        [tmp, dlpd(id), d2lpd(id)] = fn(args{:});
        lpd = lpd + sum(tmp);
        if isinf(lpd), return, end
      case 4
        [tmp, dlpd(id), d2lpd(id)] = fn(args{:});
        lpd = lpd + sum(tmp);
        if any(isinf(tmp))
            idx = id(isinf(tmp));
            info(ninfo+1:ninfo+length(idx)) = idx;
            ninfo = ninfo + length(idx);
        end
    end
end
if nargout>3
    info = info(1:ninfo);
end

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

%@test:1
% Build a 14-parameter prior (2 of each distribution family) with add(),
% then verify that density at the mode dominates density at random draws.
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',    'mean', 0.4, 'std', 0.2, '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
    % Density at the prior mode dominates draws.
    lpdstar = Prior.density(Prior.p5);
    LPD = NaN(10000,1);
    parfor i = 1:10000
        x = Prior.draw();
        LPD(i) = Prior.density(x);
    end
    t(1) = true;
catch
    t(1) = false;
end

if t(1)
    t(2) = all(LPD<=lpdstar);
end
T = all(t);
%@eof:1

%@test:2
% Build a 14-parameter prior (2 of each distribution family) with add(),
% then verify that conditional marginal moments match the declared mean/std.
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',    'mean', 0.4, 'std', 0.2, '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
    n = length(Prior);
    mu = NaN(n,1);
    std = NaN(n,1);
    for i=1:n
        % Define conditional density (it's also a marginal since priors are orthogonal)
        f = @(x) exp(Prior.densities(substitute(Prior.p5, i, x)));
        % TODO: Check the version of Octave we use (integral is available as a compatibility wrapper in latest Octave version)
        m = integral(f, Prior.p3(i), Prior.p4(i));
        density = @(x) f(x)/m; % rescaling is required since the probability mass depends on the conditioning.
        mu(i) = integral(@(x) x.*density(x), Prior.p3(i), Prior.p4(i));
        std(i) = sqrt(integral(@(x) ((x-mu(i)).^2).*density(x), Prior.p3(i), Prior.p4(i)));
    end
    t(1) = true;
catch
    t(1) = false;
end

if t(1)
    t(2) = all(abs(mu-.4)<1e-6);
    t(3) = all(abs(std-.2)<1e-6);
end
T = all(t);
%@eof:2

%@test:3
% Regression test: location parameter (p3) shift invariance.
% For every distribution that uses p3 as a location shift in density.m
% (Gamma, InvGamma1, InvGamma2, Weibull), evaluating at (x + delta) with
% (p3 + delta) must give the same log-density and gradient as at x with p3.
%
% Note on multi-output calls: dprior's subsref overload returns a single
% output, so the 2-output form must use density(o, x) rather than o.density(x).

delta = 0.7;
try
    %           Gamma  InvGamma1  InvGamma2  Weibull
    p3_base = [0.5;   1.0;       1.5;       2.0  ];
    h1      = [3.0;   2.0;       2.0;       2.0  ];  % α, s, s, k
    h2      = [0.5;   5.0;       6.0;       1.0  ];  % β, ν, ν, λ

    z   = [1.0; 0.5; 0.5; 1.0];  % z = x - p3 (same for A and B)
    x_A = p3_base + z;
    x_B = x_A + delta;             % x shifted; p3 will also shift by delta

    A = dprior();
    A.add('gam',  'shape', 'gamma',      'h1', h1(1), 'h2', h2(1), 'lower_bound', p3_base(1));
    A.add('ig1',  'shape', 'inv_gamma',  'h1', h1(2), 'h2', h2(2), 'lower_bound', p3_base(2));
    A.add('ig2',  'shape', 'inv_gamma2', 'h1', h1(3), 'h2', h2(3), 'lower_bound', p3_base(3));
    A.add('weib', 'shape', 'weibull',    'h1', h1(4), 'h2', h2(4), 'lower_bound', p3_base(4));

    B = dprior();
    B.add('gam',  'shape', 'gamma',      'h1', h1(1), 'h2', h2(1), 'lower_bound', p3_base(1) + delta);
    B.add('ig1',  'shape', 'inv_gamma',  'h1', h1(2), 'h2', h2(2), 'lower_bound', p3_base(2) + delta);
    B.add('ig2',  'shape', 'inv_gamma2', 'h1', h1(3), 'h2', h2(3), 'lower_bound', p3_base(3) + delta);
    B.add('weib', 'shape', 'weibull',    'h1', h1(4), 'h2', h2(4), 'lower_bound', p3_base(4) + delta);

    lpd_A            = A.density(x_A);    % 1-output path
    [lpd_A2, dlpd_A] = density(A, x_A);  % 2-output path
    lpd_B            = B.density(x_B);
    [lpd_B2, dlpd_B] = density(B, x_B);
    t(1) = true;
catch
    t(1) = false;
end
if t(1)
    % 1-output: f(x; p3) == f(x+δ; p3+δ) for all 4 distributions
    t(2) = abs(lpd_A - lpd_B) < 1e-10;
    % 2-output path: same shift invariance holds
    t(3) = abs(lpd_A2 - lpd_B2) < 1e-10;
    % Gradient depends only on z = x - p3 and is therefore shift-invariant
    t(4) = all(abs(dlpd_A - dlpd_B) < 1e-10);
    % Analytical check for Weibull (index 4, z=1, k=2, λ=1):
    %   log f(z; k, λ) = log(k) - log(λ) + (k-1)·log(z/λ) - (z/λ)^k
    %   d/dx log f     = (k-1)/z - k·z^(k-1)/λ^k
    k_w = h1(4); lam_w = h2(4); z_w = z(4);
    expected_lpd_w  = log(k_w) - log(lam_w) + (k_w-1).*log(z_w./lam_w) - (z_w./lam_w).^k_w;
    expected_dlpd_w = (k_w-1)./z_w - k_w.*z_w.^(k_w-1)./lam_w.^k_w;
    % Isolate Weibull contribution via a single-distribution prior
    W = dprior('weib', 'shape', 'weibull', 'h1', h1(4), 'h2', h2(4), 'lower_bound', p3_base(4));
    [lpd_W, dlpd_W] = density(W, x_A(4));
    t(5) = abs(lpd_W  - expected_lpd_w)  < 1e-10;
    t(6) = abs(dlpd_W - expected_dlpd_w) < 1e-10;
end
T = all(t);
%@eof:3
