1function [theta, fxsim, neval] = rotated_slice_sampler(objective_function,theta,thetaprior,sampler_options,varargin)
2% ----------------------------------------------------------
3% ROTATED SLICE SAMPLER - with stepping out (Neal, 2003)
4% extension of the orthogonal univarite sampler (slice_sampler.m)
5% copyright M. Ratto (European Commission)
6%
7% objective_function(theta,varargin): -log of any unnormalized pdf
8% with varargin (optional) a vector of auxiliaty parameters
9% to be passed to f( ).
10% ----------------------------------------------------------
11%
12% INPUTS
13%   objective_function:       objective function (expressed as minus the log of a density)
14%   theta:                    last value of theta
15%   thetaprior:               bounds of the theta space
16%   sampler_options:          posterior sampler options
17%   varargin:                 optional input arguments to objective function
18%
19% OUTPUTS
20%   theta:       new theta sample
21%   fxsim:       value of the objective function for the new sample
22%   neval:       number of function evaluations
23%
24% SPECIAL REQUIREMENTS
25%   none
26
27% Copyright (C) 2015-2017 Dynare Team
28%
29% This file is part of Dynare.
30%
31% Dynare is free software: you can redistribute it and/or modify
32% it under the terms of the GNU General Public License as published by
33% the Free Software Foundation, either version 3 of the License, or
34% (at your option) any later version.
35%
36% Dynare is distributed in the hope that it will be useful,
37% but WITHOUT ANY WARRANTY; without even the implied warranty of
38% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
39% GNU General Public License for more details.
40%
41% You should have received a copy of the GNU General Public License
42% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
43
44theta=theta(:);
45npar = length(theta);
46neval = zeros(npar,1);
47W1=[];
48if isfield(sampler_options,'WR')
49    W1 = sampler_options.WR;
50end
51if ~isempty(sampler_options.mode)
52    mm = sampler_options.mode;
53    n = length(mm);
54    for j=1:n
55        distance(j)=sqrt(sum((theta-mm(j).m).^2));
56    end
57    [m, im] = min(distance);
58
59    r=im;
60    V1 = mm(r).m;
61    jj=0;
62    for j=1:n
63        if j~=r
64            jj=jj+1;
65            tmp=mm(j).m-mm(r).m;
66            %tmp=mm(j).m-theta;
67            V1(:,jj)=tmp/norm(tmp);
68        end
69    end
70    resul=randperm(n-1,n-1);
71    V1 = V1(:,resul);
72
73    %V1 = V1(:, randperm(n-1));
74    %     %d = chol(mm(r).invhess);
75    %     %V1 = transpose(feval(sampler_options.proposal_distribution, transpose(mm(r).m), d, npar));
76    %
77    %     V1=eye(npar);
78    %     V1=V1(:,randperm(npar));
79    %     for j=1:2,
80    %         V1(:,j)=mm(r(j)).m-theta;
81    %         V1(:,j)=V1(:,j)/norm(V1(:,j));
82    %     end
83    %     % Gram-Schmidt
84    %     for j=2:npar,
85    %         for k=1:j-1,
86    %             V1(:,j)=V1(:,j)-V1(:,k)'*V1(:,j)*V1(:,k);
87    %         end
88    %         V1(:,j)=V1(:,j)/norm(V1(:,j));
89    %     end
90    %     for j=1:n,
91    %         distance(j)=sqrt(sum((theta-mm(j).m).^2));
92    %     end
93    %     [m, im] = min(distance);
94    %     if im==r,
95    %         fxsim=[];
96    %         return,
97    %     else
98    %         theta1=theta;
99    %     end
100else
101    V1 = sampler_options.V1;
102end
103npar=size(V1,2);
104
105for it=1:npar
106    theta0 = theta;
107    neval(it) = 0;
108    xold  = 0;
109    % XLB   = thetaprior(3);
110    % XUB   = thetaprior(4);
111    tb=sort([(thetaprior(:,1)-theta)./V1(:,it) (thetaprior(:,2)-theta)./V1(:,it)],2);
112    XLB=max(tb(:,1));
113    XUB=min(tb(:,2));
114    if isempty(W1)
115        W = (XUB-XLB); %*0.8;
116    else
117        W = W1(it);
118    end
119
120    % -------------------------------------------------------
121    % 1. DRAW Z = ln[f(X0)] - EXP(1) where EXP(1)=-ln(U(0,1))
122    %    THIS DEFINES THE SLICE S={x: z < ln(f(x))}
123    % -------------------------------------------------------
124
125    fxold = -feval(objective_function,theta,varargin{:});
126    %I have to be sure that the rotation is for L,R or for Fxold, theta(it)
127    neval(it) = neval(it) + 1;
128    Z = fxold + log(rand(1,1));
129    % -------------------------------------------------------------
130    % 2. FIND I=(L,R) AROUND X0 THAT CONTAINS S AS MUCH AS POSSIBLE
131    %    STEPPING-OUT PROCEDURE
132    % -------------------------------------------------------------
133    u = rand(1,1);
134    L = max(XLB,xold-W*u);
135    R = min(XUB,L+W);
136
137    %[L R]=slice_rotation(L, R, alpha);
138    while(L > XLB)
139        xsim = L;
140        theta = theta0+xsim*V1(:,it);
141        fxl = -feval(objective_function,theta,varargin{:});
142        neval(it) = neval(it) + 1;
143        if (fxl <= Z)
144            break
145        end
146        L = max(XLB,L-W);
147    end
148    while(R < XUB)
149        xsim = R;
150        theta = theta0+xsim*V1(:,it);
151        fxr = -feval(objective_function,theta,varargin{:});
152        neval(it) = neval(it) + 1;
153        if (fxr <= Z)
154            break
155        end
156        R = min(XUB,R+W);
157    end
158    % ------------------------------------------------------
159    % 3. SAMPLING FROM THE SET A = (I INTERSECT S) = (LA,RA)
160    % ------------------------------------------------------
161    fxsim = Z-1;
162    while (fxsim < Z)
163        u = rand(1,1);
164        xsim = L + u*(R - L);
165        theta = theta0+xsim*V1(:,it);
166        fxsim = -feval(objective_function,theta,varargin{:});
167        neval(it) = neval(it) + 1;
168        if (xsim > xold)
169            R = xsim;
170        else
171            L = xsim;
172        end
173    end
174end
175
176% if ~isempty(sampler_options.mode),
177%     dist1=sqrt(sum((theta-mm(r).m).^2));
178%     if dist1>distance(r),
179%         theta=theta1;
180%         fxsim=[];
181%     end
182% end
183end
184