function [samples,seed] = ars1d(fun, dims, opts)
% ARS1D   1-D derivative-free adaptive rejection sampler
%
%   [samples,seed] = ars1d(fun, dims, opts)
%
%  Efficient 1-D sampler for log-convex density functions. While sampling
%  , computes and refines a piecewise-exponential approximation to the 
%  density function. Note that this MEX interface does not save this
%  approximation between calls; it's best to sample as many variates as you
%  think you'll need up front.
%
%  The fun argument is the density function you wish to sample. It should be
%  be log concave and non-negative everywhere (unless it returns the log of the
%  density; see below).
%
%  The dims argument is optional; if blank, the sampler returns one sample.
%  Otherwise, the argument works the same as the argument in single-argument
%  calls to MATLAB's rand function.
%
%  The optional opts argument is a struct containing arguments for the sampler.
%  If opts is not specified, default options are used. Calling this function
%  with no arguments returns the default opts struct. Fields in opts include:
%
%    x_min			(default: -inf)
%
%     Specifies the left bound of support for the density function.
%     May be infinite.
%
%    x_max			(default: inf)
%
%     Specifies the right bound of support for the density function.
%     May be infinite.
%
%    function_returns_log	(default: false)
%
%     Some distributions may be more numerically stable when expressed as the
%     log of the density function instead of the function itself. Set
%     this field to true iff fun returns the log of the density function.
%
%    seed			(optional)
%
%     The sampler uses a different random number generator than MATLAB's---
%     one provided by the Boost C++ library. Each call to the sampler must
%     seed this separate generator. If this field is missing from the opts
%     struct, a random seed value is generated from MATLAB's random number
%     generator. This way, to ensure identical samples, it's only necessary
%     to control MATLAB's random seed. If you wish to set the sampler's random
%     seed manually, however, specify a floating point number here in [0,1].
%
%    hint1, hint2, hint3	(default for all: -inf)
%
%  Though pains have been taken to avoid it, this code may encounter
%  numerical contingencies that prevent it from sampling the density
%  function. If this happens, the samples result will contain a single NaN
%  value. Sometimes these difficulties can be overcome with more helpful
%  hints; other times, as in the case of betapdf(x,1.000001,2), the
%  resolution of doubles is insufficient for the sampler to locate points
%  around the mode of the density function.
%
%  The extra return value "seed" will contain the random seed used to
%  initialize the sampler's random generator.
%
%  The ars1d.m file is a wrapper for a MEX function, and it performs all of
%  the error checking for the routine. If you're sure your arguments are OK,
%  you can call mex_ars1d directly. With mex_ars1d, all of the optional
%  arguments and parameters for ars1d are mandatory, including the seed field
%  in opts.

% If no arguments, return the default options struct
if nargin == 0
  samples = struct;
  samples.x_min =			-inf;
  samples.x_max =			inf;
  samples.function_returns_log =	false;
  samples.hint1 =			-inf;
  samples.hint2 =			-inf;
  samples.hint3 =			-inf;
  return;
end

% If no options, use the default options struct. Otherwise, sanity check
% the options.
if nargin < 3
  opts = ars1d;
else
  assert(opts.x_min < opts.x_max, 'x_min must be less than x_max in opts');
  assert(isfield(opts,'function_returns_log'), ...
	 'opts missing required field "function_returns_log"');
  assert(islogical(opts.function_returns_log), ...
	 'opts.function_returns_log must be a logical value');
  assert(isfield(opts,'hint1'), 'opts missing required field "hint1"');
  assert(isfield(opts,'hint2'), 'opts missing required field "hint2"');
  assert(isfield(opts,'hint3'), 'opts missing required field "hint3"');
end

% If no random seed, create a random seed
if ~isfield(opts,'seed')
  opts.seed = rand;
end

% If no dimensions, we just want one sample.
if nargin < 2
  dims = 1;
end

% Furnish the random seed we used
if nargout > 1
  seed = opts.seed;
end

% Call the sampler
samples = mex_ars1d(fun, dims, opts);