function [pars, post] = get_nested_post(pfile, npars, wp, samps)

% function [pars, post] = get_nested_post(pfile, npars, wp, samps)
%
% This function is for use with the output of the MultiNest code to produce
% posteriors from the output samples. As an input it takes in the filename
% of the MutliNest output *-.txt file, which should have 2+npars columns:
% the first is the sample prior mass weigthed likelihood normalised by the 
% evidence, the second is -2* log likelihood, and the rest are the 
% parameter values for each of the samples. If the value of npars is not 
% the same as the number in file an error will occur.
%
% wp is a vector containing the parameter(s) for which you want to create
% the posterior e.g. wp=1 then you will get a posterior only on the 1st
% parameter, or wp=[1 3] will give you a 2D posterior on the 1st and 3rd
% parameters. For the moment this can only contain a maximum of three
% parameters
%
% samps specifies the griding of the output posterior.

% load pfile
d = load(pfile);

s = size(d);

% check that npars is correct
if npars ~= s(2)-2
    error('Error... npars is not correct for the given file!');
end

lwp = length(wp);

% check the length of the vector containing which parameters
if lwp > npars || max(wp) > npars || lwp > 3 || lwp < 1
    error('Error... output parameters are out of range!');
end

% allocate memory for output posterior
if lwp == 1
    post = zeros(samps,1);
else
    post = zeros(samps*ones(1,lwp));
end

% set parameter vectors (using the maximum and minumum values)
partemp = zeros(samps+1, lwp);
pars = zeros(samps, lwp);

for i=1:lwp
    partemp(:,i) = linspace(min(d(:,wp(i)+2)), max(d(:,wp(i)+2)), samps+1);
    dp = (partemp(2,i)-partemp(1,i))/2;
    pars(:,i) = partemp(1:end-1,i) + dp;
end

% get the posterior values in a grid by summing the values output by
% MultiNest within each grid point
for i=1:samps
    if lwp == 1
        pos1 = wp(1)+2;
        post(i) = sum(d(d(:,pos1) > partemp(i,1) & ...
            d(:,pos1) < partemp(i+1,1),1));
    elseif lwp > 1
        pos1 = wp(1)+2;
        pos2 = wp(2)+2;
    
        if lwp == 2
            for j=1:samps
                post(i,j) = sum(d(d(:,pos1) > partemp(i,1) & ...
                    d(:,pos1) < partemp(i+1,1) & d(:,pos2) > ...
                    partemp(j,2) & d(:,pos2) < partemp(j+1,2),1));
            end
        elseif lwp ==3
            pos3 = wp(3)+2;
            
            for j=1:samps
                for k=1:samps
                    post(i,j,k) = sum(d(d(:,pos1) > ...
                        partemp(i,1) & d(:,pos1) < ...
                        partemp(i+1,1) & d(:,pos2) > ...
                        partemp(j,2) & d(:,pos2) < ...
                        partemp(j+1,2) & d(:,pos3) > ...
                        partemp(k,3) & d(:,pos3) < ...
                        partemp(k+1,3),1));
                end
            end
        end
    end
end

if lwp == 1
    plot(pars, post);
elseif lwp == 2
    surf(pars(:,1), pars(:,2), post');
    lighting phong;
end