function [s,d] = myparamstruct(this,d,penalty,varargin)
% MYPARAMSTRUCT  [Not a public function] Parse structure with parameter estimation specs.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox. 
% -Copyright (c) 2007-2012 Jaromir Benes.

isfromstruct = ~any(strcmp(varargin,'model'));

%**************************************************************************

list = fieldnames(d).';
nlist = length(list);
remove = false(1,nlist);
for i = 1 : nlist
    if isempty(d.(list{i}))
        remove(i) = true;
    end
end
d = rmfield(d,list(remove));
list(remove) = [];

[assignpos,stdcorrpos,assignval,stdcorrval,type] = ...
    mynameposition(this,list);

s = struct();

% Parameters to estimate.
s.plist = list(type == 4);
s.assignpos  = assignpos(type == 4);
s.stdcorrpos = stdcorrpos(type == 4);
np = length(s.plist);

s.p0 = nan(1,np);
s.pl = nan(1,np);
s.pu = nan(1,np);
s.prior = cell(1,np);
s.priorindex = false(1,np);

validbounds = true(1,np);
withinbounds = true(1,np);
doparameters();
dochkbounds();

% Remove parameter fields and return a struct with non-parameter fields.
d = rmfield(d,s.plist);

% Nested functions.

%**************************************************************************
    function doparameters()
        for ii = 1 : np
            name = s.plist{ii};
            spec = d.(name);
            if isnumeric(spec)
                spec = num2cell(spec);
            end
            % Starting value.
            if isfromstruct && ~isempty(spec) && isnumericscalar(spec{1})
                p0 = spec{1};
            else
                p0 = NaN;
            end
            % If not supplied by the user, use the currently assigned value.
            if isnan(p0)
                if ~isnan(assignpos(ii))
                    p0 = assignval(1,ii,:);
                else
                    p0 = stdcorrval(1,ii,:);
                end
            end
            % Lower bound.
            if length(spec) > 1 && isnumericscalar(spec{2})
                pl = spec{2};
            else
                pl = -Inf;
            end
            % Upper bound.
            if length(spec) > 2  && isnumericscalar(spec{3})
                pu = spec{3};
            else
                pu = Inf;
            end
            % Check that the lower bound is really lower than the upper bound.
            if pl >= pu
                validbounds(ii) = false;
                continue
            end
            % Check that the starting values in within the bounds.
            if p0 < pl || p0 > pu
                withinbounds(ii) = false;
                continue
            end
            
            % Prior distribution function (function_handle), or penalty function
            % ([weight] or [weight,pbar]).
            isprior = false;
            prior = [];
            if length(spec) > 3 && ~isempty(spec{4})
                isprior = true;
                if isa(spec{4},'function_handle')
                    % The 4th entry is a prior distribution function handle.
                    prior = spec{4};
                elseif isnumeric(spec{4}) && penalty > 0
                    % The 4th entry is a penalty function, compute the total weight including
                    % the `'penalty='` option.
                    totalweight = spec{4}(1)*penalty;
                    if length(spec{4}) == 1
                        % Only the weight specified. The centre of penalty function is then set
                        % identical to the starting value.
                        pbar = p0;
                    else
                        % Both the weight and the centre specified.
                        pbar = spec{4}(2);
                    end
                    if isnan(pbar)
                        if ~isnan(assignpos(ii))
                            pbar = assignval(1,ii,:);
                        else
                            pbar = stdcorrval(1,ii,:);
                        end
                    end
                    % Convert penalty function to a normal prior: w*(p - pbar)^2 == 1/2*((p -
                    % pbar)/sgm)^2 => sgm = 1/sqrt(2*w).
                    sgm = 1/sqrt(2*totalweight);
                    prior = logprior.normal(pbar,sgm);
                end
            end
            
            s.p0(ii) = p0;
            s.pl(ii) = pl;
            s.pu(ii) = pu;
            s.prior{ii} = prior;
            s.priorindex(ii) = isprior;
            
        end
    end
% doparameters().

%**************************************************************************
    function dochkbounds()
        if any(~validbounds)
            utils.error(class(this), ...
                ['Lower and upper bounds for this parameter ', ...
                'are inconsistent: ''%s''.'], ....
                s.plist{~validbounds});
        end
        
        if any(~withinbounds)
            utils.error(class(this), ...
                ['Starting value for this parameter is ', ...
                'outside the specified bounds: ''%s''.'], ....
                s.plist{~withinbounds});
        end
    end
% dochkbounds().

end