text:
function varargout = warpvars_vbmc(varargin)
%WARPVARS Linear and nonlinear transformation of variables.
%
%  TRINFO = PDFTRANS(NVARS,LB,UB) returns the transformation structure
%  TRINFO for a problem with NVARS dimensions and lower/upper bounds
%  respectively LB and UB. LB and UB are either scalars or row arrays that
%  can contain real numbers and Inf's
%  The ordering LB <= UB needs to hold coordinate-wise.
%
%  Variables with lower or upper bounds are transformed via a log transform.
%  Variables with both lower and upper bounds are transformed via a logit
%  transform.
%
%  Y = TRANSVARS(X,'dir',TRINFO) performs direct transform of constrained
%  variables X into unconstrained variables Y according to transformation
%  encoded in structure TRINFO. X must be a N x NVARS array, where N is the
%  number of input data and NVARS is the number of dimensions.
%
%  X = TRANSVARS(Y,'inv',TRINFO) performs inverse transform of unconstrained
%  variables Y into constrained variables X.
%
%  P = TRANSVARS(Y,'prob',TRINFO) returns probability multiplier for the
%  original pdf evaluated at f^{-1}(Y), that is | df^{-1}(y) / dy |.
%
%  LP = TRANSVARS(Y,'logprob',TRINFO) returns log probability term for the
%  original log pdf evaluated at f^{-1}(Y).

%  Author: Luigi Acerbi
%  e-mail: luigi.acerbi@gmail.com

if nargin < 3
    error('TRANSVARS requires a minimum of three input arguments.');
end

