Skip to content
Snippets Groups Projects
Verified Commit 09e14f3f 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.

(cherry picked from commit 1912f677)
parent 4fa32017
No related branches found
No related tags found
1 merge request!1815WIP Cherry-picks for 4.6
...@@ -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