function nmin = nfminbnd(fun,n1,n2,itermax,varargin)
%NFMINBND Private routine for optimizing integer complexity /
%regularisation parameters
fmin = inf;
if n2 - n1 < itermax % small interval, try all values
for n = n1:n2
ff = feval(fun,n,varargin{:});
if ff < fmin
nmin = n;
fmin = ff;
end
end
else % large interval, minimize
n0 = round((n2-n1)/2); % start by bounds and half way
N = [n2 n1 n0]; % points
F = [feval(fun,n2,varargin{:}) feval(fun,n1,varargin{:}) ...
feval(fun,n0,varargin{:})]; % function values
ftot = sum(F); % just to have a second stopping criterion
ftotmin = inf; % initially inf
[ff,j] = min(F); % initial minimum
nmin = N(j); % at this point
iter = 3; % that took already 3 function evaluations!
% now loop as long as we improve for at
% most itermax iterations
while (iter < itermax) & ((ff < fmin) | (ftot < ftotmin))
fmin = ff; % store our minimum function
ftotmin = sum(F); % and our set of function values
X = [N.^2; N; [1 1 1]]; % try quadratic approximation
a = prinv(X')*F'; % its coefficients
if a(1) <= 0 & a(2) < 0 % wrong !! point to maximum
nnew = nmin+1; % go right of minimum
elseif a(1) <= 0 & a(2) >= 0 % wrong !!
nnew = nmin-1; % go left of minimum
else
nnew = round(-a(2)/(2*a(1))); % new minimum of parabola
end
if nnew <= n1 % if it is out of bound
nnew = min(N) + 1; % put it close to what we have
elseif nnew >= n2
nnew = max(N) - 1;
end
if any(N == nnew) % we found a point we have
j = find(N == nnew);
F(j) = ff; % this will force a break
else
fnew = feval(fun,nnew,varargin{:}); % new point, evaluate
[fmax,j] = max(F); % find the worst point we have
if fnew >= fmax % new point is worse
iter = itermax; % force end of search
else % new point is better
F(j) = fnew; % take it instead
N(j) = nnew;
end
end
[ff,j] = min(F); % minimum
ftot = sum(F); % total sum
nmin = N(j); % best point
end
end
return