add submission zipsubmission

1 files changed, 296 insertions, 0 deletions

diff --git a/Submission files/cprnd.m b/Submission files/cprnd.m
new file mode 100644
index 0000000..4183ac3
--- /dev/null
+++ b/Submission files/cprnd.m
@@ -0,0 +1,296 @@
+function [X S] = cprnd(N,A,b,options)
+%CPRND Draw from the uniform distribution over a convex polytope.
+%   X = cprnd(N,A,b) Returns an N-by-P matrix of random vectors drawn
+%   from the uniform distribution over the interior of the polytope
+%   defined by Ax <= b. A is a M-by-P matrix of constraint equation
+%   coefficients. b is a M-by-1 vector of constraint equation
+%   constants.
+%    
+%   cprnd(N,A,b,options) allows various options to be specified.
+
+%   'method'    Specifies the algorithm. One of the strings 'gibbs',
+%               'hitandrun' (the default), and 'achr'. The default
+%               algorithm 'hitandrun' is a vanilla hit-and-run sampler
+%               [1]. 'gibbs' specifies a Gibbs sampler, 'achr' is the
+%               Adaptive Centering Hit-and-Run algorithm of [2]. 
+%
+%   'x0'        Vector x0 is a starting point which should be interior 
+%               to the polytope. If x0 is not supplied CPRND uses the
+%               Chebyshev center as the initial point.
+%    
+%   'isotropic' Perform an adaptive istropic transformation of the 
+%               polytope. The values 0, 1 and 2 respectively turn
+%               off the transformation, construct it during a runup
+%               period only, or continuously update the
+%               tranformation throughout sample production. The
+%               transformation makes sense only for the Gibbs and
+%               Hit-and-run samplers (ACHR is invariant under
+%               linear transformations).
+%   
+%   'discard'   Number of initial samples (post-runup) to discard.
+%
+%   'runup'     When the method is gibbs or hitandrun and the
+%               isotropic transformation is used, this is the
+%               number of initial iterations of the algorithm in
+%               the untransformed space. That is a sample of size
+%               runup is generated and its covariance used as the
+%               basis of a transformation. 
+%               When the method is achr runup is the number of
+%               conventional hitandrun iterations. See [2].
+%
+%   'ortho'     Zero/one flag. If turned on direction vectors for
+%               the hitandrun algorithm are generated in random
+%               orthonormal blocks rather than one by
+%               one. Experimental and of dubious utility.
+%   
+%   'quasi' Allows the user to specify a quasirandom number generator
+%               (such as 'halton' or 'sobol').  Experimental and of
+%               dubious utility.
+%
+%      
+
+%   By default CPRND employs the hit-and-run sampler which may
+%   exhibit marked serial correlation and very long convergence times.
+%
+% References
+% [1] Kroese, D.P. and Taimre, T. and Botev, Z.I., "Handbook of Monte
+% Carlo Methods" (2011), pp. 240-244.
+% [2] Kaufman, David E. and Smith, Robert L., "Direction Choice for
+% Accelerated Convergence in Hit-and-Run Sampling", Op. Res. 46,
+% pp. 84-95.
+%
+% Copyright (c) 2011-2012 Tim J. Benham, School of Mathematics and Physics,
+%                University of Queensland.
+
+    function y = stdize(z)
+        y = z/norm(z);
+    end
+
+    p = size(A,2);                         % dimension
+    m = size(A,1);                         % num constraint ineqs
+    x0 = [];
+    runup = [];                              % runup necessary to method
+    discard = [];                            % num initial pts discarded
+    quasi = 0;
+    method = 'achr';
+    orthogonal = 0;
+    isotropic = [];
+    
+    % gendir generates a random unit (direction) vector.
+    gendir = @() stdize(randn(p,1));
+
+    % Alternative function ogendir forces directions to come in
+    % orthogonal bunches.
+    Ucache = {};
+    function u = ogendir()
+        if length(Ucache) == 0
+            u = stdize(randn(p,1));
+            m = null(u');               % orthonormal basis for nullspace
+            Ucache = mat2cell(m',ones(1,p-1));
+        else
+            u = Ucache{end}';
+            Ucache(end) = [];
+        end
+    end
+    
+    % Check input arguments
+    
+    if m < p+1                             
+        % obv a prob here
+        error('cprnd:obvprob',['At least ',num2str(p+1),' inequalities ' ...
+                            'required']);
+    end
+
+    if nargin == 4
+        if isstruct(options)
+
+            if isfield(options,'method')
+                method = lower(options.method);
+                switch method
+                  case 'gibbs'
+                  case 'hitandrun'
+                  case 'achr'
+                  otherwise
+                    error('cprnd:badopt',...
+                          ['The method option takes only the ' ...
+                           'values "gibbs", "hitandrun", and "ACHR".']);
+                end
+            end
+            
+            if isfield(options,'isotropic')
+                % Controls application of isotropic transformation,
+                % which seems to help a lot.
+                isotropic = options.isotropic;
+            end
+
+            if isfield(options,'discard')
+                % num. of samples to discard
+                discard = options.discard;
+            end
+
+            if isfield(options,'quasi')
+                % Use quasi random numbers, which doesn't seem to
+                % help much.
+                quasi = options.quasi;
+                if quasi && ~ischar(quasi), quasi='halton'; end
+                if ~strcmp(quasi,'none')
+                    qstr = qrandstream(quasi,p,'Skip',1);
+                    gendir = @() stdize(norminv(qrand(qstr,1),0,1)');
+                end
+            end
+
+            if isfield(options,'x0')
+                % initial interior point
+                x0 = options.x0;
+            end
+
+            if isfield(options,'runup')
+                % number of points to generate before first output point
+                runup = options.runup;
+            end
+            
+            if isfield(options,'ortho')
+                % Generate direction vectors in orthogonal
+                % groups. Seems to help a little.
+                orthogonal = options.ortho;
+            end
+            
+        else
+            x0 = options;                   % legacy support
+        end
+    end
+
+    % Default and apply options
+    
+    if isempty(isotropic)
+        if ~strcmp(method,'achr')
+            isotropic = 2;
+        else
+            isotropic = 0;
+        end
+    end
+    
+    if orthogonal
+        gendir = @() ogendir();
+    end
+    
+    % Choose a starting point x0 if user did not provide one.
+    if isempty(x0)
+        x0 = chebycenter(A,b); % prob. if p==1?
+    end
+    
+    % Default the runup to something arbitrary.
+    if isempty(runup)
+        if strcmp(method,'achr')
+            runup = 10*(p+1);
+        elseif isotropic > 0
+            runup = 10*p*(p+1);
+        else 
+            runup = 0;
+        end
+    end
+
+    % Default the discard to something arbitrary
+    if isempty(discard)
+        if strcmp(method,'achr')
+            discard = 25*(p+1);
+        else
+            discard = runup;
+        end
+    end
+
+    X = zeros(N+runup+discard,p);
+
+    n = 0;                                  % num generated so far
+    x = x0;
+
+    % Initialize variables for keeping track of sample mean, covariance
+    % and isotropic transform.
+    M = zeros(p,1);                         % Incremental mean.
+    S2 = zeros(p,p);                        % Incremental sum of
+                                            % outer prodcts.
+    S = eye(p); T = eye(p); W = A;
+    
+    while n < N+runup+discard
+        y = x;
+        
+        % compute approximate stochastic transformation
+        if isotropic>0
+            if n == runup || (isotropic > 1 && n > runup)
+                T = chol(S,'lower');
+                W = A*T;
+            end
+            y = T\y;
+        end
+        
+        switch method
+            
+          case 'gibbs'
+            % choose p new components
+            for i = 1:p
+                % Find points where the line with the (p-1) components x_i
+                % fixed intersects the bounding polytope.
+                e = circshift(eye(p,1),i-1);
+                z = W*e;
+                c = (b - W*y)./z;
+                tmin = max(c(z<0)); tmax = min(c(z>0));
+                % choose a point on that line segment
+                y = y + (tmin+(tmax-tmin)*rand)*e;
+            end
+            
+          case 'hitandrun'
+            % choose a direction
+            u = gendir();
+            % determine intersections of x + ut with the polytope
+            z = W*u;
+            c = (b - W*y)./z;
+            tmin = max(c(z<0)); tmax = min(c(z>0));
+            % choose a point on that line segment
+            y = y + (tmin+(tmax-tmin)*rand)*u;
+            
+          case 'achr'
+            % test whether in runup or not
+            if n < runup
+                % same as hitandrun
+                u = gendir();
+            else
+                % choose a previous point at random
+                v = X(randi(n),:)';
+                % line sampling direction is from v to sample mean
+                u = (v-M)/norm(v-M);
+            end
+            % proceed as in hit and run
+            z = A*u;
+            c = (b - A*y)./z;
+            tmin = max(c(z<0)); tmax = min(c(z>0));
+            % Choose a random point on that line segment
+            y = y + (tmin+(tmax-tmin)*rand)*u;
+        end
+        
+        if isotropic>0
+            x = T*y;
+        else
+            x = y;
+        end
+        
+        X(n+1,:)=x';
+        n = n + 1;
+        
+        % Incremental mean and covariance updates
+        delta0 = x - M;       % delta new point wrt old mean
+        M = M + delta0/n;     % sample mean
+        delta1 = x - M;       % delta new point wrt new mean
+        if n > 1
+            S2 = S2 + (n-1)/(n*n)*delta0*delta0'...
+                 + delta1*delta1';
+            S0 = S;
+            S = S2/(n-1);           % sample covariance
+        else 
+            S = eye(p);
+        end
+        
+    end
+        
+    X = X((discard+runup+1):(N+discard+runup),:); 
+
+end