%% Transform variables
if nargin == 3 && (isstruct(varargin{3}) || ischar(varargin{2}))

    Tol = sqrt(eps);    % Small number

    action = varargin{2};
    trinfo = varargin{3};

    if isempty(action)
        error('The transformation direction cannot be empty. Allowed values are direct (''dir'' or ''d'') and inverse (''inv'' or ''i'').');
    end

    if isempty(trinfo)
        % Empty TRINFO - consider as identity transformation

        x = varargin{1};

        switch lower(action(1))
            case {'d','i'}
                varargout{1} = x;
            case 'p'
                varargout{1} = ones(size(x,1),1);
            case 'l'
                varargout{1} = zeros(size(x,1),1);
            case {'m','f','g'}
                error('TRINFO is empty.');
            otherwise
                error(['Unkwnown transformation direction ''' action '''. Allowed values are direct (''dir'' or ''d'') and inverse (''inv'' or ''i'').']);
        end
    else

        scale = [];
        if isfield(trinfo,'scale') && ~isempty(trinfo.scale) && any(trinfo.scale ~= 1)
            scale = trinfo.scale;
        end

        if ~isfield(trinfo,'R_mat'); trinfo.R_mat = []; end

        nvars = numel(trinfo.lb_orig);  % Number of variables

        switch lower(action(1))
        %% DIRECT TRANSFORM
            case 'd'    % Direct transform
                x = varargin{1};
                y = x;
                a = trinfo.lb_orig;
                b = trinfo.ub_orig;
                mu = trinfo.mu;
                delta = trinfo.delta;

                % Unbounded scalars (possibly center and rescale)
                idx = trinfo.type == 0;
                if any(idx)
                    y(:,idx) = bsxfun(@rdivide,bsxfun(@minus,x(:,idx),mu(idx)),delta(idx));
                end

                % Lower bounded scalars
                idx = trinfo.type == 1;
                if any(idx)
                    y(:,idx) = log(bsxfun(@minus, x(:,idx), a(idx)));
                end

                % Upper bounded scalars
                idx = trinfo.type == 2;
                if any(idx)
                    y(:,idx) = log(bsxfun(@minus, b(idx), x(:,idx)));
                end

                % Lower and upper bounded scalars
                idx = trinfo.type == 3;
                if any(idx)
                    z = bsxfun(@rdivide, bsxfun(@minus, x(:,idx), a(idx)), ...
                        b(idx) - a(idx));
                    y(:,idx) = log(z./(1-z));
                    y(:,idx) = bsxfun(@rdivide,bsxfun(@minus,y(:,idx),mu(idx)),delta(idx));
                end

                % Lower and upper bounded scalars with Beta CDF transform
                idx = trinfo.type == 4;
                if any(idx)
                    alpha = trinfo.alpha;
                    beta = trinfo.beta;
                    for ii = find(idx)
                        % z = betacdf((x(:,ii) - a(ii)) / (b(ii) - a(ii)),alpha(ii),beta(ii));
                        z = min(max(eps,betacdf((x(:,ii) - a(ii)) / (b(ii) - a(ii)),alpha(ii),beta(ii))),1-eps);
                        y(:,ii) = log(z./(1-z));
                    end
                end

                % Lower and upper bounded scalars with Kumaraswamy CDF transform
                idx = trinfo.type == 5;
                if any(idx)
                    alpha = trinfo.alpha;
                    beta = trinfo.beta;
                    for ii = find(idx)
                        % z = kumarcdf((x(:,ii) - a(ii)) / (b(ii) - a(ii)),alpha(ii),beta(ii));
                        % p = min(max(eps,kumarcdf(z,alpha(ii),beta(ii))),1-eps);
                        % y(:,ii) = log(p./(1-p));
                        z = (x(:,ii) - a(ii)) / (b(ii) - a(ii));
                        % p = 1 - (1 - z.^alpha(ii)).^beta(ii);
                        y(:,ii) = log1p(-(1 - z.^alpha(ii)).^beta(ii)) - beta(ii)*log1p(-z.^alpha(ii));
                    end
                     y(:,idx) = bsxfun(@rdivide,bsxfun(@minus,y(:,idx),mu(idx)),delta(idx));
                end

                % Lower and upper bounded scalars with Kumaraswamy-logistic-power transform
                idx = trinfo.type == 6;
                if any(idx)
                    alpha = trinfo.alpha;
                    beta = trinfo.beta;
                    mu = trinfo.mu;
                    gamma = trinfo.gamma;
                    for ii = find(idx)
                        z = (x(:,ii) - a(ii)) / (b(ii) - a(ii));
                        y(:,ii) = log1p(-(1 - z.^alpha(ii)).^beta(ii)) - beta(ii)*log1p(-z.^alpha(ii)) - mu(ii);
                        y(:,ii) = sign(y(:,ii)).*abs(y(:,ii)).^gamma(ii);
                    end
                end

                % Lower and upper bounded scalars with nonparametric CDF transform
                idx = trinfo.type == 7;
                if any(idx)
                    xspace = trinfo.xspace;
                    pspace = trinfo.pspace;
                    for ii = find(idx)
                        y(:,ii) = norminv(interp1(xspace(ii,:),pspace(ii,:),x(:,ii)));
                    end
                end

                % Lower and upper bounded scalars with GMM CDF transform
                idx = trinfo.type == 8;
                if any(idx)
                    for ii = find(idx)
                        gmm = trinfo.gmm{ii};
                        z = (x(:,ii) - a(ii)) / (b(ii) - a(ii));
                        z = gmm.lambda*z + (1-gmm.lambda)*(gmm1cdf(z,gmm.w,gmm.Mu,gmm.Sigma)-gmm.Min)./gmm.Norm;
                        y(:,ii) = norminv(z);
                        % y(:,ii) = logiinv(z);
                    end
                end

                % Lower and upper bounded scalars with Kumaraswamy-logit transform
                idx = trinfo.type == 9;
                if any(idx)
                    alpha = trinfo.alpha;
                    beta = trinfo.beta;
                    for ii = find(idx)
                        z = (x(:,ii) - a(ii)) / (b(ii) - a(ii));
                        u = z.^(alpha(ii));

                        % Small u (close to zero)
                        uzero = u < sqrt(eps);
                        y(uzero,ii) = log(beta(ii)) + alpha(ii)*log(z(uzero)) + 0.5*u(uzero)*(beta(ii)+1) + 0.5*((beta(ii)-1)^2/4-beta(ii))*u(uzero).^2;

                        % Large u (close to 1)
                        uone = u > (1 - sqrt(eps));
                        y(uone,ii) = log(1./(1 - u(uone)).^beta(ii) - 1);

                        % Other values
                        y(~uzero & ~uone,ii) = log1p(-(1 - z(~uzero & ~uone).^alpha(ii)).^beta(ii)) - beta(ii)*log1p(-z(~uzero & ~uone).^alpha(ii));
                    end
                end

                % Unbounded with logistic-Kumaraswamy-logit transform
                idx = trinfo.type == 10;
                if any(idx)
                    alpha = trinfo.alpha;
                    beta = trinfo.beta;
                    for ii = find(idx)
                        z = (x(:,ii)-mu(:,ii)) ./ delta(:,ii);
                        %z = exp(z)./(exp(z)+1);
                        %y(:,ii) = log1p(-(1 - z.^alpha(:,ii)).^beta(:,ii)) - beta(:,ii).*log1p(-z.^alpha(:,ii));
                        s = 1./(1+exp(-z));
                        u = s.^alpha(ii);

                        % Small u (near zero)
                        uzero = u < Tol;
                        %y(uzero,ii) = log(beta(ii)) + alpha(ii)*log(s(uzero)) + beta(ii)*u(uzero);
                        y(uzero,ii) = log(beta(ii)) + alpha(ii)*z(uzero) + beta(ii)*u(uzero);

                        % Large u (near one)
                        uone = u > (1 - Tol);
                        w = 1./(1 + exp(z(uone)));
                        y(uone,ii) = -log(alpha(ii)) + beta(ii)*z(uone) - (alpha(ii)*w).^beta(ii);

                        % All other cases
                        urest = ~uzero & ~uone;
                        y(urest,ii) = log1p(-(1 - s(urest).^alpha(ii)).^beta(ii)) - beta(ii).*log1p(-s(urest).^alpha(ii));
                    end
                end

                % Lower and upper bounded scalars with Kumaraswamy-logit inverse CDF transform
                idx = trinfo.type == 11;
                if any(idx)
                    alpha = trinfo.alpha;
                    beta = trinfo.beta;
                    for ii = find(idx)
                        % z = kumarinv((x(:,ii) - a(ii)) / (b(ii) - a(ii)),alpha(ii),beta(ii));
                        % p = min(max(eps,kumarcdf(z,alpha(ii),beta(ii))),1-eps);
                        % y(:,ii) = log(p./(1-p));
                        z = (x(:,ii) - a(ii)) / (b(ii) - a(ii));
                        % p = 1 - (1 - z.^alpha(ii)).^beta(ii);
                        % y(:,ii) = 1./alpha(ii).*log1p(-(1 - z).^(1./beta(ii))) - log1p(-(1-(1-z).^(1./beta(ii))).^(1./alpha(ii)));
                        u = (1./beta(ii)).*log1p(-z);

                        uidx = u < log(eps)/2;
                        % u(uidx) = 1 - 1./alpha(ii).*exp(u(uidx));
                        y(uidx,ii) = -1./alpha(ii).*exp(u(uidx)) + log(alpha(ii)) - u(uidx);
                        % u(~uidx) = (1 - exp(u(~uidx))).^(1/alpha(ii));
                        y(~uidx,ii) = (1/alpha(ii)).*log1p(-exp(u(~uidx))) - log1p(-(1 - exp(u(~uidx))).^(1/alpha(ii)));
                        usmall = abs(u) < sqrt(eps);
                        y(usmall,ii) = (1/alpha(ii)).*log(-u(usmall)) - log1p(-u(usmall)/alpha(ii));

                        % y(:,ii) = log(expm1(-1./alpha(ii).*log1p(-(1-z).^(1./beta(ii)))));
                        y(:,ii) = bsxfun(@rdivide,bsxfun(@minus,y(:,ii),mu(ii)),delta(ii));
                    end
                end

                % Lower and upper bounded scalars (cumulative normal)
                idx = trinfo.type == 12;
                if any(idx)
                    z = bsxfun(@rdivide, bsxfun(@minus, x(:,idx), a(idx)), ...
                        b(idx) - a(idx));
                    y(:,idx) = -sqrt(2).*erfcinv(2*z);
                    y(:,idx) = bsxfun(@rdivide,bsxfun(@minus,y(:,idx),mu(idx)),delta(idx));
                end

                % Lower and upper bounded scalars (Student-t, nu = 4)
                idx = trinfo.type == 13;
                if any(idx)
                    z = bsxfun(@rdivide, bsxfun(@minus, x(:,idx), a(idx)), ...
                        b(idx) - a(idx));
                    aa = sqrt(4*z.*(1-z));
                    q = cos(acos(aa)/3)./aa;
                    y(:,idx) = sign(z - 0.5).*(2.*sqrt(q-1));
                    y(:,idx) = bsxfun(@rdivide,bsxfun(@minus,y(:,idx),mu(idx)),delta(idx));
                end

                % Rotate output
                if ~isempty(trinfo.R_mat); y = y*trinfo.R_mat; end

                % Rescale output
                if ~isempty(scale); y = bsxfun(@rdivide,y,scale); end

                varargout{1} = y;

            %% INVERSE TRANSFORM
            case 'i'    % Inverse transform
                y = varargin{1};
                % Rescale input
                if ~isempty(scale); y = bsxfun(@times,y,scale); end

                % Rotate input
                if ~isempty(trinfo.R_mat); y = y*trinfo.R_mat'; end

                x = y;
                a = trinfo.lb_orig;
                b = trinfo.ub_orig;
                mu = trinfo.mu;
                delta = trinfo.delta;

                % Unbounded scalars (possibly unscale and uncenter)
                idx = trinfo.type == 0;
                if any(idx)
                    x(:,idx) = bsxfun(@plus,bsxfun(@times,y(:,idx),delta(idx)),mu(idx));
                end

                % Lower bounded scalars
                idx = trinfo.type == 1;
                if any(idx)
                    x(:,idx) = bsxfun(@plus, exp(y(:,idx)), a(idx));
                end

                % Upper bounded scalars
                idx = trinfo.type == 2;
                if any(idx)
                    x(:,idx) = bsxfun(@minus, b(idx), exp(y(:,idx)));
                end

                % Lower and upper bounded scalars
                idx = trinfo.type == 3;
                if any(idx)
                    x(:,idx) = bsxfun(@plus,bsxfun(@times,y(:,idx),delta(idx)),mu(idx));
                    x(:,idx) = bsxfun(@plus, a(:,idx), bsxfun(@times, ...
                        b(idx)-a(idx), 1./(1+exp(-x(:,idx)))));
                end

                % Lower and upper bounded scalars with Beta CDF transform
                idx = trinfo.type == 4;
                if any(idx)
                    alpha = trinfo.alpha;
                    beta = trinfo.beta;
                    for ii = find(idx)
                        z = 1./(1+exp(-y(:,ii)));
                        x(:,ii) = a(ii) + (b(ii)-a(ii))*betainv(z,alpha(ii),beta(ii));
                    end
                end

                % Lower and upper bounded scalars with Kumaraswamy CDF transform
                idx = trinfo.type == 5;
                if any(idx)
                    alpha = trinfo.alpha;
                    beta = trinfo.beta;
                    y(:,idx) = bsxfun(@plus,bsxfun(@times,y(:,idx),delta(idx)),mu(idx));
                    for ii = find(idx)
                        % z = 1./(1+exp(-y(:,ii)));
                        % x(:,ii) = a(ii) + (b(ii)-a(ii))*kumarinv(z,alpha(ii),beta(ii));
                        z = exp(-y(:,ii))./(1+exp(-y(:,ii)));
                        x(:,ii) = a(ii) + (b(ii)-a(ii))*(1-z.^(1/beta(ii))).^(1/alpha(ii));
                    end
                end

                % Lower and upper bounded scalars with Kumaraswamy-logistic-power transform
                idx = trinfo.type == 6;
                if any(idx)
                    alpha = trinfo.alpha;
                    beta = trinfo.beta;
                    mu = trinfo.mu;
                    gamma = trinfo.gamma;
                    for ii = find(idx)
                        z = sign(y(:,ii)).*abs(y(:,ii)).^(1/gamma(ii)) + mu(ii);
                        z = exp(-z)./(1+exp(-z));
                        x(:,ii) = a(ii) + (b(ii)-a(ii))*(1-z.^(1/beta(ii))).^(1/alpha(ii));
                    end
                end

                % Lower and upper bounded scalars with nonparametric CDF transform
                idx = trinfo.type == 7;
                if any(idx)
                    xspace = trinfo.xspace;
                    pspace = trinfo.pspace;
                    for ii = find(idx)
                        x(:,ii) = interp1(pspace(ii,:),xspace(ii,:),normcdf(y(:,ii)));
                    end
                end

                % Lower and upper bounded scalars with GMM CDF transform
                idx = trinfo.type == 8;
                if any(idx)
                    for ii = find(idx)
                        gmm = trinfo.gmm{ii};
                        z = normcdf(y(:,ii));
                        % z = logicdf(y(:,ii));
                        for j = 1:size(z,1)
                            z(j) = tgmminv(z(j),gmm);
                        end
                        x(:,ii) = a(ii) + (b(ii)-a(ii)).*z;
                    end
                end

                % Lower and upper bounded scalars with Kumaraswamy-logit transform
                idx = trinfo.type == 9;
                if any(idx)
                    alpha = trinfo.alpha;
                    beta = trinfo.beta;
                    for ii = find(idx)
                        z = exp(-y(:,ii))./(1+exp(-y(:,ii)));   % 1 - logistic(z)
                        x(:,ii) = a(ii) + (b(ii)-a(ii))*(1-z.^(1/beta(ii))).^(1/alpha(ii));
                    end
                end

                % Unbounded scalars with logistic-Kumaraswamy-logit transform
                idx = trinfo.type == 10;
                if any(idx)
                    alpha = trinfo.alpha;
                    beta = trinfo.beta;
                    for ii = find(idx)
                        % z = exp(-y(:,ii))./(1+exp(-y(:,ii)));   % 1 - logistic(z)
                        % % z = (1-z.^(1/beta(ii))).^(1/alpha(ii));
                        % z = 1./alpha(:,ii).*log1p(-z.^(1./beta(:,ii))) - log1p(-(1-z.^(1./beta(:,ii))).^(1./alpha(:,ii)));

                        % Small u (~ zero)
                        uzero = y(:,ii) < log(Tol);
                        u = 1./(1 + exp(-y(uzero,ii)));
                        %x(uzero,ii) = 1/alpha(ii)*(-log1p(exp(-y(uzero,ii)))-log(beta(ii))) - log1p(-(u/beta(ii)).^(1/alpha(ii)));
                        x(uzero,ii) = 1/alpha(ii)*(y(uzero,ii)-log(beta(ii))) - log1p(-(u/beta(ii)).^(1/alpha(ii)));

                        % Large u (~ one)
                        uone = y(:,ii) > -log(Tol);
                        w = 1./(1 + exp(y(uone,ii))).^(1/beta(ii));
                        x(uone,ii) = 1/alpha(ii)*log1p(-w) + 1/beta(ii)*y(uone,ii) + log(alpha(ii));

                        % All other cases
                        urest = ~uzero & ~uone;
                        z = 1./(1+exp(y(urest,ii)));
                        x(urest,ii) = 1./alpha(:,ii).*log1p(-z.^(1./beta(:,ii))) - log1p(-(1-z.^(1./beta(:,ii))).^(1./alpha(:,ii)));

                        x(:,ii) = x(:,ii).*delta(:,ii) + mu(:,ii);
                    end
                end

                % Lower and upper bounded scalars with Kumaraswamy inverse CDF transform
                idx = trinfo.type == 11;
                if any(idx)
                    alpha = trinfo.alpha;
                    beta = trinfo.beta;
                    y(:,idx) = bsxfun(@plus,bsxfun(@times,y(:,idx),delta(idx)),mu(idx));
                    for ii = find(idx)
                        z = 1./(1+exp(-y(:,ii)));
                        x(:,ii) = a(ii) + (b(ii)-a(ii))*(1 - (1-z.^(alpha(ii))).^beta(ii));
                    end
                end

                % Lower and upper bounded scalars (cumulative normal)
                idx = trinfo.type == 12;
                if any(idx)
                    x(:,idx) = bsxfun(@plus,bsxfun(@times,y(:,idx),delta(idx)),mu(idx));
                    x(:,idx) = bsxfun(@plus, a(:,idx), bsxfun(@times, ...
                        b(idx)-a(idx), 0.5 * erfc(-x(:,idx) ./ sqrt(2))));
                end

                % Lower and upper bounded scalars (Student-t, nu = 4)
                idx = trinfo.type == 13;
                if any(idx)
                    x(:,idx) = bsxfun(@plus,bsxfun(@times,y(:,idx),delta(idx)),mu(idx));
                    t2 = x(:,idx).^2;
                    f = 0.5 + 3/8*x(:,idx)./sqrt(1 + t2/4).*(1 - t2./(1 + t2/4)/12);
                    x(:,idx) = bsxfun(@plus, a(:,idx), bsxfun(@times, ...
                        b(idx)-a(idx), f));
                end

                % Force to stay within bounds
                a(isfinite(a)) = a(isfinite(a)) + eps(a(isfinite(a)));
                b(isfinite(b)) = b(isfinite(b)) - eps(b(isfinite(b)));
                x = bsxfun(@min,bsxfun(@max,x,a),b);
                varargout{1} = x;

            %% PDF (OR LOG PDF) CORRECTION
            case {'p','l','g'}  % pdf (or log pdf) correction
                y = varargin{1};
                % Rescale input
                if ~isempty(scale); y = bsxfun(@times,y,scale); end

                % Rotate input
                if ~isempty(trinfo.R_mat); y = y*trinfo.R_mat'; end

                logpdf_flag = strcmpi(action(1),'l');
                if logpdf_flag
                    p = zeros(size(y,1),nvars);
                else
                    p = ones(size(y,1),nvars);
                end
                grad_flag = strcmpi(action(1),'g');

                a = trinfo.lb_orig;
                b = trinfo.ub_orig;
                mu = trinfo.mu;
                delta = trinfo.delta;

                % Unbounded scalars
                idx = trinfo.type == 0;
                if any(idx)
                    p(:,idx) = repmat(log(delta(idx)),[size(p,1),1]);
                end

                % Lower or upper bounded scalars
                idx = trinfo.type == 1 | trinfo.type == 2;
                if any(idx)
                    p(:,idx) = y(:,idx);
                end

                % Lower and upper bounded scalars
                idx = trinfo.type == 3;
                if any(idx)
                    y(:,idx) = bsxfun(@plus,bsxfun(@times,y(:,idx),delta(idx)),mu(idx));
                    z = -log1p(exp(-y(:,idx)));
                    p(:,idx) = bsxfun(@plus, log(b(idx)-a(idx)), -y(:,idx) + 2*z);
                    p(:,idx) = bsxfun(@plus, p(:,idx), log(delta(idx)));
                end

                % Lower and upper bounded scalars with Beta CDF transform
                idx = trinfo.type == 4;
                if any(idx)
                    for ii = find(idx)
                        alpha = trinfo.alpha;
                        beta = trinfo.beta;

                        z = -log1p(exp(-y(:,ii)));
                        x = min(max(eps,betainv(1./(1+exp(-y(:,ii))),alpha(ii),beta(ii))),1-eps);
                        logbeta = (alpha(ii)-1)*log(x) + (beta(ii)-1)*log1p(-x) ...
                            + gammaln(alpha(ii)+beta(ii)) - gammaln(alpha(ii)) - gammaln(beta(ii));
                        p(:,ii) = log(b(ii)-a(ii)) -logbeta -y(:,ii) + 2*z;

                        if any(~isfinite(p))
                            fprintf('aaaa!');
                        end
                    end
                end

                % Lower and upper bounded scalars with Kumaraswamy CDF transform
                idx = trinfo.type == 5;
                if any(idx)
                    for ii = find(idx)
                        alpha = trinfo.alpha;
                        beta = trinfo.beta;

                        y(:,ii) = bsxfun(@plus,bsxfun(@times,y(:,ii),delta(ii)),mu(ii));

                        z = -log1p(exp(-y(:,ii)));
                        % x = kumarinv(1./(1+exp(-y(:,ii))),alpha(ii),beta(ii));

                        %p = 1./(1+exp(-y(:,ii)));
                        %x = (-((-(p-1)).^(1/beta(ii))-1)).^(1/alpha(ii));
                        %log(-(p-1)) = -y(,::) - log1p(exp(-y(:,ii))));

                        %(1/beta(ii)).*log(-(p-1))

                        % x = min(max(eps,kumarinv(1./(1+exp(-y(:,ii))),alpha(ii),beta(ii))),1-eps);
                        u = 1./(1 + exp(-y(:,ii)));
                        logf = (1-1/alpha(ii))*log1p(-(1-u).^(1/beta(ii))) + (1-1/beta(ii))*(-y(:,ii)+z) ...
                            + log(alpha(ii)) + log(beta(ii));

                        % Special case for very small u
                        %idx_small = u < 1e12;
                        %logf(idx_small) = (1-1/alpha(ii))*(-y(idx_small,ii) + z - log(beta(ii))) + (1-1/beta(ii))*(-y(idx_small,ii)+z(idx_small)) ...
                        %    + log(alpha(ii)) + log(beta(ii));

                        %logf = (alpha(ii)-1)*log(x) + (beta(ii)-1)*log1p(-x.^alpha(ii)) ...
                        %    + log(alpha(ii)) + log(beta(ii));
                        p(:,ii) = log(b(ii)-a(ii)) -logf -y(:,ii) + 2*z;
                        p(:,ii) = bsxfun(@plus, p(:,ii), log(delta(ii)));

                        if any(~isfinite(p))
                            p(~isfinite(p)) = -Inf;
                            fprintf('aaaa!');
                        end
                    end
                end

                % Lower and upper bounded scalars with Kumaraswamy-logistic-power transform
                idx = trinfo.type == 6;
                if any(idx)
                    for ii = find(idx)
                        alpha = trinfo.alpha;
                        beta = trinfo.beta;
                        mu = trinfo.mu;
                        gamma = trinfo.gamma;

                        yl = sign(y(:,ii)).*abs(y(:,ii)).^(1/gamma(ii)) + mu(ii);

                        z = -log1p(exp(-yl));
                        u = 1./(1 + exp(-yl));
                        logf = (1-1/alpha(ii))*log1p(-(1-u).^(1/beta(ii))) + (1-1/beta(ii))*(-yl+z) ...
                            + log(alpha(ii)) + log(beta(ii));
                        logf = logf + log(gamma(ii)) + (gamma(ii)-1)*log(abs(yl-mu(ii)));
                        p(:,ii) = log(b(ii)-a(ii)) -logf -yl + 2*z;
                    end
                end

                % Lower and upper bounded scalars with nonparametric CDF transform
                idx = trinfo.type == 7;
                if any(idx)
                    xspace = trinfo.xspace;
                    pspace = trinfo.pspace;
                    for ii = find(idx)
                        z = -0.5*log(2*pi) -0.5*y(:,ii).^2;

                        yinv = normcdf(y(:,ii));
                        [~,pos] = histc(yinv,pspace(ii,:));
                        dx = [-log(diff(xspace(ii,:))),Inf];
                        logf = dx(pos);

                        p(:,ii) = log(b(ii)-a(ii)) -logf(:) + z;
                    end
                end

                % Lower and upper bounded scalars with GMM CDF transform
                idx = trinfo.type == 8;
                if any(idx)
                    for ii = find(idx)
                        gmm = trinfo.gmm{ii};

                        z = -0.5*log(2*pi) -0.5*y(:,ii).^2;
                        % z = -y(:,ii) - 2*log1p(exp(-y(:,ii)));

                        yinv = normcdf(y(:,ii));
                        for j = 1:size(yinv,1); yinv(j) = tgmminv(yinv(j),gmm); end

                        logf = log(gmm.lambda + (1-gmm.lambda)*gmm1pdf(yinv,gmm.w,gmm.Mu,gmm.Sigma)./gmm.Norm);

                        p(ii) = log(b(ii)-a(ii)) - logf + z;

                        if any(~isfinite(p))
                            p(~isfinite(p)) = -Inf;
                            fprintf('aaaa!');
                        end
                    end
                end

                % Lower and upper bounded scalars with Kumaraswamy-logit transform
                idx = trinfo.type == 9;
                if any(idx)
                    for ii = find(idx)
                        alpha = trinfo.alpha;
                        beta = trinfo.beta;
                        nf = (b(ii)-a(ii))/alpha(ii)/beta(ii);
                        k = 1./(1+exp(-y(:,ii)));

                        logk = -log1p(exp(-y(:,ii)));
                        log1mk = logk - y(:,ii);
                        logz = 1/alpha(ii)*log1p(-(1-k).^(1/beta(ii)));

                        % k small (close to zero)
                        kzero = k < sqrt(eps);
                        log1mk(kzero) = log1p(-k(kzero));
                        logz(kzero) = 1/alpha(ii) * (logk(kzero)-log(beta(ii))) + 1/alpha(ii)*log1p(0.5*(1-1/beta(ii)*k(kzero)));

                        %z = (1-(1-k).^(1/beta(ii))).^(1/alpha(ii));
                        % 1 - z^alpha = (1-k).^(1/beta(ii))
                        p(:,ii) = log(nf) + (1/beta(ii)-1) * log1mk + (1-alpha(ii))*logz -y(:,ii)+2*logk;
                    end
                end

                % Unbounded scalars with logistic-Kumaraswamy-logit transform
                idx = trinfo.type == 10;
                if any(idx)
                    for ii = find(idx)
                        alpha = trinfo.alpha;
                        beta = trinfo.beta;

                        lnf = log(delta(:,ii)./alpha(:,ii)./beta(:,ii));

                        if 0

                            u = 1./(1+exp(y(:,ii)));
                            logu = -log1p(exp(y(:,ii)));
                            z = (1-u.^(1./beta(:,ii))).^(1./alpha(:,ii));
                            p(:,ii) = lnf + y(:,ii) + (1+1./beta(:,ii)).*logu - log1p(-u.^(1./beta(:,ii))) - log1p(-z);

                        else

                            p(:,ii) = lnf -y(:,ii) - 2*log1p(exp(-y(:,ii)));

                            % Small u (~ zero)
                            uzero = y(:,ii) < log(Tol);
                            u = 1./(1 + exp(-y(uzero,ii)));
                            % p(uzero,ii) = p(uzero,ii) - log(u/beta(ii)) + (1/beta(ii)-1)*log1p(-u);
                            p(uzero,ii) = p(uzero,ii) - y(uzero,ii) + log(beta(ii)) + (1/beta(ii)-1)*log1p(-u);
                            t = (1 - (1-u).^(1/beta(ii))).^(1/alpha(ii));
                            p(uzero,ii) = p(uzero,ii) - log1p(-t);

                            % Large u (~ one)
                            uone = y(:,ii) > -log(Tol);
                            w = 1./(1 + exp(y(uone,ii)));
                            p(uone,ii) = p(uone,ii) - log1p(-w.^(1/beta(ii))) + y(uone,ii) + log(alpha(ii));

                            % All other cases
                            urest = ~uzero & ~uone;
                            u = 1./(1 + exp(-y(urest,ii)));
                            p(urest,ii) = p(urest,ii) - log1p(-(1-u).^(1/beta(ii))) + (1/beta(ii)-1)*log1p(-u);
                            t = (1 - (1-u).^(1/beta(ii))).^(1/alpha(ii));
                            p(urest,ii) = p(urest,ii) - log1p(-t);

                        end

                    end
                end

                % Lower and upper bounded scalars with inverse Kumaraswamy logit CDF transform
                idx = trinfo.type == 11;
                if any(idx)
                    for ii = find(idx)
                        alpha = trinfo.alpha;
                        beta = trinfo.beta;

                        y(:,ii) = bsxfun(@plus,bsxfun(@times,y(:,ii),delta(ii)),mu(ii));

                        y_large = y(:,ii) > -log(eps)/2;

                        z = 1./(1+exp(-y(:,ii)));
                        u = 1 - (1-z.^(alpha(ii))).^beta(ii);

                        u(y_large) = 1 - alpha(ii)*exp(-beta(ii)*y(y_large,ii))./(1+alpha(ii)*beta(ii)*exp(-y(y_large,ii)));

%                        z = -log1p(exp(-y(:,ii)));
                        %u = 1./(1 + exp(-y(:,ii)));
                        %u = a(ii) + (b(ii)-a(ii))*(1 - (1-u.^(alpha(ii))).^beta(ii));
%                        logf = log(alpha(ii)) + log(beta(ii)) + (alpha(ii)-1).*log(u) + (beta(ii)-1) .* log1p(-u.^alpha(ii));
                         logf = (1-1/alpha(ii))*log1p(-(1-u).^(1/beta(ii))) + (1-1/beta(ii))*log1p(-u) ...
                             + log(alpha(ii)) + log(beta(ii));
                         logf(y_large) = (1-1/alpha(ii))*log1p(-alpha(ii)*exp(-y(y_large,ii))./(1+alpha(ii)*exp(-y(y_large,ii)))) ...
                             + (1-1/beta(ii))*log(alpha(ii)*exp(-beta(ii)*y(y_large,ii))./(1+alpha(ii)*beta(ii)*exp(-y(y_large,ii)))) ...
                             + log(alpha(ii)) + log(beta(ii));

%                         logf = (1-1/alpha(ii))*log1p(-(1-u).^(1/beta(ii))) + (1-1/beta(ii))*(-y(:,ii)+z) ...
%                             + log(alpha(ii)) + log(beta(ii));

                        w = log1p(-z);
                        w(y_large) = -y(y_large,ii) - log1p(exp(-y(y_large,ii)));

                        p(:,ii) = log(b(ii)-a(ii)) + logf + log(z) + w; %-y(:,ii) + 2*z;
                        p(:,ii) = bsxfun(@plus, p(:,ii), log(delta(ii)));

                        if any(~isfinite(p))
                            p(~isfinite(p)) = -Inf;
                            fprintf('aaaa!');
                        end
                    end
                end

                % Lower and upper bounded scalars (cumulative normal)
                idx = trinfo.type == 12;
                if any(idx)
                    y(:,idx) = bsxfun(@plus,bsxfun(@times,y(:,idx),delta(idx)),mu(idx));
                    z = -0.5*log(2*pi) - 0.5*y(:,idx).^2;
                    p(:,idx) = bsxfun(@plus, log(b(idx)-a(idx)), z);
                    p(:,idx) = bsxfun(@plus, p(:,idx), log(delta(idx)));
                end

                % Lower and upper bounded scalars (Student-t, nu = 4)
                idx = trinfo.type == 13;
                if any(idx)
                    y(:,idx) = bsxfun(@plus,bsxfun(@times,y(:,idx),delta(idx)),mu(idx));
                    z = log(3/8) - 5/2*log1p(y(:,idx).^2/4);
                    p(:,idx) = bsxfun(@plus, log(b(idx)-a(idx)), z);
                    p(:,idx) = bsxfun(@plus, p(:,idx), log(delta(idx)));
                end

                %if ~isempty(trinfo.R_mat) && lower(action(1)) == 'g'
                %    p = p*(trinfo.R_mat*diag();
                %end

                % Scale transform
                if ~isempty(scale) && ~grad_flag
                    p = bsxfun(@plus,p,log(scale));
                end

                if ~grad_flag; p = sum(p,2); end
                if ~logpdf_flag; p = exp(p); end

                varargout{1} = p;

            %% FIRST DERIVATIVE WRT TRANSFORMATION PARAMETERS (ignores final rotation and scaling)
            case {'f'}
                y = varargin{1};
                nvars = numel(trinfo.lb_orig);

                % Rescale input
                if ~isempty(scale); y = bsxfun(@times,y,scale); end

                % Rotate input
                if ~isempty(trinfo.R_mat); y = y*trinfo.R_mat'; end

                a = trinfo.lb_orig;
                b = trinfo.ub_orig;
                delta = trinfo.delta;

                % Lower and upper bounded scalars with Kumaraswamy-logit transform
                % and unbounded scalars with logistic-Kumaraswamy-logit transform
                idx = (trinfo.type == 9 | trinfo.type == 10);
                if any(idx)
                    for ii = find(idx)
                        alpha = trinfo.alpha;
                        beta = trinfo.beta;

                        k = 1./(1+exp(-y(:,ii)));
                        talpha = 1 - (1-k).^(1/beta(ii));
                        logt = 1/alpha(ii) .* log1p(-(1-k).^(1/beta(ii)));

                        p(:,ii) = talpha.*beta(ii).*logt./(1-talpha)./k;
                        p(:,ii+nvars) = -log1p(-talpha)./k;
                    end
                end

                varargout{1} = p;

            %% MIXED DERIVATIVE WRT TRANSFORMATION PARAMETERS of FIRST DERIVATIVE
            % (ignores final rotation and scaling)
            case {'m'}
                y = varargin{1};
                nvars = numel(trinfo.lb_orig);

                % Rescale input
                if ~isempty(scale); y = bsxfun(@times,y,scale); end

                % Rotate input
                if ~isempty(trinfo.R_mat); y = y*trinfo.R_mat'; end

                a = trinfo.lb_orig;
                b = trinfo.ub_orig;
                delta = trinfo.delta;

                % Lower and upper bounded scalars with Kumaraswamy-logit transform
                % and unbounded scalars with logistic-Kumaraswamy-logit transform
                idx = (trinfo.type == 9 | trinfo.type == 10);
                if any(idx)
                    for ii = find(idx)
                        alpha = trinfo.alpha;
                        beta = trinfo.beta;

                        k = 1./(1+exp(-y(:,ii)));
                        talpha = 1 - (1-k).^(1/beta(ii));
                        logt = 1/alpha(ii) .* log1p(-(1-k).^(1/beta(ii)));
                        t = talpha.^(1/alpha(ii));

                        if trinfo.type(ii) == 9
                            nf = - 1./(b(ii) - a(ii)) ./ t;
                        else
                            nf = (t-1)./delta(ii);
                        end
                        den = 1 ./ k.^2 ./ (talpha-1) .* nf;
                        p(:,ii) = talpha.*beta(ii).*((-k.*(1+alpha(ii).*logt)) + talpha.*(1 + (1-k).*(-1+alpha(ii).*beta(ii).*logt))) ...
                            .* den ./ (talpha - 1);
                        p(:,ii+nvars) = talpha .* alpha(ii) .* (1 + (1-k).*(-1 - y(:,ii) - log1p(exp(-y(:,ii))))) ...
                            .* den;
                    end
                end

                varargout{1} = p;

            otherwise
                error(['Unkwnown transformation direction ''' action '''. Allowed values are direct (''dir'' or ''d'') and inverse (''inv'' or ''i'').']);
        end
    end

else
%% Create transform

    nvars = varargin{1};
    lb = varargin{2}(:)';
    ub = varargin{3}(:)';
    if nargin > 3
        plb = varargin{4}(:)';
        pub = varargin{5}(:)';
    else
        plb = []; pub = [];
    end
    if nargin > 5
        bounded_type = varargin{6};
    else
        % Default bounded type is logit
        bounded_type = 3;
    end

    % Empty LB and UB are Infs
    if isempty(lb); lb = -Inf; end
    if isempty(ub); ub = Inf; end

    % Empty plausible bounds equal hard bounds
    if isempty(plb); plb = lb; end
    if isempty(pub); pub = ub; end

    % Convert scalar inputs to row vectors
    if isscalar(lb); lb = lb*ones(1,nvars); end
    if isscalar(ub); ub = ub*ones(1,nvars); end
    if isscalar(plb); plb = plb*ones(1,nvars); end
    if isscalar(pub); pub = pub*ones(1,nvars); end

    % Check that the order of bounds is respected
    assert(all(lb <= plb & plb < pub & pub <= ub), ...
        'Variable bounds should be LB <= PLB < PUB <= UB for all variables.');

    % Transform to log coordinates
    trinfo.lb_orig = lb;
    trinfo.ub_orig = ub;

    trinfo.type = zeros(1,nvars);
    for i = 1:nvars
        if isfinite(lb(i)) && isinf(ub(i)); trinfo.type(i) = 1; end
        if isinf(lb(i)) && isfinite(ub(i)); trinfo.type(i) = 2; end
        if isfinite(lb(i)) && isfinite(ub(i)) && lb(i) < ub(i); trinfo.type(i) = bounded_type; end
    end

    % Centering (at the end of the transform)
    trinfo.mu = zeros(1,nvars);
    trinfo.delta = ones(1,nvars);

    % Get transformed PLB and PUB
    plb = warpvars_vbmc(plb,'d',trinfo);
    pub = warpvars_vbmc(pub,'d',trinfo);

    % Center in transformed space
    for i = 1:nvars
        if isfinite(plb(i)) && isfinite(pub(i))
            trinfo.mu(i) = 0.5*(plb(i)+pub(i));
            trinfo.delta(i) = (pub(i)-plb(i));
        end
    end

    varargout{1} = trinfo;

end

end

%--------------------------------------------------------------------------
function y = tgmminv(p,gmm)
    if p <= 0; p = eps; elseif p >= 1; p = 1 - eps; end    % Correct bounds
    fun = @(z) tgmmfzero(z,p,gmm.w,gmm.Mu,gmm.Sigma,gmm.Min,gmm.iMax,gmm.Norm,gmm.lambda);
    y = qfzero(fun,[-eps,1+eps]);
    if y <= 0; y = eps; elseif y >= 1; y = 1 - eps; end    % Correct bounds
end


text:
function a = reshape(b,varargin)
% Overloaded reshape function for file_array objects
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

%
% $Id: reshape.m 7147 2017-08-03 14:07:01Z spm $

if length(struct(b))~=1, error('Can only reshape simple file_array objects.'); end

args = [];
for i=1:length(varargin)
    args = [args varargin{i}(:)'];
end
if prod(args)~=prod(b.dim)
    error('To RESHAPE the number of elements must not change.');
end
a = b;
a.dim = args;

author: spm


text:
function srtlist = fuf (folders,reclev,opt)
%FUF recursively retrieves files under specified folder(s).
%
%   Syntax: srtlist = fuf (folders,reclev,opt)
%
%   INPUT ARGUMENTS
%       - folders       can be a character array or a cell array of strings. Valid strings
%                       are matlab path;
%       - reclev        {1} optional scalar/numeric array  of  booleans: 1 stands for recursive
%                       search, whilst 0 limit the search to the specified folder.
%                       reclev can be either a scalar, or a 1 X length(folders) array.
%                       Use the 2nd option to specify the search rule on a folder-basis.
%                       Default is 1, i.d., each folder is searched for recursively.
%       - opt           {'normal'} | 'detail'  is an optional field: when it's set to 'detail',
%                       FUF returns the full path of files searched for, while, when set to
%                       normal, it returns file names only.
%
%
%   OUTPUT ARGUMENT
%       - srtlist   is a sorted cell array of strings containing the name or the full path
%                       of files recursively found under the given folders.
%
%   FUF applies the wildcard * for searching files through all the directories beneath the given folder(s)
%   if the reclev parameter is set to 1. If you want to have different wildcards for folders belonging
%   to the same tree root, you need to pass them separately as input arguments.
%
%   REMARKS
%       - Folders don't have to belong to the same directory, therefore this function doesn't
%         have to be invoked from a particular directory.
%       - Matlab partial paths are not valid input arguments.
%       - Wildcard * can be used to narrow the search.
%       - To get the full path of a given file, FUF doesn't use the which command,
%         therefore it normally works also with java function that are not loaded.
%       - After completion, the working directory is set to the current directory at the time
%         of the function call.
%
%   EXAMPLES
%
%   To retrieve all files under the folder Utils:
%
%       >>  dir('C:\matlabR12\work\Interface\Utils')
%
%           .                       Contents.m
%           ..                      my_interface_evalmcw.m
%
%   To retrieve recursively the full path of all .m files starting with Cont under the folder work:
%
%       >> fuf('C:\matlabR12\work\Con*.m','detail')
%
%           'C:\matlabR12\work\Laboratorio_modellistica\DTT\Contents.m'
%           'C:\matlabR12\work\Laboratorio_modellistica\SDF\Contents.m'
%           'C:\matlabR12\work\PREPOSTGUIS\Contents.m'
%           'C:\matlabR12\work\PREPOSTGUIS\Utils\Contents.m'
%           'C:\matlabR12\work\Pavia\Contents.m'
%           'C:\matlabR12\work\TwoLe_front_end\Contents.m'
%           [1x74 char]
%           'C:\matlabR12\work\TwoLe_front_end\My_Classes\@Sector\Contents.m'
%           'C:\matlabR12\work\TwoLe_front_end\Utils\Contents.m'
%           'C:\matlabR12\work\Utils\Contents.m'
%
%   To retrieve recursively under the folder PREPOSTGUIS all .fig files and non recursively, under
%   the folder MatlabR12 all .txt files:
%
%       >> fuf({'C:\matlabR12\work\PREPOSTGUIS\*.fig','C:\matlabR12\*.txt'},[1,0],'detail')
%
%               'C:\matlabR12\license.txt'
%               'C:\matlabR12\work\PREPOSTGUIS\Private\MV_Manager.fig'
%               'C:\matlabR12\work\PREPOSTGUIS\R_h.fig'
%               'C:\matlabR12\work\PREPOSTGUIS\TSA.fig'
%               'C:\matlabR12\work\PREPOSTGUIS\Visual_3D.fig'

%                                         -$-$-$-
%
%        Author:    Francesco di Pierro                  Reasearch Assistant
%                                                    Center for Water Systems (CWS)
%                                                Dep. of Engineering and Computer Science
%                                                        University of Exeter
%                                                    e-mail: <F.Di-Pierro@ex.ac.uk>
%
%                                         -$-$-$-

%--------------------CHECK INPUT ARGUMENT TYPE AND SET DEFAULT VALUES------------------------%

%argument cheking and parsing default parameters

error(nargchk(1,3,nargin));
error(nargoutchk(0,1,nargout));
INITIAL_DIR = pwd;
if nargin==1
    if ischar(folders)
        folders = cellstr(folders);
    elseif ~iscellstr(folders)
        error('The 1st argument to the function fuf must be either a string or a cellstring!')
    end
    reclev = ones(size(folders));
    opt = 'normal';
elseif nargin==2
    if ischar(folders)
        folders = cellstr(folders);
    elseif ~iscellstr(folders)
        error('The 1st argument to the function fuf must be either a string or a cellstring!')
    end
    if isnumeric(reclev)
        if ~all(ismember(reclev,[0,1]))
            error('The 2nd argument to the function  must be a scalar/vector of booleans!')
        end
        if prod(size(reclev))==1
            reclev = reclev*ones(size(folders));
        elseif prod(size(reclev))~=prod(size(folders))
            error('The 2nd argument to the function must be a either a scalar or a 1 X length(1st argument) array!')
        end
        opt = 'normal';
    elseif ~ischar(reclev) | ~any(strcmp(reclev,{'normal','detail'}))
        error('Mismatched argument specification; the 3rd argument can be either set to ''normal'' or ''detail''');
    else
        opt = reclev;
        reclev = ones(size(folders));
    end
else
    if ischar(folders)
        folders = cellstr(folders);
    elseif ~iscellstr(folders)
        error('The 1st argument to the function fuf must be either a string or a cellstring!')
    end
    if isnumeric(reclev)
        if ~all(ismember(reclev,[0,1]))
            error('The 2nd argument to the function must be a scalar/vector of booleans!')
        end
        if prod(size(reclev))==1
            reclev = reclev*ones(size(folders));
        elseif prod(size(reclev))~=prod(size(folders))
            error('The 2nd argument to the function must be a either a scalar or a 1 X length(1st argument) array!')
        end
    else
        error('The 2nd argument to the function must be a scalar/vector of booleans!')
    end
    if ~ischar(opt) | ~any(strcmp(opt,{'normal','detail'}))
        error('Mismatched argument specification: Tte 3rd argument to the function must be either set to ''normal'' or ''detail''');
    end
end

%scan folders searching for incorrect folder names and partialpaths!!!
EXIT = 0;
for i=1:length(folders)
    [d,f,e] = fileparts(folders{i});
    if (~isdir(d) & ~isempty(e))
        warning(['"',d,'" is not a valid folder name!'])
        EXIT = 1;
    elseif (~isdir(fullfile(d,f,'')) & isempty(e))
        warning(['"',fullfile(d,f,''),'" is not a valid folder name!'])
        EXIT = 1;
    elseif (isempty(dir(d)) & ~isempty(e)) | (isempty(dir(fullfile(d,f,''))) & isempty(e))
        warning('Matlab PARTIALPATHS not allowed!')
        EXIT = 1;
    else                                                %this cheks the very unlike event arising
        sub = dir(pwd);                                 %when the partial path provided is a directory
        [subel{1:length(sub)}] = deal(sub.name);        %under the current one: this, in fact, is the
        [subtype{1:length(sub)}] = deal(sub.isdir);     %only situation where the dir command handles
        testel1 = strcmp(subel,fullfile(d,f,''));       %partialpaths!!!
        testel2 = strcmp(subel,d);
        if ~isempty(subtype(testel1)) | ~isempty(subtype(testel2))
            warning('Matlab PARTIALPATHS not allowed!')
            EXIT = 1;
        end
    end
end
if EXIT
    disp('')
    error(strvcat('One or more not valid folder names encountered! Function aborted!'))
end

%--------------------------------CORE FUNCTION--------------------------------%

sorted_list = [];                                           %initialize the output list

sorted_list = rec(sorted_list,folders,reclev,opt);          %call the function

sorted_list = sortrows(sorted_list);                        %sort the list

if nargout, srtlist = sorted_list; else disp (' '),disp(sorted_list), end

cd(INITIAL_DIR);                                            %and set current directory back to the initial one

%--------------------------------RECURSIVE FUNCTION--------------------------------%

function sorted_list = rec(sorted_list,folders,reclev,opt)

for i=1:length(folders)
    %first build the new search condition made of all files satistying the
    %search condition and all directories under the current one
    val = []; val1 = []; val2 = []; cnval1 = {}; cnval2 = {}; cdval = {}; cdva2 = {};   %initialize helper variables
    [pth,fname,ext] = fileparts(folders{i});
    cd(pth);                                            %move to the directory: isdir only recognizes directories on the Matlab search path or the current one!
    if isdir(fname)
        wild = '';
        cd(fname)
        pth = fullfile(pth,fname);
        val = dir;
    else
        %get the filenames satisfying the search condition
        val1 = dir(folders{i});
        wild = [fname,ext];
        if isempty(val1)
            cnval1 = [];    cdval1 = [];
        else
            [cnval1{1:length(val1)}] = deal(val1.name); [cdval1{1:length(val1)}] = deal(val1.isdir);
        end
        %and the directory/ies under the current one
        val2 = dir;
        wo = logical(zeros(size(val2)));
        for k=1:length(val2)
            if val2(k).isdir
                wo(k) = 1;
            end
        end
        val2(~wo) = [];
        %and build the new search structure
        [cnval2{1:length(val2)}] = deal(val2.name); [cdval2{1:length(val2)}] = deal(val2.isdir);
        cnval = [cnval1,cnval2];    cdval = [cdval1,cdval2];
        [val(1:length(cnval)).name] = deal(cnval{:}); [val(1:length(cdval)).isdir] = deal(cdval{:});
    end
    for j=1:length(val)
        if (val(j).isdir) & not(strcmp(val(j).name,'.')) & not(strcmp(val(j).name,'..'))    %if the jth object under the ith folder is a valid folder name(directory)...
            new_path = fullfile(pth,val(j).name);                                           %set the current directory to that one:
            cd(new_path);
            if reclev(i)                            %recursively call the function keeping in mind the search condition
                sorted_list = rec(sorted_list,{fullfile(new_path,wild)},reclev(i),opt);
            end
            cd ..                                   %get back to the previous directory
        elseif not(strcmp(val(j).name,'.')) & not(strcmp(val(j).name,'..')) %if the jth object under the ith folder is a valid file name
            if strcmp(opt,'detail')
                sorted_list =[sorted_list; cellstr(fullfile(pwd,val(j).name))];
            else
                sorted_list =[sorted_list; cellstr(val(j).name)];
            end
        end
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function [mbar] = psi2mbar(psi)
% Convert units of pressure from pounds per square inch to millibar.
% Chad Greene 2012
mbar = psi*68.9476;


text:
clear
close all
clc;

addpath(genpath('../libs'))
path_to_matconvnet = '/home/skong2/scratch/matconvnet-1.0-beta23_modifiedDagnn';

run(fullfile(path_to_matconvnet, 'matlab', 'vl_setupnn'));
addpath(genpath(fullfile('dependencies', 'matconvnet','examples')));
%% read matconvnet model
load('imdb_toydata_v3_from_mnist.mat');
imdb.path = './toydata_v3';
imdb.path_to_dataset = './toydata_v3';
% set GPU
gpuId = 4;
gpuDevice(gpuId);
flagSaveFig = true; % {true false} whether to store the result

saveFolder = 'main007_instSeg_v1_absEucMM';
modelName = 'softmax_net-epoch-83.mat';

netbasemodel = load( fullfile('./exp', saveFolder, modelName) );
netbasemodel = netbasemodel.net;

netbasemodel.layers(136).block = rmfield(netbasemodel.layers(136).block, 'ignoreAverage');
netbasemodel.layers(135).block = rmfield(netbasemodel.layers(135).block, 'ignoreAverage');

netbasemodel = dagnn.DagNN.loadobj(netbasemodel);
%% modify layers, e.g. removing
netbasemodel.meta.normalization.averageImage = imdb.meta.meanvalue;
netbasemodel.meta.normalization.imageSize = [imdb.meta.height, imdb.meta.width, 1, 1];
%% 1st mean-shift grouping loop
loopNum = 20;

sName_l2norm = 'res7_l2norm';
netbasemodel.vars(netbasemodel.layers(netbasemodel.getLayerIndex(sName_l2norm)).outputIndexes).precious = 1;

% keepLayerName = sprintf('obj_instSeg_reg');
% netbasemodel.layers(netbasemodel.getLayerIndex(keepLayerName)).block.lastLayerName = 'res7_cosSim';
% keepLayerName = sprintf('obj_instSeg_MM');
% netbasemodel.layers(netbasemodel.getLayerIndex(keepLayerName)).block.lastLayerName = 'res7_cosSim';

rmLayerName = 'obj_instSeg_MM';
if ~isnan(netbasemodel.getLayerIndex(rmLayerName))
    netbasemodel.removeLayer(rmLayerName); % remove layer
end
rmLayerName = 'obj_instSeg_reg';
if ~isnan(netbasemodel.getLayerIndex(rmLayerName))
    netbasemodel.removeLayer(rmLayerName); % remove layer
end
rmLayerName = 'res7_cosSim';
if ~isnan(netbasemodel.getLayerIndex(rmLayerName))
    netbasemodel.removeLayer(rmLayerName); % remove layer
end

gt_name =  sprintf('gt_ins');
weight_for_losses = {'obj_instSeg_reg', 1, 'obj_instSeg_MM', 1};
for loopIdx = 1:loopNum
    [netbasemodel, sName, sName_l2norm] = addOneLoop_forMeanShiftGrouping(netbasemodel, sName_l2norm, loopIdx);

    keepLayerName = sprintf('loop%d_meanshift_S_is_XX', loopIdx);
    netbasemodel.layers(netbasemodel.getLayerIndex(keepLayerName)).block.analyzeGradient = false;

    keepLayerName = sprintf('loop%d_meanshift_Y_l2norm', loopIdx);
    netbasemodel.vars(netbasemodel.layers(netbasemodel.getLayerIndex(keepLayerName)).outputIndexes).precious = 1;

    rmLayerName = sprintf('loop%d_meanshift_cosSim', loopIdx);
    netbasemodel.removeLayer(rmLayerName); % remove layer
end
%%

netbasemodel.move('gpu');
% netMat.mode = 'test';
netbasemodel.mode = 'normal';
testList = find(imdb.set==2);
rawCountsAll = {};

semanticMaskMat = [];
instanceMaskMat = [];
imgMat = [];
predInstanceMaskMat_loop0 = [];
predInstanceMaskMat = {};
for loopIdx = 1:loopNum
    predInstanceMaskMat{loopIdx} = [];
end
st = 0;
for i = 1:5%length(testList)
    curBatch = fullfile(imdb.path_to_dataset, imdb.imgList{testList(i)});
    datamat = load(curBatch);

    datamat.semanticMaskMat = reshape(datamat.semanticMaskMat, [size(datamat.semanticMaskMat,1), size(datamat.semanticMaskMat,2), 1, size(datamat.semanticMaskMat,3)]);
    datamat.instanceMaskMat = reshape(datamat.instanceMaskMat, [size(datamat.instanceMaskMat,1), size(datamat.instanceMaskMat,2), 1, size(datamat.instanceMaskMat,3)]);

    imFeed = bsxfun(@minus, datamat.imgMat, imdb.meta.meanvalue);
    inputs = {'input', gpuArray(imFeed)};
    netbasemodel.eval(inputs) ;
    res7_l2norm = gather(netbasemodel.vars(netbasemodel.layers(netbasemodel.getLayerIndex('res7_l2norm')).outputIndexes).value);
    predInstanceMaskMat_loop0(:,:,:,st+1:st+size(datamat.imgMat,4)) = res7_l2norm;

    for loopIdx = 1:loopNum
        curLayerName = sprintf('loop%d_meanshift_Y_l2norm', loopIdx);
        predInstanceMaskMat{loopIdx}(:,:,:,st+1:st+size(datamat.imgMat,4)) = ...
            gather(netbasemodel.vars(netbasemodel.layers(netbasemodel.getLayerIndex(curLayerName)).outputIndexes).value);
    end

    imgMat(:,:,:,st+1:st+size(datamat.imgMat,4)) = datamat.imgMat;
    semanticMaskMat(:,:,:,st+1:st+size(datamat.imgMat,4)) = datamat.semanticMaskMat;
    instanceMaskMat(:,:,:,st+1:st+size(datamat.imgMat,4)) = datamat.instanceMaskMat;
    fprintf('\t%d/%d\n', i, length(testList));

    st = st + size(datamat.imgMat,4);
end
%% grouping results to visualize
num = 25; % the total number of example to show
M = 5; % the number of rows in the panel, each row shows N images as below
N = 5; % the number of columns in the panel, each column shows M images as above
% h = 28;
% w = 28;
panelSZ = round(64/2)*2; % the size (height/width) of the small square image
% stepSize = (panelSZ-h)/2;

% image
showcaseSamples = stitchImage2panel(imgMat, panelSZ, M, N, 1, num);

% ground-truth semantic mask
showcaseSemanticMask= semanticMaskMat(:,:,:,1:num);
showcaseSemanticMask = reshape(showcaseSemanticMask, [numel(showcaseSemanticMask(:,:,1,1)), num]);
showcaseSemanticMask = showStitch(showcaseSemanticMask, [panelSZ panelSZ], M, N, 'whitelines', 'linewidth', 5);

% ground-truth instance mask
showcaseInstanceMask= instanceMaskMat(:,:,:,1:num);
showcaseInstanceMask = reshape(showcaseInstanceMask, [numel(showcaseInstanceMask(:,:,1,1)), num]);
showcaseInstanceMask = showStitch(showcaseInstanceMask, [panelSZ panelSZ], M, N, 'whitelines', 'linewidth', 5);

predInstanceMaskMat_loop0 = stitchImage2panel(predInstanceMaskMat_loop0, panelSZ, M, N, 1, num);
predInstanceMaskMat_loop0 = predInstanceMaskMat_loop0 - min(predInstanceMaskMat_loop0(:));
predInstanceMaskMat_loop0 = predInstanceMaskMat_loop0 ./ max(predInstanceMaskMat_loop0(:));

predInstanceMask = {};
for loopIdx = 1:loopNum
    predInstanceMask{loopIdx} = stitchImage2panel(predInstanceMaskMat{loopIdx}, panelSZ, M, N, 1, num);
    predInstanceMask{loopIdx} = predInstanceMask{loopIdx} - min(predInstanceMask{loopIdx}(:));
    predInstanceMask{loopIdx} = predInstanceMask{loopIdx} ./ max(predInstanceMask{loopIdx}(:));
end

%% save the results
figFolder = './figFolder';
if ~isdir(figFolder)
    mkdir(figFolder);
end
a = strfind(modelName,'epoch');
path_to_save = fullfile(figFolder, [strrep(saveFolder,'/','') '_' modelName(a(1):end-4)]);
if ~isdir(path_to_save)
    mkdir(path_to_save);
end

% summary vidualization
imgFig2 = figure(1);
subWindowH = 2;
subWindowW = 2;
set(imgFig2, 'Position', [100 100 1400 900]) % [1 1 width height]
windowID = 1;
subplot(subWindowH, subWindowW, windowID); imagesc(showcaseSamples); axis off image; title('random samples'); colorbar; windowID = windowID + 1;
subplot(subWindowH, subWindowW, windowID); imagesc(showcaseSemanticMask); axis off image; title('showcase SemanticMask'); caxis([0, 11]); colorbar; windowID = windowID + 1;
subplot(subWindowH, subWindowW, windowID); imagesc(showcaseInstanceMask); axis off image; title('showcase instanceMask'); colorbar; windowID = windowID + 1;
subplot(subWindowH, subWindowW, windowID); imagesc(predInstanceMaskMat_loop0); axis off image; title('predInstanceMask'); colorbar; windowID = windowID + 1;
if flagSaveFig
    export_fig( sprintf('%s/fig00_visualization.jpg', path_to_save) );
end

imgFig2 = figure(2);
subWindowH = 2;
subWindowW = 3;
windowID = 1;
set(imgFig2, 'Position', [100 100 1400 900]) % [1 1 width height]
for loopIdx = 1:6
    subplot(subWindowH, subWindowW, windowID);
    imagesc(predInstanceMask{loopIdx}); axis off image; title(sprintf('loop%d',loopIdx)); caxis([0, 11]); colorbar;
    windowID = windowID + 1;
end
if flagSaveFig
    export_fig( sprintf('%s/fig01_visualization_looping1.jpg', path_to_save) );
end

imgFig3 = figure(3);
subWindowH = 2;
subWindowW = 3;
windowID = 1;
set(imgFig3, 'Position', [100 100 1400 900]) % [1 1 width height]
for loopIdx = 7:12
    subplot(subWindowH, subWindowW, windowID);
    imagesc(predInstanceMask{loopIdx}); axis off image; title(sprintf('loop%d',loopIdx)); caxis([0, 11]); colorbar;
    windowID = windowID + 1;
end
if flagSaveFig
    export_fig( sprintf('%s/fig01_visualization_looping2.jpg', path_to_save) );
end

imgFig4 = figure(4);
subWindowH = 2;
subWindowW = 3;
windowID = 1;
set(imgFig4, 'Position', [100 100 1400 900]) % [1 1 width height]
for loopIdx = 13:18
    subplot(subWindowH, subWindowW, windowID);
    imagesc(predInstanceMask{loopIdx}); axis off image; title(sprintf('loop%d',loopIdx)); caxis([0, 11]); colorbar;
    windowID = windowID + 1;
end
if flagSaveFig
    export_fig( sprintf('%s/fig01_visualization_looping3.jpg', path_to_save) );
end

if flagSaveFig
    showcaseSemanticMask = uint8(showcaseSemanticMask * 255/11);
    showcaseInstanceMask = showcaseInstanceMask*255/5;

    imwrite(showcaseSamples, fullfile(path_to_save, sprintf('fig02_showcaseSamples.bmp')));
    imwrite(showcaseSemanticMask, fullfile(path_to_save, sprintf('fig03_showcaseSemanticMask.bmp')));
    imwrite(showcaseInstanceMask, fullfile(path_to_save, sprintf('fig04_showcaseInstanceMask.bmp')));
    imwrite(predInstanceMaskMat_loop0, fullfile(path_to_save, sprintf('fig05_predInstanceMask-loop0.bmp')));
    imwrite(predInstanceMask{1}, fullfile(path_to_save, sprintf('fig06_predInstanceMask-loop1.bmp')));
    imwrite(predInstanceMask{2}, fullfile(path_to_save, sprintf('fig07_predInstanceMask-loop2.bmp')));
    imwrite(predInstanceMask{3}, fullfile(path_to_save, sprintf('fig08_predInstanceMask-loop3.bmp')));
    imwrite(predInstanceMask{4}, fullfile(path_to_save, sprintf('fig09_predInstanceMask-loop4.bmp')));
    imwrite(predInstanceMask{5}, fullfile(path_to_save, sprintf('fig10_predInstanceMask-loop5.bmp')));

    save(sprintf('%s/results.mat', path_to_save), 'imgMat', 'instanceMaskMat', 'semanticMaskMat', 'predInstanceMaskMat');
end

%% leaving blank
%{
%}

author: aimerykong


text:
function p00_plot ( prob )

%*****************************************************************************80
%
%% P00_PLOT plots the data for any of the tests.
%
%  Discussion:
%
%    For now we assume that the data dimension is 2, so that we are simply
%    creating a single X-Y plot.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the problem index.
%
  prob_num = p00_prob_num;

  if ( prob < 1 || prob_num < prob )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_PLOT - Fatal error!\n' );
    fprintf ( 1, '  Values of PROB must be between 1 and %d.\n', prob_num );
    error ( 'P00_PLOT - Fatal error!' );
  end

  data_num = p00_data_num ( prob );

  dim_num = p00_dim_num ( prob );

  p = p00_data ( prob, dim_num, data_num );

  x = p(1,:);
  y = p(2,:);

  plot ( x, y, 'Linewidth', 2 );
  hold on
  plot ( x, y, 'r.', 'Markersize', 25 )
  grid
  xlabel ( '<---X--->' )
  ylabel ( '<---Y--->' )
  title ( sprintf ( 'TEST\\_INTERP Data Set %d', prob ) );
  hold off

  filename = sprintf ( 'p%02d_plot.png', prob );
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Created graphics file "%s".\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function dfadjustlayout(dffig,showctrl)
%ADJUSTLAYOUT Adjust layout of buttons and graph in figure window

%   $Revision: 1.1.6.5 $  $Date: 2004/01/24 09:35:12 $
%   Copyright 2003-2004 The MathWorks, Inc.

% Get some measurements
fpos = get(dffig,'Position');
fwidth = fpos(3);
fheight = max(1,fpos(4));

% Adjust selection controls at top
hsel = getappdata(dffig,'selectioncontrols');
lmargin = 0;
if ~isempty(hsel)
   % If there are no selection controls, don't try to compute their height
   emax = zeros(1,4);
   allpos = get(hsel, 'Position');
   allpos = vertcat(allpos{:});
   allextent = get(hsel, 'Extent');
   allextent = vertcat(allextent{:});
   maxheight = max(allextent(2:end,4));   % max height of non-frame controls
end
for j=2:2:length(hsel)
   % Adjust label
   p = allpos(j,:);
   e = allextent(j,:);
   p(1) = lmargin + 5;
   p(2) = fheight-1.45*maxheight;
   p(3) = e(3);
   p(4) = maxheight;
   lmargin = p(1) + p(3);
   set(hsel(j),'Position',p);
   emax = max(emax,e);
   labelwidth = e(3);

   % Adjust drop-down
   p = get(hsel(j+1),'Position');
   e = get(hsel(j+1),'Extent');
   p(1) = lmargin + 5;
   p(2) = fheight-1.25*maxheight;
   p(3) = max(e(3),2.25*labelwidth);
   p(4) = maxheight;
   lmargin = p(1) + p(3) + 20;
   set(hsel(j+1),'Position',p);
   emax = max(emax,e);
end
if ~isempty(hsel)
   % Get position of the base of the frame containing these controls
   framebase = max(1, fheight-1.6*maxheight);
   p = [1, framebase, fwidth, 1.6*maxheight];
   set(hsel(1),'Position',p);
else
   % Say the base of this frame is right at the top of the figure
   framebase = fheight;
end

% Adjust all button positions below selection controls
hbuttons = getappdata(dffig,'buttoncontrols');
if isempty(hbuttons)
   % If there are no buttons, say their base is right at the frame
   buttonbase = framebase;
else
   % If there are buttons, compute the height of their base
   nbuttons = length(hbuttons);
   extents = get(hbuttons,'Extent');
   extents = vertcat(extents{:});
   bheight = 1.5 * extents(1,4);  % 1.5 * text height
   gutter = bheight/4;            % between buttons
   margin = bheight/2;            % around text within button
   bwidth = extents(:,3)' + 2*margin;
   totalwidth = sum(bwidth) + gutter*(nbuttons-1);
   startpos = max(0, (fwidth/2) - (totalwidth/2));
   bleft = startpos + [0, cumsum(bwidth + gutter)];
   buttonbase = max(1, framebase-bheight-gutter);
   for j=1:nbuttons
      pos = [bleft(j), buttonbase, bwidth(j), bheight];
      set(hbuttons(j), 'Position',pos);
   end
end

% Position the axes in the remaining area
ax = get(dffig,'CurrentAxes');
axbase = 0.11*fheight;
p1 = max(1, [.13*fwidth,  axbase, .775*fwidth, .9*buttonbase-axbase]);
set(ax,'Units','pixels','Position',p1);
set(ax,'Units','normalized');

if nargin<2
   showctrl = dfgetset('showaxlimctrl');
end
if isequal(showctrl,'on')
   dfaxlimctrl(dffig,'off');
   dfaxlimctrl(dffig,'on');
end

author: zouchuhang

repo: LayoutNet


text:
function colorM = visualRefColormap(maxDose,refDose,numRows,baseStr,topStr)
%Returns the 'visualReference' colormap based on a 'base' colormap baseStr and
%a 'top' colormap topStr.  The base colormap is stretched from 0 to refDose.  Between
%refDose and maxDose (if there are such doses), the colormap follows the colormap
%defined by topStr.  See CERRColormap to see valid colormap strings.
%
%JOD, 16 Jan 03.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%We want to rescale the colormap to an arbitrary number of bins.
%Given the 'visual reference dose', and the fact that we always want a colormap
%with a length of n, then the number of bins up to the visual reference dose are:
%zero dose should map to an index of one.
%D_vr should map to the highest index in the colormap, cM.
%Let m be the number of rows in cM.

%Get base and top unrescaled colormaps
baseM = CERRColorMap(baseStr);
topM  = CERRColorMap(topStr);

doseV = linspace(0,maxDose,numRows);

s = size(baseM,1);
t = size(topM,1);

indexBase = (doseV/refDose) * (s - 1);   %Scale dose from 0 to (s-1)

highIndex = [(doseV/refDose) > 1];

%truncate
indexBase(highIndex) = [];

%now interpolate into the colormap:
baseRescaledM = zeros(length(indexBase),3);
baseRescaledM(:,1) = interp1([0:s-1],baseM(:,1),indexBase)';
baseRescaledM(:,2) = interp1([0:s-1],baseM(:,2),indexBase)';
baseRescaledM(:,3) = interp1([0:s-1],baseM(:,3),indexBase)';

colorM = baseRescaledM;

if any(highIndex)

  topRows = numRows - length(indexBase);
  topRescaledM = zeros(topRows,3);
  indexTop = ((doseV - refDose)/(max(doseV(:)) - refDose)) * (t - 1);  %Scale doses to (t-1)
  indexTop(~highIndex) = [];

  %and the corresponding colormap is
  topRescaledM(:,1) = interp1([0:t-1],topM(:,1),indexTop)';
  topRescaledM(:,2) = interp1([0:t-1],topM(:,2),indexTop)';
  topRescaledM(:,3) = interp1([0:t-1],topM(:,3),indexTop)';

  colorM = [baseRescaledM; topRescaledM];

end

%fini

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function jed = epoch_to_jed_macedonian ( )

%*****************************************************************************80
%
%% EPOCH_TO_JED_MACEDONIAN: epoch of the Macedonian calendar as a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real JED, the Julian Ephemeris Date of the epoch.
%
  jed = 1607708.5;

  return
end


text:
function binomial_test ( )

%*****************************************************************************80
%
%% BINOMIAL_TEST tests BINOMIAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BINOMIAL_TEST:\n' );
  fprintf ( 1, '  A demonstration of the binomial method\n' );
  fprintf ( 1, '  for option valuation.\n' );

  s0 = 2.0;
  e = 1.0;
  r = 0.05;
  sigma = 0.25;
  t1 = 3.0;
  m = 256;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The asset price at time 0, S0    = %f\n', s0 );
  fprintf ( 1, '  The exercise price         E     = %f\n', e );
  fprintf ( 1, '  The interest rate          R     = %f\n', r );
  fprintf ( 1, '  The asset volatility       SIGMA = %f\n', sigma );
  fprintf ( 1, '  The expiry date            T1    = %f\n', t1 );
  fprintf ( 1, '  The number of intervals    M     = %d\n', m );

  c = binomial ( s0, e, r, sigma, t1, m );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The option value is %f\n', c );

  return
end


text:
function [dat] = set_x(dat,inp,ind,fInd)

%   data=set_x(data,x) sets the X matrix
%   data=set_x(data,x,ind) sets the X matrix for specific indexes

if nargin==2 | (isempty(ind)&isempty(fInd)),
  dat.myX=inp;
  if ~(size(inp,1)==length(dat.index))
    dat.index  = [1:size(inp,1)];
  end
  if ~(size(inp,2)==length(dat.findex) )
    dat.findex = [1:size(inp,2)];
  end
  return;
end

if nargin<3 | isempty(ind),
    ind=[1:size(inp,1)];
end
if nargin<4|isempty(fInd),
    fInd=[1:size(inp,2)];
end

if isempty(dat.myX)&nargin>2
    global X;
    dat.myX=X;
    disp('warning -- instantiated some values of x when using globals');
end

dat.myX(ind,fInd)=inp;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%%
%   Author: Yash Bansod
%
% GitHub: <https://github.com/YashBansod>
%
% Demonstrates the calculation of Jacobian for a 3DOF 2D Planar Robot

%% Clear the environment and the command line
clear;
clc;
close all;

%% Calculation of Jacobian for 3DOF 2D Planar Robot
syms t1 t2 t3 L1 L2 L3
orgX = 0;
orgY = 0;

X(1)= orgX+ L1*cos(t1);
Y(1)= orgY+ L1*sin(t1);
X(2)= X(1)+ L2*cos(t1 + t2);
Y(2)= Y(1)+ L2*sin(t1 + t2);
X(3)= X(2)+ L3*cos(t1 + t2 + t3);
Y(3)= Y(2)+ L3*sin(t1 + t2 + t3);

Jacob = [diff(X(3),t1) , diff(X(3), t2) , diff(X(3), t3);
            diff(Y(3), t1), diff(Y(3), t2) , diff(Y(3), t3)]

author: YashBansod

repo: Robotics-Planning-Dynamics-and-Control


text:
function hmm = states_supdate(hmm,hmm_noisy,rho,update)
% update==1, W; update==2, Omega; update==3, sigma; update==4, alpha
K = length(hmm_noisy.state);
Sind = hmm.train.Sind==1;
regressed = sum((hmm.train.Sind==1),1)>0;

for k = 1:K
    if any(update==1) && isfield(hmm_noisy.state(1),'W') && ~isempty(hmm_noisy.state(1).W.Mu_W)
        hmm.state(k).W.Mu_W = (1-rho) * hmm.state(k).W.Mu_W + ...
            rho * hmm_noisy.state(k).W.Mu_W;
        hmm.state(k).W.iS_W = (1-rho) * hmm.state(k).W.iS_W + ...
            rho * hmm_noisy.state(k).W.iS_W;
        if length(size(hmm.state(k).W.S_W))==2 && ... % full, uniquefull, uniqueAR or 1 dimension
               size(hmm.state(k).W.S_W,1)==size(hmm.state(k).W.S_W,2)
            hmm.state(k).W.S_W = inv(hmm.state(k).W.iS_W);
        else % diag or shareddiag
            for n = 1:size(hmm.state(k).W.S_W,1)
                hmm.state(k).W.S_W(n,Sind(:,n),Sind(:,n)) = ...
                    inv(permute(hmm.state(k).W.iS_W(n,Sind(:,n),Sind(:,n)),[2 3 1]));
            end
        end
    elseif any(update==2) && isfield(hmm_noisy.state(1),'Omega')
        hmm.state(k).Omega.Gam_rate = (1-rho) * hmm.state(k).Omega.Gam_rate + ...
            rho * hmm_noisy.state(k).Omega.Gam_rate;
        hmm.state(k).Omega.Gam_shape = (1-rho) * hmm.state(k).Omega.Gam_shape + ...
            rho * hmm_noisy.state(k).Omega.Gam_shape;
        if strcmp(hmm_noisy.train.covtype,'full') || ...
                strcmp(hmm_noisy.train.covtype,'uniquefull') || ...
                strcmp(hmm_noisy.train.covtype,'sharedfull')
            hmm.state(k).Omega.Gam_irate(regressed,regressed) = ...
                inv(hmm.state(k).Omega.Gam_rate(regressed,regressed));
        end
    elseif any(update==3) && isfield(hmm_noisy.state(1),'sigma')
        hmm.state(k).sigma.Gam_rate = (1-rho) * hmm.state(k).sigma.Gam_rate + ...
            rho * hmm_noisy.state(k).sigma.Gam_rate;
        hmm.state(k).sigma.Gam_shape = (1-rho) * hmm.state(k).sigma.Gam_shape + ...
            rho * hmm_noisy.state(k).sigma.Gam_shape;
    elseif any(update==4) && isfield(hmm_noisy.state(1),'alpha')
        hmm.state(k).alpha.Gam_rate = (1-rho) * hmm.state(k).alpha.Gam_rate + ...
            rho * hmm_noisy.state(k).alpha.Gam_rate;
        hmm.state(k).alpha.Gam_shape = (1-rho) * hmm.state(k).alpha.Gam_shape + ...
            rho * hmm_noisy.state(k).alpha.Gam_shape;
    elseif any(update==5) && isfield(hmm_noisy.state(1),'beta')
        hmm.state(k).beta.Gam_rate = (1-rho) * hmm.state(k).beta.Gam_rate + ...
            rho * hmm_noisy.state(k).beta.Gam_rate;
        hmm.state(k).beta.Gam_shape = (1-rho) * hmm.state(k).beta.Gam_shape + ...
            rho * hmm_noisy.state(k).beta.Gam_shape;
    end
end

if any(update==2) && isfield(hmm_noisy,'Omega')
    hmm.Omega.Gam_rate = (1-rho) * hmm.Omega.Gam_rate + ...
        rho * hmm_noisy.Omega.Gam_rate;
    hmm.Omega.Gam_shape = (1-rho) * hmm.Omega.Gam_shape + ...
        rho * hmm_noisy.Omega.Gam_shape;
    if strcmp(hmm_noisy.train.covtype,'uniquefull') || ...
            strcmp(hmm_noisy.train.covtype,'sharedfull')
        hmm.Omega.Gam_irate(regressed,regressed) = ...
            inv(hmm.Omega.Gam_rate(regressed,regressed));
    end
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function dat= proc_commonMedianReference(dat, refChans, rerefChans)
%PROC_COMMONMEDIANREFERENCE - rereference signals to common median reference
%
%Synopsis:
% dat= proc_commonMedianReference(dat, <refChans, rerefChans>)
%
%
%Arguments:
%      dat        - data structure of continuous or epoched data
%      refChans   - channels used as average reference, see util_chanind for format,
%                   default util_scalpChannels(dat)
%      rerefChans - those channels are rereferenced, default refChans
%
%Returns:
%      dat        - updated data structure
%
% SEE util_scalpChannels, util_chanind

% Author: Benjamin Blankertz
dat = misc_history(dat);

if ~exist('refChans','var') || isempty(refChans),
  refChans= util_scalpChannels(dat);
end
if ~exist('rerefChans','var') || isempty(rerefChans), rerefChans= refChans; end

misc_checkType(dat, 'STRUCT(x clab)');
misc_checkType(refChans, 'CELL{CHAR}|CHAR');
misc_checkType(rerefChans, 'CELL{CHAR}|CHAR');

rc= util_chanind(dat, refChans);
rrc= util_chanind(dat, rerefChans);
car= median(dat.x(:,rc,:), 2);
%% this might consume too much memory:
%car= repmat(car, [1 length(rrc) 1]);
%dat.x(:,rrc,:)= dat.x(:,rrc,:) - car;

for cc= rrc,
  dat.x(:,cc,:)= dat.x(:,cc,:) - car;
  dat.clab{cc}= [dat.clab{cc} ' car'];
end

author: bbci


text:
clear all;
rand('seed',0);
randn('seed',0);
fprintf('test mexFistaPathCoding\n');
p=100;
n=1000;
% generate a DAG
G=sprand(p,p,0.05);
G=mexRemoveCyclesGraph(G);
fprintf('\n');

% generate a data matrix
X=randn(n,p);
X=X-repmat(mean(X),[size(X,1) 1]);
X=mexNormalize(X);
Y=randn(n,2);
Y=Y-repmat(mean(Y),[size(Y,1) 1]);
Y=mexNormalize(Y);
W0=zeros(size(X,2),size(Y,2));

% input graph
graph.weights=G;
graph.stop_weights=zeros(1,p);
graph.start_weights=10*ones(1,p);

% FISTA parameters
param.num_threads=-1; % all cores (-1 by default)
param.verbose=true;   % verbosity, false by default
param.lambda=0.005; % regularization parameter
param.it0=1;      % frequency for duality gap computations
param.max_it=100; % maximum number of iterations
param.L0=0.01;
param.tol=1e-4;
param.precision=10000000;
param.pos=false;

fprintf('Square Loss + convex path penalty\n');
param.loss='square';
param.regul='graph-path-conv';
tic
[W1 optim_info]=mexFistaPathCoding(Y,X,W0,graph,param);
t=toc;
fprintf('mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n',mean(optim_info(1,:)),mean(optim_info(3,:)),t,mean(optim_info(4,:)));
num=mexCountConnexComponents(graph.weights,W1(:,1));
fprintf('Num of connected components: %d\n',num);

fprintf('\n');
fprintf('Square Loss + non-convex path penalty\n');
param.loss='square';
param.regul='graph-path-l0';
param.lambda=0.0001; % regularization parameter
param.ista=true;
tic
[W2 optim_info]=mexFistaPathCoding(Y,X,W0,graph,param);
t=toc;
num=mexCountConnexComponents(graph.weights,W2(:,1));
fprintf('Num of connected components: %d\n',num);

fprintf('\n');
fprintf('Note that for non-convex penalties, continuation strategies sometimes perform better:\n');
tablambda=param.lambda*sqrt(sqrt(sqrt(2))).^(20:-1:0);
lambda_orig=param.lambda;
tic
W2=W0;
for ii = 1:length(tablambda)
   param.lambda=tablambda(ii);
   param.verbose=false;
   [W2]=mexFistaPathCoding(Y,X,W2,graph,param);
end
param.verbose=true;
param.lambda=lambda_orig;
[W2 optim_info]=mexFistaPathCoding(Y,X,W2,graph,param);
t=toc;
num=mexCountConnexComponents(graph.weights,W2(:,1));
param.ista=false;
fprintf('Num of connected components: %d\n',num);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [M_, S_]=Prior2Posterior(M,Q,M_Q,S,G,S_G)

M_=M+S*Q'*inv(Q*S*Q')*(M_Q-Q*M);
S_=S+(S*G')*( inv(G*S*G')*S_G*inv(G*S*G') -inv(G*S*G') )*(G*S);


