Skip to content
Snippets Groups Projects
Verified Commit 1912f677 authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

MATLAB compatibility fix: automatic broadcasting was introduced in R2016b

For earlier versions, either use bsxfun or handle special cases differently.
parent 8fff9911
Branches
Tags
No related merge requests found
...@@ -27,7 +27,11 @@ nterms = size(z,2); ...@@ -27,7 +27,11 @@ nterms = size(z,2);
for k=1:size(y,1) for k=1:size(y,1)
ytmp = squeeze(y(k,:,:)); ytmp = squeeze(y(k,:,:));
yres = ytmp(end,:) - sum(ytmp(1:end-1,:)); yres = ytmp(end,:) - sum(ytmp(1:end-1,:));
if ~isoctave && matlab_ver_less_than('9.1') % Automatic broadcasting was introduced in MATLAB R2016b
w = bsxfun(@rdivide, abs(ytmp(1:end-1,:)), sum(abs(ytmp(1:end-1,:))))
else
w = abs(ytmp(1:end-1,:))./sum(abs(ytmp(1:end-1,:))); w = abs(ytmp(1:end-1,:))./sum(abs(ytmp(1:end-1,:)));
end
% ytmp(1:end-1,:) = ytmp(1:end-1,:) + repmat(yres,[nterms-1 1])/(nterms-1); % ytmp(1:end-1,:) = ytmp(1:end-1,:) + repmat(yres,[nterms-1 1])/(nterms-1);
ytmp(1:end-1,:) = ytmp(1:end-1,:) + repmat(yres,[nterms-1 1]).*w; ytmp(1:end-1,:) = ytmp(1:end-1,:) + repmat(yres,[nterms-1 1]).*w;
y(k,:,:) = ytmp; y(k,:,:) = ytmp;
......
...@@ -9,7 +9,7 @@ function g = best_1978(a ,b) ...@@ -9,7 +9,7 @@ function g = best_1978(a ,b)
% OUTPUTS % OUTPUTS
% - g [double] n*1 vector, gamma variates. % - g [double] n*1 vector, gamma variates.
% Copyright (C) 2006-2018 Dynare Team % Copyright (C) 2006-2020 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -44,7 +44,12 @@ while mm ...@@ -44,7 +44,12 @@ while mm
X(index) = bb(index)+Y(index); % x X(index) = bb(index)+Y(index); % x
id1 = index(X(index)<0); % Reject. id1 = index(X(index)<0); % Reject.
id2 = setdiff(index, id1); id2 = setdiff(index, id1);
if numel(id2) ~= 0 || isoctave || ~matlab_ver_less_than('9.1')
% If id2=[], LHS of the .* has size [0,0], while RHS is [0,1].
% Since there is no automatic broadcast in MATLAB < R2016b, skip the
% statement in that case.
Z(id2) = 64.0*(W(id2).^3).*(rand(length(id2),1).^2); % d Z(id2) = 64.0*(W(id2).^3).*(rand(length(id2),1).^2); % d
end
id3 = id2(Z(id2)>1.0-2.0*Y(id2).*Y(id2)./X(id2)); % Reject. id3 = id2(Z(id2)>1.0-2.0*Y(id2).*Y(id2)./X(id2)); % Reject.
id4 = id3(log(Z(id3))>2.0*(bb(id3).*log(X(id3)./bb(id3))-Y(id3))); % Reject. id4 = id3(log(Z(id3))>2.0*(bb(id3).*log(X(id3)./bb(id3))-Y(id3))); % Reject.
index = [id1, id4]; index = [id1, id4];
......
...@@ -9,7 +9,7 @@ function g = knuth(a, b) ...@@ -9,7 +9,7 @@ function g = knuth(a, b)
% OUTPUTS % OUTPUTS
% - g [double] n*1 vector, gamma variates. % - g [double] n*1 vector, gamma variates.
% Copyright (C) 2006-2018 Dynare Team % Copyright (C) 2006-2020 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -38,7 +38,15 @@ while mm ...@@ -38,7 +38,15 @@ while mm
X(index) = Y(index).*bb(index) + a(index) - 1; X(index) = Y(index).*bb(index) + a(index) - 1;
id1 = index(X(index)<=0); % Rejected draws. id1 = index(X(index)<=0); % Rejected draws.
id2 = setdiff(index, id1); id2 = setdiff(index, id1);
if numel(id2) == 0 && ~isoctave && matlab_ver_less_than('9.1')
% The LHS of the > comparison in the "else" branch has size = [0, 1]
% while the RHS has size = [0, 0].
% Automatic broadcasting was only introduced in MATLAB R2016b, so we
% must handle this case separately to avoid an error.
id3 = [];
else
id3 = id2(rand(length(id2), 1)>(1+Y(id2).*Y(id2)).*exp((a(id2)-1).*(log(X(id2))-log(a(id2)-1))-bb(id2).*Y(id2))); % Rejected draws. id3 = id2(rand(length(id2), 1)>(1+Y(id2).*Y(id2)).*exp((a(id2)-1).*(log(X(id2))-log(a(id2)-1))-bb(id2).*Y(id2))); % Rejected draws.
end
index = [id1, id3]; index = [id1, id3];
mm = length(index); mm = length(index);
end end
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment