function IRF = read_irf_bundle(files)
% Reads matrix-style IRFs:
% horizon: Hx1 (or 1xH)
% point/lo/hi: HxN (or NxH)

h  = readmatrix(files.h);
pt = readmatrix(files.pt);

IRF.h  = h(:);
IRF.pt = pt;

IRF.lo68 = readmatrix(files.lo68);
IRF.hi68 = readmatrix(files.hi68);
IRF.lo90 = readmatrix(files.lo90);
IRF.hi90 = readmatrix(files.hi90);

H = numel(IRF.h);

% If horizons are in columns (NxH), transpose everything to HxN
if size(IRF.pt,1) ~= H && size(IRF.pt,2) == H
    IRF.pt   = IRF.pt';
    IRF.lo68 = IRF.lo68';
    IRF.hi68 = IRF.hi68';
    IRF.lo90 = IRF.lo90';
    IRF.hi90 = IRF.hi90';
end

end