text:
function p10_story ( )

%*****************************************************************************80
%
%% P10_STORY prints the "story" for problem 10.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  100 uniformly random X values between -2 and 5 were selected, and\n' );
  fprintf ( 1, '  the formula Y = 2 + 5 * X + 10 * N(0,1) was evaluated, where N(0,1)\n' );
  fprintf ( 1, '  represents random normal values with 0 mean and unit variance.\n' );

  return
end


text:
% OtTableDelete Delete the current record.
%
% Usage:
%    OtTableDelete(htable)

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-04-11 Created (Johan Meijdam, VORtech)

function iret = OtTableDelete(htable)
   iret = Ot('table.delete', htable);
end


text:
%kkgconst 'Generate Object Containing Constant Value Data'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kgconst.pane file
%
% Parameters:
% Integer: wsize 'Width ', default: 1: 'Width size of constant value data'
% Integer: hsize 'Height ', default: 1: 'Height size of constant value data'
% Integer: dsize 'Depth ', default: 1: 'Depth size of constant value data'
% Integer: tsize 'Time ', default: 1: 'Time size of constant value data'
% Integer: esize 'Elements ', default: 1: 'Element dimension size of constant value data'
% OutputFile: o 'Output', required: 'Output file containing constant value data'
%
% Example: o = kkgconst( {'wsize',1;'hsize',1;'dsize',1;'tsize',1;'esize',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kgconst - Generate Object Containing Constant Value Data
%
%  DESCRIPTION
% .I kgconst
% creates a data object with a constant valued \fBvalue segment" of dimension
% Width * Height * Depth * Time * Elements.
% The data type of the value segment can be
% \fBbit, byte, unsigned byte, short, unsigned short, integer, unsigned integer,
% long, unsigned long, float, double, complex" and \fBdouble complex\fP.
% The real, and imaginary constant levels are specified using the level
% parameters (real, imag).
%
%
%
%  EXAMPLES
% kgconst -o const -wsize 256 -hsize 256 -esize 3 -type integer -real 10
%
% Creates a data object called const of size width = 256, height = 256
% where each vector contains 3 elements.  The data storage type is integer,
% and the constant is 10.
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkgconst(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kkgconst(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'wsize', 1;'hsize', 1;'dsize', 1;'tsize', 1;'esize', 1;'o', '__output'};
maxval={2,2,2,2,2,0};
minval={2,2,2,2,2,0};
istoggle=[1,1,1,1,1,0];
was_set=istoggle * 0;
paramtype={'Integer','Integer','Integer','Integer','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kgconst"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [SenErrDef, IniErrDef, ObsErrDef, ini_pva, wander, dt]=sys_path(dir_name,record)
rstat=[134 26 345 476 5456;654 7345 8345 9678 140];       %control how random the errors are generated
%rstat=rand(2,5)*30000;
dt=1/100;
D2R=pi/180;

%%%imu error parameters (discrete time param. defined at 100Hz)
%Accelerometers (defined in m-sec^2)
SenErrDef(1).A=1;
SenErrDef(1).B=3e-006;
SenErrDef(1).C=1;
SenErrDef(1).D=1e-3/sqrt(0.01);
SenErrDef(1).sP=sqrt(5.1e-007);
SenErrDef(1).tparam=[];

SenErrDef(2)=SenErrDef(1);
SenErrDef(3)=SenErrDef(1);

%Gyroscopes (defined in rad/sec)
SenErrDef(4).A=0.999980000199999;
SenErrDef(4).B=1.9999800001658e-006;
SenErrDef(4).C=1;
SenErrDef(4).D=0.0006/sqrt(0.01);
SenErrDef(4).sP=sqrt(1e-007);
SenErrDef(4).tparam=[];  %first param=sP of temperature scale factor (RC), second=B of temperature RW

SenErrDef(5)=SenErrDef(4);
SenErrDef(6)=SenErrDef(4);

% SenErrDef(6).A=1;
% SenErrDef(6).B=0;
% SenErrDef(6).C=1;
% SenErrDef(6).D=0;
% SenErrDef(6).sP=0;
% SenErrDef(1)=SenErrDef(6);
% SenErrDef(2)=SenErrDef(6);
% SenErrDef(3)=SenErrDef(6);
% SenErrDef(4)=SenErrDef(6);
% SenErrDef(5)=SenErrDef(6);

%%initial err defs
IniErrDef.pos_sP=diag([1/6e6,1/6e6,1]);
IniErrDef.vel_sP=eye(3)*0.01;
IniErrDef.att_sP=eye(3)*(10*pi/180);
IniErrDef.wander_sP=eye(2)*2;

%3D simulation in NED
%initial nav values
vel_b=zeros(3,1);
att_n=[-3;5;110]*D2R;
pos_n=[39.88864*D2R;32.78002*D2R;1102]; %51.081755-114.1360749

%%observation err defs
[Rn, Re, g, sL, cL, WIE_E]=geoparam_v000(pos_n);
Sn=Rn+pos_n(3);
Se=(Re+pos_n(3))*cos(pos_n(1));
ObsErrDef.sR=diag([0.5/Sn, 0.5/Se, 1.5]); %position error std (m)

%%%%%%%%%%% Generate path %%%%%%%%%%%%%%%%
if (record)
    %%motion definitions
    mot_def(1,:)=[1 0 0 0 0 5];
    mot_def(2,:)=[5 0 0 0 2 10];
    mot_def(3,:)=[1 0 0 0 0 5];
    mot_def(4,:)=[3 0 0 90*D2R 0 15];
    mot_def(5,:)=[1 0 0 0 0 5];
    mot_def(6,:)=[3 0 0 90*D2R 0 15];
    mot_def(7,:)=[1 0 0 0 0 5];
    mot_def(8,:)=[3 0 0 -90*D2R 0 15];
    mot_def(9,:)=[1 0 0 0 0 10];
    mot_def(10,:)=[3 0 0 -90*D2R 0 15];
    mot_def(11,:)=[1 0 0 0 0 5];
    mot_def(12,:)=[5 0 0 0 0 10];

    PathGen_v002(dir_name, [pos_n, vel_b, att_n], mot_def, [1 1/dt;1 0.1], 0, [], [], [], rstat(:,1));

    % %%%Add error
    % %Imu errors
    AddIMUErr_v000([dir_name 'mimu.bin'], [dir_name 'imu.bin'], [dir_name 'imuerr.bin'], 7, 2:7, SenErrDef, [], rstat(:,2));
    AddObsErr_v000([dir_name 'gps.bin'], 7, [dir_name 'obs_pos.bin'], [5 6 7], ObsErrDef, rstat(:,3));

    %AddIMUErr_v000([dir_name 'mimu.bin'], [dir_name 'imu.bin'], [dir_name 'imuerr.bin'], 7, 2:7, [], [], rstat(:,2));
    %AddObsErr_v000([dir_name 'gps.bin'], 7, [dir_name 'obs_pos.bin'], [5 6 7], [], rstat(:,3));
end

%initialization errors
randn('state',rstat(:,4));
pos_err=IniErrDef.pos_sP*randn(3,1);
vel_err=IniErrDef.vel_sP*randn(3,1);
att_err=IniErrDef.att_sP*randn(3,1);
Cerr=euler2dcm_v000(att_err);
Cbg=euler2dcm_v000(att_n);
Cbg_err=Cerr'*Cbg;
att_ini=dcm2euler_v000(Cbg_err);

ini_pva=[pos_n+pos_err vel_b+vel_err att_ini];  %%Attitude will be initialized using acc data
%wander=[sin(att_ini(3));cos(att_ini(3))];
wander=[sin(105*pi/180);cos(105*pi/180)];

repo: ekfmonoslam


text:
%% FAST Algorithm for Corner Detection
%
% In this demo, we will:
%
% * understand the basics of FAST algorithm
% * find corners using OpenCV functionalities for FAST algorithm.
%
% Sources:
%
% * <https://docs.opencv.org/3.2.0/df/d0c/tutorial_py_fast.html>
%

%% Theory
%
% We saw several feature detectors and many of them are really good. But when
% looking from a real-time application point of view, they are not fast
% enough. One best example would be SLAM (Simultaneous Localization and
% Mapping) mobile robot which have limited computational resources.
%
% As a solution to this, FAST (Features from Accelerated Segment Test)
% algorithm was proposed by Edward Rosten and Tom Drummond in their paper
% "Machine learning for high-speed corner detection" in 2006 (later revised it
% in 2010). A basic summary of the algorithm is presented below. Refer to the
% original paper for more details (All the images are taken from original
% paper).
%
%% Feature Detection using FAST
%
% * Select a pixel $p$ in the image which is to be identified as an interest
%   point or not. Let its intensity be $I_p$.
% * Select appropriate threshold value $t$.
% * Consider a circle of 16 pixels around the pixel under test (See the image
%   below)
%
% <<https://docs.opencv.org/3.2.0/fast_speedtest.jpg>>
%
% * Now the pixel $p$ is a corner if there exists a set of $n$ contiguous
%   pixels in the circle (of 16 pixels) which are all brighter than $I_p + t$,
%   or all darker than $I_p - t$ (Shown as white dash lines in the above
%   image). $n$ was chosen to be 12.
% * A *high-speed test* was proposed to exclude a large number of non-corners.
%   This test examines only the four pixels at 1, 9, 5 and 13 (First 1 and 9
%   are tested if they are too brighter or darker. If so, then checks 5 and
%   13). If $p$ is a corner, then at least three of these must all be brighter
%   than $I_p + t$ or darker than $I_p - t$. If neither of these is the case,
%   then $p$ cannot be a corner. The full segment test criterion can then be
%   applied to the passed candidates by examining all pixels in the circle.
%
% This detector in itself exhibits high performance, but there are several
% weaknesses:
%
% * It does not reject as many candidates for |n < 12|.
% * The choice of pixels is not optimal because its efficiency depends on
%   ordering of the questions and distribution of corner appearances.
% * Results of high-speed tests are thrown away.
% * Multiple features are detected adjacent to one another.
%
% First 3 points are addressed with a machine learning approach. Last one is
% addressed using non-maximal suppression.
%
%% Machine Learning a Corner Detector
%
% * Select a set of images for training (preferably from the target
%   application domain)
% * Run FAST algorithm in every images to find feature points.
% * For every feature point, store the 16 pixels around it as a vector. Do it
%   for all the images to get feature vector $P$.
% * Each pixel (say $x$) in these 16 pixels can have one of the following
%   three states:
%
% <<https://docs.opencv.org/3.2.0/fast_eqns.jpg>>
%
% * Depending on these states, the feature vector $P$ is subdivided into 3
%   subsets, $P_d$, $P_s$, $P_b$.
% * Define a new boolean variable, $K_p$, which is true if $p$ is a corner and
%   false otherwise.
% * Use the ID3 algorithm (decision tree classifier) to query each subset
%   using the variable $K_p$ for the knowledge about the true class. It
%   selects the $x$ which yields the most information about whether the
%   candidate pixel is a corner, measured by the entropy of $K_p$.
% * This is recursively applied to all the subsets until its entropy is zero.
% * The decision tree so created is used for fast detection in other images.
%
%% Non-maximal Suppression
%
% Detecting multiple interest points in adjacent locations is another problem.
% It is solved by using Non-maximum Suppression.
%
% * Compute a score function, $V$ for all the detected feature points. $V$ is
%   the sum of absolute difference between $p$ and 16 surrounding pixels
%   values.
% * Consider two adjacent keypoints and compute their $V$ values.
% * Discard the one with lower $V$ value.
%
%% Summary
%
% It is several times faster than other existing corner detectors.
%
% But it is not robust to high levels of noise. It is dependent on a
% threshold.
%

%% FAST Feature Detector in OpenCV
%
% It is called as any other feature detector in OpenCV. If you want, you can
% specify the threshold, whether non-maximum suppression to be applied or not,
% the neighborhood to be used etc.
%
% For the neighborhood, three flags are defined, |TYPE_5_8|, |TYPE_7_12|, and
% |TYPE_9_16|. Below is a simple code on how to detect and draw the FAST
% feature points.
%

%%
% load source image
img = cv.imread(fullfile(mexopencv.root(),'test','blox.jpg'), 'Grayscale',true);

%%
% create FAST object with default values, and find keypoints
fast = cv.FastFeatureDetector();
kp1 = fast.detect(img);
fprintf('%d keypoints with NonmaxSuppression\n', numel(kp1));

%%
% disable |NonmaxSuppression| and find keypoints
fast.NonmaxSuppression = false;
kp2 = fast.detect(img);
fprintf('%d keypoints without NonmaxSuppression\n', numel(kp2));

%%
% draw keypoints from both cases, and show results
out1 = cv.drawKeypoints(img, kp1, 'Color',[255 0 0]);
out2 = cv.drawKeypoints(img, kp2, 'Color',[255 0 0]);
subplot(121), imshow(out1), title('FAST w/ NonmaxSuppression')
subplot(122), imshow(out2), title('FAST w/o NonmaxSuppression')

%% Additional Resources
%
% * Edward Rosten and Tom Drummond, "Machine learning for high speed corner
%   detection" in 9th European Conference on Computer Vision, vol. 1, 2006,
%   pp. 430-443.
% * Edward Rosten, Reid Porter, and Tom Drummond, "Faster and better: a
%   machine learning approach to corner detection" in IEEE Trans. Pattern
%   Analysis and Machine Intelligence, 2010, vol 32, pp. 105-119.
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [p1,p2]=rubberbandbox(varargin)
% Function to draw a rubberband box and return the start and end points
% Usage: [p1,p2]=rubberbox;     uses current axes
%        [p1,p2]=rubberbox(h);  uses axes refered to by handle, h
% Based on an idea of Sandra Martinka's Rubberline
% Written/Edited by Bob Hamans (B.C.Hamans@student.tue.nl)
% 02-04-2003

%Check for optional argument
switch nargin
case 0
  h=gca;
case 1
  h=varargin{1};
  axes(h);
otherwise
  disp('Too many input arguments.');
end

% Get current user data
cudata=get(gcf,'UserData');
hold on;
% Wait for left mouse button to be pressed
k=waitforbuttonpress;
p1=get(h,'CurrentPoint');       %get starting point
p1=p1(1,1:2);                   %extract x and y
lh=plot(p1(1),p1(2),'-r');      %plot starting point
% Save current point (p1) data in a structure to be passed
udata.p1=p1;
udata.h=h;
udata.lh=lh;
% Set gcf object properties 'UserData' and call function 'wbmf' on mouse motion.
set(gcf,'UserData',udata,'WindowButtonMotionFcn','wbmf','DoubleBuffer','on');
k=waitforbuttonpress;

% Get data for the end point
p2=get(h,'Currentpoint');       %get end point
p2=p2(1,1:2);                   %extract x and y
set(gcf,'UserData',cudata,'WindowButtonMotionFcn','','DoubleBuffer','off'); %reset UserData, etc..
delete(lh);

function wbmf %window motion callback function
utemp=get(gcf,'UserData');
ptemp=get(utemp.h,'CurrentPoint');
ptemp=ptemp(1,1:2);
% Use 5 point to draw a rectangular rubberband box
set(utemp.lh,'XData',[ptemp(1),ptemp(1),utemp.p1(1),utemp.p1(1),ptemp(1)],'YData',[ptemp(2),utemp.p1(2),utemp.p1(2),ptemp(2),ptemp(2)]);


text:
function [nodeRho, elemRho, cur, area, v] = curvaturedensity(node,elem, normal, gamma)

n = normal(elem(:,1),:) + normal(elem(:,2),:) + normal(elem(:,3),:);

v12 = node(elem(:,2),:) - node(elem(:,1),:);
v13 = node(elem(:,3),:) - node(elem(:,1),:);

v = cross(v12,v13,2);

area = 0.5*sqrt(sum(v.^2,2));

v = v./[2*area,2*area,2*area];

elemRho = abs(9 - sum(n.^2,2))./area;
elemRho = elemRho +sqrt(eps);

for i = 1:12
    nodeRho = accumarray(elem(:),repmat(elemRho.*area,3,1),[size(node,1),1]);
%     showsolution(node,elem,nodeRho);
%     axis equal;
%     colorbar;
%     pause(0.1);
    b = accumarray(elem(:),[area;area;area],[size(node,1),1]);
    nodeRho =nodeRho./b;
    elemRho = (nodeRho(elem(:,1)) + nodeRho(elem(:,2)) + nodeRho(elem(:,3)))/3;
end

cur = nodeRho;
elemRho = (elemRho/max(elemRho)).^(gamma);
nodeRho = (nodeRho/max(nodeRho)).^(gamma);

repo: ifem