From 241e6d754e71e9784c49fc70f218269b07290b64 Mon Sep 17 00:00:00 2001 From: YanzhenXiangRobotics Date: Wed, 10 May 2023 23:18:37 +0200 Subject: Revert mistakenly pushed to master This reverts commit c5cdd03fc0588f7d34d8b97713caa6d19e5becea. --- templates/cprnd.m | 296 ------------------------------------------------------ 1 file changed, 296 deletions(-) delete mode 100755 templates/cprnd.m (limited to 'templates/cprnd.m') diff --git a/templates/cprnd.m b/templates/cprnd.m deleted file mode 100755 index 4183ac3..0000000 --- a/templates/cprnd.m +++ /dev/null @@ -1,296 +0,0 @@ -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 -- cgit v1.2.1