text:
function classifier = lapsvmp(options,data)
% {lapsvmp} trains a Laplacian SVM classifier in the primal.
%
%      classifier = lapsvmp(options,data)
%
%      options: a structure with the following fields
%               options.gamma_A: regularization parameter (ambient norm)
%               options.gamma_I: regularization parameter (intrinsic norm)
%
%               [optional fields]
%               options.Cg: {0,1} i.e. train with Newton's method or PCG
%                           (default=0)
%               options.MaxIter: maximum number of iterations (default=200)
%               options.Hinge: {0,1} i.e. train a LapSVM (1) or LapRLSC (0)
%                              (default=1)
%               options.UseBias: {0,1} i.e. use or not a bias (default=0)
%               options.InitAlpha: if it's 0, the initial weights are null;
%                                  if it's <0, they are randomly taken;
%                                  otherwise it is the initial vector
%                                  (default=0)
%               options.InitBias: the initial bias (default=0)
%               options.NewtonLineSearch: {0,1} i.e. use or not exact line
%                                         search with Newton's method
%                                         (default=0).
%               options.NewtonCholesky: {0,1} i.e. use or not Cholesky
%                                       factorization and rank 1 updates of
%                                       the Heassian (default=0 for LapSVM.
%                                       If set to 1, the  Hessian must be
%                                       positive definite).
%               options.CgStopType: {0,1,2,3,4,5,6,7,-1}
%                                   the data-based stopping criterion of cg
%                                   only (default=0), where:
%                                   0: do not stop
%                                   1: stability stop
%                                   2: validation stop
%                                   3: stability & validation stop
%                                   4: gradient norm (normalized)
%                                   5: prec. gradient norm (normalized)
%                                   6: mixed gradient norm (normalized)
%                                   7: relative objective function change
%                                  -1: debug (it saves many stats)
%               options.CgStopIter: number of cg iters after which check
%                                   the stopping condition.
%               options.CgStopParam: the parameter for the selected
%                                    CgStopType. If CgStopType is:
%                                    0: ignored
%                                    1: percentage [0,1] of tolerated
%                                       different decisions between two
%                                       consecutive checks
%                                    2: percentage [0,1] of error rate
%                                       decrement required  between two
%                                       consecutive checks
%                                    3: the two params above (i.e. it is an
%                                       array of two elements)
%                                    4: the minimum gradient norm
%                                    5: the minimum prec. gradient norm
%                                    6: the minimum mixed gradient norm
%                                    7: relative objective function change
%                                       between two consecutive checks
%                                    -1: ignored
%                                    see the code for default values.
%               options.Verbose: {0,1} (default=1)
%
%      data: a structure with the following fields
%            data.X: a N-by-D matrix of N D-dimensional training examples
%            data.K: a N-by-N kernel Gram matrix of N training examples
%            data.Y: a N-by-1 label vector in {-1,0,+1}, where 0=unlabeled
%                    (it is in {0,+1} in the case of One-Class SVM/LapSVM.
%            data.L: a N-by-N matrix of the Laplacian
%
%            [other fields]
%            data.Kv: a V-by-N kernel Gram matrix of V validation examples,
%                     required if CgStopType is 2 or 3 (validation check)
%            data.Yv: a V-by-1 vector of {-1,+1} labels for the validation
%                     examples, required if CgStopType is 2 or 3
%                     (validation check)
%
%      classifier: structure of the trained classifier (see the
%                  'saveclassfier' function).
%
% Author: Stefano Melacci (2012)
%         mela@dii.unisi.it
%         * the One-Class extension is joint work with Salvatore Frandina,
%           salvatore.frandina@gmail.com
%         * the original code structure was based on the primal SVM code of
%           Olivier Chapelle, olivier.chapelle@tuebingen.mpg.de

n=length(data.Y);
nn=nnz(data.Y==-1);

% initializing option structure with default values
if ~isfield(options,'Verbose'),           options.Verbose=1; end
if ~isfield(options,'Hinge'),             options.Hinge=1; end
if ~isfield(options,'Cg'),                options.Cg=0; end
if ~isfield(options,'MaxIter'),           options.MaxIter=200; end
if ~isfield(options,'UseBias'),           options.UseBias=0; end
if ~isfield(options,'InitAlpha'),         options.InitAlpha=false; end
if ~isfield(options,'InitBias'),          options.InitBias=0; end
if ~isfield(options,'NewtonLineSearch'),  options.NewtonLineSearch=0; end
if ~isfield(options,'NewtonCholesky'),    options.NewtonCholesky=0; end
if ~isfield(options,'CgStopType'),        options.CgStopType=0; end
switch options.CgStopType
    case 0 % none
        if ~isfield(options,'CgStopParam') || ...
            isempty(options.CgStopParam), options.CgStopParam=0; end
        if ~isfield(options,'CgStopIter'),
            options.CgStopIter=options.MaxIter+1; end
    case 1 % stability
        if ~isfield(options,'CgStopParam') || ...
            isempty(options.CgStopParam), options.CgStopParam=0.015; end
        if ~isfield(options,'CgStopIter') || isempty(options.CgStopIter),
            options.CgStopIter=round(sqrt(n)/2); end
    case 2 % validation
        v=length(data.Yv);
        if ~isfield(options,'CgStopParam') || ...
            isempty(options.CgStopParam), options.CgStopParam=1/v; end
        if ~isfield(options,'CgStopIter') || isempty(options.CgStopIter),
            options.CgStopIter=round(sqrt(n)/2); end
    case 3 % stability & validation
        v=length(data.Yv);
        if ~isfield(options,'CgStopParam') || ...
            isempty(options.CgStopParam), options.CgStopParam=...
                                                           [1/v,0.015];
        end
        if ~isfield(options,'CgStopIter'),
            options.CgStopIter=round(sqrt(n)/2); end
    case 4 % gradient norm
        if ~isfield(options,'CgStopParam') || ...
            isempty(options.CgStopParam), options.CgStopParam=1e-8; end
        if ~isfield(options,'CgStopIter') || isempty(options.CgStopIter),
            options.CgStopIter=1; end
    case 5 % preconditioned gradient norm
        if ~isfield(options,'CgStopParam') || ...
            isempty(options.CgStopParam), options.CgStopParam=1e-8; end
        if ~isfield(options,'CgStopIter') || isempty(options.CgStopIter),
            options.CgStopIter=1; end
    case 6 % mixed gradient norm
        if ~isfield(options,'CgStopParam') || ...
            isempty(options.CgStopParam), options.CgStopParam=1e-8; end
        if ~isfield(options,'CgStopIter') || isempty(options.CgStopIter),
            options.CgStopIter=1; end
    case 7 % relative objective function decrease
        if ~isfield(options,'CgStopParam') || ...
            isempty(options.CgStopParam), options.CgStopParam=1e-6; end
        if ~isfield(options,'CgStopIter') || isempty(options.CgStopIter),
            options.CgStopIter=1; end
    case -1 % debug
        if ~isfield(options,'CgStopParam') || ...
            isempty(options.CgStopParam), options.CgStopParam=0; end
        if ~isfield(options,'CgStopIter') || isempty(options.CgStopIter),
            options.CgStopIter=options.MaxIter+1; end
    otherwise
        error('Invalid CgStopType.');
end
if nn==0, oc=1; else oc=0; end

% initial alpha vector
if length(options.InitAlpha)>1
    alpha=options.InitAlpha;
else
    if options.InitAlpha==0
        alpha=[];
    else
        alpha=randn(n,1);
        alpha=alpha./norm(alpha);
    end
end

% checking common error conditions
if oc==1 && ~options.UseBias
    error('One-Class SVM requires the UseBias option to be turned on.');
end
if options.UseBias && options.LaplacianNormalize
    error(['The current implementation does not support a normalized ' ...
           'Laplacian when the UseBias options is turned on.']);
end

% initial bias
b=options.InitBias;

switch options.Cg
    case 0
        % newton
        [alpha,b,t,sec,lsiters]=newton(options,data,alpha,b,oc);
        stats=[];
    case 1
        % pcg
        [alpha,b,t,sec,lsiters,stats]=pcg(options,data,alpha,b,oc);
    otherwise
        error('Invalid solver specified in the field .Cg');
end

svs=find(alpha~=0);
classifier=saveclassifier('lapsvmp',svs,alpha(svs), ...
                          data.X(svs,:),b,options,sec,t,lsiters,stats);

function [alpha,b,t,sec,lsiters] = newton(options,data,alpha,b,oc)
% {newton} trains the classifier using the Newton's method.

tic
n=length(data.Y);
labeled=data.Y~=0;
l=nnz(labeled);
gamma_A=options.gamma_A;
gamma_I=options.gamma_I;

% initial seeding
if isempty(alpha)
    alpha=zeros(n,1);
    Kalpha=zeros(n,1);
else
    Kalpha=data.K*alpha;
end

t=0;
lr=0;
sv=false(n,1);
if nargout>4, lsiters=zeros(options.MaxIter,1); end

if gamma_I~=0, LK=data.L*data.K; end

while 1

    if options.Hinge
        sv_prev=sv;
        hloss=sparse([],[],[],n,1,l);
        hloss(labeled)=1-data.Y(labeled).*(Kalpha(labeled,:)+b);
        sv=hloss>0;
        nsv=nnz(sv);
    else
        sv_prev=sv;
        sv=labeled;
        nsv=l;
    end

    if options.Verbose
        if ~options.Hinge
            hloss=sparse([],[],[],n,1,l);
            hloss(labeled)=1-data.Y(labeled).*(Kalpha(labeled,:)+b);
        end
        if gamma_I~=0
            obj=(gamma_A*alpha'*Kalpha+sum(hloss(sv).^2)+...
                gamma_I*Kalpha'*data.L*Kalpha+oc*b)/2;
        else
            obj=(gamma_A*alpha'*Kalpha+sum(hloss(sv).^2)+oc*b)/2;
        end

        fprintf('[t=%d] obj=%f nev=%d lr=%.4f\n', [t full(obj) nsv lr]);
    end

    % goal conditions
    if t>=options.MaxIter, break, end
    if isequal(sv_prev,sv), break, end

    t=t+1;

    IsvK=sparse([],[],[],n,n,nsv*n);
    IsvK(sv,:)=data.K(sv,:);

    % computing new alphas
    onev=ones(1,n);

    if gamma_I==0 % SVM (sparse solution)
        if options.UseBias
            alpha_new=zeros(n,1);
            alpha_b_new=[0,onev(1:nsv);onev(1:nsv)',...
                         gamma_A*speye(nsv)+IsvK(sv,sv)]\ ...
                         [oc/(2*gamma_A);data.Y(sv)];
            alpha_new(sv)=alpha_b_new(2:end);
            b_new=alpha_b_new(1);
        else
            alpha_new=zeros(n,1);
            alpha_new(sv)=(gamma_A*speye(nsv)+IsvK(sv,sv))\data.Y(sv);
            b_new=0;
        end

    else % LapSVM

        % inversion by factorization
        if options.NewtonCholesky
            if t==1
                % compute the Cholesky factorization of the Hessian
                if options.UseBias
                    sumKsv=sum(data.K(sv,:));
                    hess=chol([nsv,sumKsv;sumKsv',...
                               data.K*(gamma_A*speye(n)+IsvK+gamma_I*LK)]);
                    alpha_b_new=hess\(hess'\([sum(data.Y(sv))-oc/2; ...
                                              data.K(:,sv)*data.Y(sv)]));
                    alpha_new=alpha_b_new(2:end);
                    b_new=alpha_b_new(1);
                else
                    hess=chol(data.K*(gamma_A*speye(n)+IsvK+gamma_I*LK));
                    alpha_new=hess\(hess'\(data.K(:,sv)*data.Y(sv)));
                    b_new=0;
                end
                LK=[];
            else
                % update the Cholesky factorization of the Hessian
                sv_diff=~(sv&sv_prev);
                sv_add=find(sv&sv_diff)';
                sv_rem=find(sv_prev&sv_diff)';
                if options.UseBias
                    if ~isempty(sv_add)
                        for i=sv_add
                            hess=cholupdate(hess,[1;data.K(:,i)],'+');
                        end
                    end
                    if ~isempty(sv_rem)
                        for i=sv_rem
                            hess=cholupdate(hess,[1;data.K(:,i)],'-');
                        end
                    end
                    alpha_b_new=hess\(hess'\([sum(data.Y(sv))-oc/2; ...
                                              data.K(:,sv)*data.Y(sv)]));
                    alpha_new=alpha_b_new(2:end);
                    b_new=alpha_b_new(1);
                else
                    if ~isempty(sv_add)
                        for i=sv_add
                            hess=cholupdate(hess,data.K(:,i),'+');
                        end
                    end
                    if ~isempty(sv_rem)
                        for i=sv_rem
                            hess=cholupdate(hess,data.K(:,i),'-');
                        end
                    end
                    alpha_new=hess\(hess'\(data.K(:,sv)*data.Y(sv)));
                    b_new=0;
                end
            end
        else
            % inversion without factorization
            IsvY=sparse([],[],[],n,1,nsv);
            IsvY(sv)=data.Y(sv);
            if options.UseBias
                alpha_b_new=([0,onev;sv,...
                             gamma_A*speye(n)+IsvK+gamma_I*LK])\ ...
                             [oc/(2*gamma_A);IsvY];
                alpha_new=alpha_b_new(2:end);
                b_new=alpha_b_new(1);
            else
                alpha_new=(gamma_A*speye(n)+IsvK+gamma_I*LK)\IsvY;
                b_new=0;
            end

        end
    end

    % step
    if options.NewtonLineSearch && (options.Hinge || nnz(alpha)>0)
        step=alpha_new-alpha;
        step_b=b_new-b;

        [lr,Kalpha,lsi]=linesearch(data,labeled,step,step_b,Kalpha,b,...
                                   gamma_A,gamma_I,[],[],options.Hinge,...
                                   oc);
        alpha=alpha+lr*step;
        b=b+lr*step_b;
        lsiters(t)=lsi;
    else
        alpha=alpha_new;
        b=b_new;
        lr=1;
        lsiters(t)=0;
        Kalpha=data.K*alpha;
    end
end

if nargout>4, lsiters=lsiters(1:t); end
sec=toc;

function [alpha,b,t,sec,lsiters,stats] = pcg(options,data,alpha,b,oc)
% {pcg} trains the classifier using preconditioned conjugate gradient.

tic
n=length(data.Y);
labeled=data.Y~=0;
unlabeled=data.Y==0;
l=nnz(labeled);
u=n-l;

if isfield(data,'Yv')
    v=length(data.Yv);
end
gamma_A=options.gamma_A;
gamma_I=options.gamma_I;

if isempty(alpha)
    alpha=zeros(n,1);
    Kalpha=zeros(n,1);
    if gamma_I~=0, LKalpha=zeros(n,1); else LKalpha=[]; end
    go=data.Y-b*labeled;
    obj0=(sum((1-data.Y(labeled)*b*options.UseBias).^2)+oc*b)/2;
    if options.UseBias, go_b=sum(data.Y(labeled)-b)-oc/2;
    else go_b=0; b=0; end
else
    Kalpha=data.K*alpha;
    if gamma_I~=0, LKalpha=data.L*Kalpha; else LKalpha=[]; end
    out=sparse([],[],[],n,1,l);
    out(labeled)=Kalpha(labeled)+b;
    if options.Hinge
        sv=false(n,1);
        sv(labeled)=(data.Y(labeled).*out(labeled)<1);
    else
        sv=labeled;
    end
    go=-gamma_A*alpha;
    if gamma_I~=0, go=go-gamma_I*LKalpha; end
    obj0=(Kalpha'*(-g0)+sum((out(sv)-data.Y(sv)).^2)+oc*b)/2;
    go(sv)=go(sv)-(out(sv)-data.Y(sv));

    if options.UseBias, go_b=sum(data.Y(sv)-Kalpha(sv)-b)-oc/2;
    else go_b=0; b=0; end

end
d=go; % initial search direction
d_b=go_b;
Kgo=data.K*go;
Kstep=Kgo;

t=0;
stats=[];

switch options.CgStopType
    case 4, ng0=sqrt(sum(Kgo.^2)+go_b^2);
    case 5, ngp0=sqrt(sum(go.^2)+go_b^2);
    case 6, ngm0=sqrt(sum(Kgo'*go)+go_b^2);
    case -1
        stats=zeros(options.MaxIter+1,5+n+1);
        ng0=sqrt(sum(Kgo.^2)+go_b^2);
        ngp0=sqrt(sum(go.^2)+go_b^2);
        ngm0=sqrt(sum(Kgo'*go)+go_b^2);
        stats(t+1,:)=[t,obj0/obj0,ng0/ng0,ngp0/ngp0,ngm0/ngm0,alpha',b];
end

if nargout>3, lsiters=zeros(options.MaxIter,1); end
valerr_prev=1;
obj_prev=obj0;
yfx_unlabeled_prev=false(u,1);

while 1
    t=t+1;

    % goal condition: maximum number of iterations
    if t>options.MaxIter, t=t-1; break, end

    % do an exact line search
    [lr,Kalpha,lsi,LKalpha]=linesearch(data,labeled,d,d_b,Kalpha,b,...
                                       gamma_A,gamma_I,Kstep,LKalpha,...
                                       options.Hinge,oc);

    % goal condition: converged to optimal solution
    if lr==0, t=t-1; break, end

    alpha=alpha+lr*d;
    b=b+lr*d_b;
    lsiters(t)=lsi;

    % compute new precgradient and objective
    out=sparse([],[],[],n,1,l);
    out(labeled)=Kalpha(labeled)+b;

    if options.Hinge
        sv=false(n,1);
        sv(labeled)=(data.Y(labeled).*out(labeled)<1);
    else
        sv=labeled;
    end

    g=gamma_A*alpha;
    if gamma_I~=0, g=g+gamma_I*LKalpha; end

    if options.Verbose || options.CgStopType==7 || options.CgStopType==-1
        obj=(Kalpha'*(g)+sum((out(sv)-data.Y(sv)).^2)+oc*b)/2;
        if options.Verbose
            fprintf('[t=%d] obj=%f nev=%d lr=%.4f\n', ...
                     [t full(obj) nnz(sv) lr]);
        end
    end

    g(sv)=g(sv)+(out(sv)-data.Y(sv));

    if options.UseBias, g_b=sum(out(sv)-data.Y(sv))+oc/2; else g_b=0; end

    gn=-g;
    gn_b=-g_b;

    % goal condition: data based
    if mod(t-1,options.CgStopIter)==(options.CgStopIter-1)
        switch options.CgStopType
            case 1 % stability
                yfx_unlabeled=sign(Kalpha(unlabeled)+b);
                diff_unlabeled=1-nnz(yfx_unlabeled==yfx_unlabeled_prev)/u;
                if diff_unlabeled<options.CgStopParam, break, end
                yfx_unlabeled_prev=yfx_unlabeled;
            case 2 % validation
                valerr=1-nnz(sign(data.Kv*alpha+b)==data.Yv)/v;
                if (valerr>(valerr_prev-options.CgStopParam)), break, end
                valerr_prev=valerr;
            case 3 % stability & validation
                valerr=1-nnz(sign(data.Kv*alpha+b)==data.Yv)/v;
                yfx_unlabeled=sign(Kalpha(unlabeled)+b);
                diff_unlabeled=1-nnz(yfx_unlabeled==yfx_unlabeled_prev)/u;
                if (valerr>(valerr_prev-options.CgStopParam(1))) && ...
                    diff_unlabeled<options.CgStopParam(2), break, end
                valerr_prev=valerr;
                yfx_unlabeled_prev=yfx_unlabeled;
            case 4 % gradient norm
                ng=sqrt(sum(gn.^2)+gn_b^2);
                if ng/ng0<options.CgStopParam, break, end
            case 5 % preconditioned gradient norm
                ngp=sqrt(sum(go.^2)+go_b^2);
                if ngp/ngp0<options.CgStopParam, break, end
            case 6 % mixed gradient norm
                ngm=sqrt(Kgo'*go+go_b^2);
                if ngm/ngm0<options.CgStopParam, break, end
            case 7 % relative objective function decrease
                if (obj_prev-obj)<options.CgStopParam, break, end
                obj_prev=obj;
        end
    end

    Kgn=data.K*gn; % multiply by the preconditioner

    % debug
    if options.CgStopType==-1
        ng=sqrt(sum(Kgn.^2)+gn_b^2);
        ngp=sqrt(sum(gn.^2)+gn_b^2);
        ngm=sqrt(Kgn'*gn+gn_b^2);
        stats(t+1,:)=[t,obj/obj0,ng/ng0,ngp/ngp0,ngm/ngm0,alpha',b];
    end

    % Polack-Ribiere update with automatic restart
    be=max(0,(Kgn'*(gn-go)+gn_b*(gn_b-go_b))/(Kgo'*go+go_b^2));

    d=be*d+gn;
    d_b=be*d_b+gn_b;

    Kstep=be*Kstep+Kgn;

    go=gn;
    go_b=gn_b;
    Kgo=Kgn;
end
sec=toc;
if nargout>4, lsiters=lsiters(1:t);
    if nargout>5 && ~isempty(stats), stats=stats(1:t,:); end
end

function [lr,Kalpha,lsi,LKalpha] = linesearch(data,labeled,step,step_b, ...
                                              Kalpha,b,gamma_A,gamma_I, ...
                                              Kstep,LKalpha,hinge,oc)
% {linesearch} does a line search in direction step.

act=step~=0; % the set of points for which alpha change (active)
if isempty(Kstep)
    Kstep=data.K(:,act)*step(act);
end

% precomputations
stepKstep=step(act)'*Kstep(act);
stepKalpha=step(act)'*Kalpha(act);
if gamma_I~=0
    KstepL=Kstep'*data.L;
    KstepLKalpha=KstepL*Kalpha;
    KstepLKstep=KstepL*Kstep;
end

out=Kalpha(labeled)+b;
outstep=Kstep(labeled)+step_b;
out_minus_Y=out-data.Y(labeled);

% breakpoints
if hinge
    sv=(1-out.*data.Y(labeled))>0;
    deltas=-out_minus_Y./outstep;
    deltas(deltas<0)=0;
    [deltas,deltas_map]=sort(deltas);
    lab=length(deltas);
    i=find(deltas>0,1);
    lsi=1;
else
    sv=true(length(out),1);
    lsi=1;
end

% intercepts
if gamma_I~=0
    left=outstep(sv)'*out_minus_Y(sv)+gamma_A*stepKalpha+...
         gamma_I*KstepLKalpha + (oc/2)*step_b;
    right=outstep(sv)'*outstep(sv)+gamma_A*stepKstep+...
          gamma_I*KstepLKstep;
else
    left=outstep(sv)'*out_minus_Y(sv)+gamma_A*stepKalpha + (oc/2)*step_b;
    right=outstep(sv)'*outstep(sv)+gamma_A*stepKstep;
end

% first minimum
zcross=-left/right;
if right<=0, zcross=0; end

if hinge && ~isempty(i)
    if sv(deltas_map(i))==true
        if 0<=zcross && zcross<deltas(i), not_got_it=0;
        else not_got_it=1; end
    else
        if 0<=zcross && zcross<=deltas(i), not_got_it=0;
        else not_got_it=1; end
    end

    while not_got_it
        % updating support vectors
        j=-2*sv(deltas_map(i))+1;
        sv(deltas_map(i))=~sv(deltas_map(i));

        % updating intercepts
        left=left+j*outstep(deltas_map(i))*out_minus_Y(deltas_map(i));
        right=right+j*outstep(deltas_map(i))*outstep(deltas_map(i));

        % computing minimum
        zcross=-left/right;

        % goal conditions
        if i==lab, break, end
        if sv(deltas_map(i))==true
            if sv(deltas_map(i+1))==true
                if deltas(i)<zcross && zcross<deltas(i+1),
                    not_got_it=0; end
            else
                if deltas(i)<zcross && zcross<=deltas(i+1),
                    not_got_it=0; end
            end
        else
            if sv(deltas_map(i+1))==true
                if deltas(i)<=zcross && zcross<deltas(i+1),
                    not_got_it=0; end
            else
                if deltas(i)<=zcross && zcross<=deltas(i+1),
                    not_got_it=0; end
            end
        end

        i=i+1;
        lsi=lsi+1;
    end
end
lr=zcross;

if lr<0, lr=0; return, end % converged
Kalpha=Kalpha+lr*Kstep;
if nargout>3 && gamma_I~=0, LKalpha=LKalpha+lr*KstepL'; end

author: viggin

repo: domain-adaptation-toolbox


text:
% SB2_FULLSTATISTICS  Compute all relevant statistics in full for SPARSEBAYES
%
% [SIGMA,MU,S_IN,Q_IN,S_OUT,Q_OUT,FACTOR,LOGML,GAMMA,...
%	BETABASIS_PHI,BETA] = ...
%  SB2_FULLSTATISTICS(LIKELIHOOD,BASIS,PHI,TARGETS,USED,ALPHA,BETA,...
%       MU,BASIS_PHI,BASIS_TARGETS,OPTIONS)
%
% OUTPUT ARGUMENTS:
%
%	SIGMA			Posterior covariance matrix for relevant bases
%	MU				Posterior mean
%	S_IN			S-factors for in-model (relevant) basis vectors
%	Q_IN			Q-factors for in-model (relevant) basis vectors
%	S_OUT			S-factors for all basis vectors
%	Q_OUT			Q-factors for all basis vectors
%	FACTOR			Q^2-S relevance factors
%	LOGML			Log-marginal-likelihood
%	GAMMA			"Well-determinedness" factors
%	BETABASIS_PHI	Cached value of BASIS'*B*PHI matrix
%	BETA			Inverse noise variance (vector of beta
%					approximations in non-Gaussian case)
%
% INPUT ARGUMENTS:
%
%	LIKELIHOOD		LIKELIHOOD structure
%	BASIS			Full NxM basis matrix
%	PHI				Current relevant basis matrix
%	TARGETS			N-vector with target output values
%	USED			Relevant basis vector indices
%	ALPHA			Hyperparameters
%	BETA			Inverse noise variance (ignored in non-Gauss case)
%	MU				Current posterior mean (for non-Gauss)
%	BASIS_PHI		Cached value of BASIS'*PHI (for Gauss only)
%	BASIS_TARGETS	Cached value of BASIS'*TARGETS (for Gauss only)
%	OPTIONS			Standard OPTIONS structure (see SB2_USEROPTIONS)
%
% NOTES:
%
% This function computes the posterior, and other, statistics for the
% SPARSEBAYES algorithm in "long-hand" fashion. i.e. it does not use the
% more effecient "iterative" updates.
%
% It is required on every iteration (unless approximating) in the
% non-Gaussian case, and on iterations in the Gaussian case where the noise
% estimate (BETA) is updated.
%
% This function is intended for internal use by SPARSEBAYES only.
%

%
% Copyright 2009, Vector Anomaly Ltd
%
% This file is part of the SPARSEBAYES library for Matlab (V2.0).
%
% SPARSEBAYES 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 2 of the License, or (at your option)
% any later version.
%
% SPARSEBAYES 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 SPARSEBAYES in the accompanying file "licence.txt"; if not, write to
% the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
% MA 02110-1301 USA
%
% Contact the author: m a i l [at] m i k e t i p p i n g . c o m
%
function [SIGMA,Mu,S_in,Q_in,S_out,Q_out,Factor,logML,Gamma,...
	  betaBASIS_PHI,beta] = ...
    SB2_FullStatistics(LIKELIHOOD,BASIS,PHI,Targets,Used,Alpha,beta,...
		       Mu,BASIS_PHI,BASIS_Targets,OPTIONS)

% Mu is only required for non-Gauss
% BASIS_PHI and BASIS_Targets are only required for Gauss

% beta (a vector for non-Gauss) is returned only for non-Gauss

MAX_POSTMODE_ITS	= 25; % More than enough

[N M_full]	= size(BASIS);
M			= size(PHI,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% COMPUTE POSTERIOR

%
% Compute full relevant statistics
%
if LIKELIHOOD.InUse==LIKELIHOOD.Gaussian
  %
  % Posterior is analytic: use linear algebra here
  %
  % Compute posterior covariance SIGMA (inverse Hessian)
  % via Cholesky decomposition
  %
  U		= chol(PHI'*PHI*beta + diag(Alpha));
  Ui	= inv(U);
  SIGMA	= Ui * Ui';
  % Posterior mean Mu
  Mu	= (SIGMA * (PHI'*Targets)) * beta;
  % Data error and likelihood
  y				= PHI * Mu;
  e				= (Targets - y);
  ED			= e'*e;
  %
  dataLikely	= (N*log(beta) - beta*ED)/2;
  %
else
  %
  % Posterior must be approximated: find the posterior mode as the basis
  % for the Laplace approximation
  %
  [Mu U beta dataLikely] = ...
      SB2_PosteriorMode(LIKELIHOOD,PHI,Targets,Alpha,Mu, ...
						MAX_POSTMODE_ITS,OPTIONS);
  % Compute covariance approximation
  %
  Ui	= inv(U);
  SIGMA	= Ui * Ui';
  % Compute posterior-mean-based outputs and errors
  %
  switch LIKELIHOOD.InUse
   case LIKELIHOOD.Bernoulli,
    y	= SB2_Sigmoid(PHI * Mu);
   case LIKELIHOOD.Poisson
    y	= exp(PHI*Mu);
  end
  e		= (Targets-y);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% COMPUTE LOG MARGINAL LIKELIHOOD

logdetHOver2	= sum(log(diag(U)));
logML			= dataLikely - (Mu.^2)'*Alpha/2 + ...
    sum(log(Alpha))/2 - logdetHOver2;
% Well-determinedness factors
DiagC	= sum(Ui.^2,2);
Gamma	= 1 - Alpha.*DiagC;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% COMPUTE THE Q & S VALUES

%
% Q: "quality" factor - related to how well the basis function contributes
% to reducing the error
%
% S: "sparsity factor" - related to how orthogonal a given basis function
% is to the currently used set of basis functions
%
if LIKELIHOOD.InUse==LIKELIHOOD.Gaussian
  %
  % Gaussian case simple: beta a scalar
  %
  betaBASIS_PHI	= beta*BASIS_PHI;
  %
  % The S_in calculation exploits the fact that BASIS is normalised,
  %  i.e. sum(BASIS.^2,1)==ones(1,M_full)
  %
  S_in		= beta - sum((betaBASIS_PHI*Ui).^2,2);
  Q_in		= beta*(BASIS_Targets - BASIS_PHI*Mu);
else
  %
  % Non-Gaussian case: beta an N-vector
  %
  betaBASIS_PHI	= BASIS'*(PHI.* (beta*ones(1,M)));
  S_in			= (beta'*(BASIS.^2))' - sum((betaBASIS_PHI*Ui).^2,2);
  Q_in			= BASIS'*e;
end
%
S_out		= S_in;
Q_out		= Q_in;
%
% S,Q with that basis excluded: equations (23)
%
S_out(Used)	= (Alpha .* S_in(Used)) ./ (Alpha - S_in(Used));
Q_out(Used)	= (Alpha .* Q_in(Used)) ./ (Alpha - S_in(Used));
%
% Pre-compute the "relevance factor" for ongoing convenience
%
Factor		= (Q_out.*Q_out - S_out);

repo: REST


text:
function y = sinh(x)
%SINH         Implements  sinh(x)  for intervals
%
%   y = sinh(x)
%
%interval standard function implementation
%

% written  12/30/98     S.M. Rump
% modified 08/31/99     S.M. Rump  complex allowed, sparse input,
%                                  major revision, improved accuracy
% modified 09/02/00     S.M. Rump  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
% modified 12/04/05     S.M. Rump  extreme values for approximate part
% modified 09/06/07     S.M. Rump  approximate std fcts removed
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/13/12     S.M. Rump  INTLAB_INTVAL_STDFCTS
%

  if issparse(x)
    [ix,jx,sx] = find(x);
    [m,n] = size(x);
    y = sparse(ix,jx,sinh(full(sx)),m,n);
    return
  end

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if x.complex
    y = ( exp(x) - exp(-x) ) / 2;
    setround(rndold)
    return
  end

  y = x;

  xinf = x.inf(:);
  xsup = x.sup(:);

  IndexInfPos = ( xinf>=0 );
  len1 = sum(IndexInfPos);
  IndexSupNeg = ( xsup<=0 );
  len2 = sum(IndexSupNeg);

  Y = sinhpos( [ xinf(IndexInfPos) ; -xsup(IndexSupNeg) ] , -1 );
  y.inf(IndexInfPos) = Y(1:len1);
  y.sup(IndexSupNeg) = -Y( len1+1 : end );

  IndexInfNeg = ~IndexInfPos;
  len1 = sum(IndexInfNeg);
  IndexSupPos = ~IndexSupNeg;
  len2 = sum(IndexSupPos);

  Y = sinhpos( [ -xinf(IndexInfNeg) ; xsup(IndexSupPos) ] , 1 );
  y.inf(IndexInfNeg) = -Y(1:len1);
  y.sup(IndexSupPos) = Y( len1+1 : len1+len2 );

  setround(rndold)

function y = sinhpos(x,rnd)
% rigorous sinh(x) for nonnegative double vector x with
% rounding corresponding to rnd

  INTLAB_STDFCTS_SINH = getappdata(0,'INTLAB_STDFCTS_SINH');

  y = x;

  % small input
  index = ( x<8 );
  if any(index)
    setround(0)
    xx = x(index);
    [f,e] = log2(xx);
    bits = 14 + min(e,0);
    xs = pow2( floor(f.*2.^bits) , e-bits );  % max. 14 bits of mantissa,
                                              % no bit below 2^-14
    d = xx - xs;                              % 0 <= d < 2^-14*x < 2^-11
    sinhxs = sinh(xs);                        % round to nearest

    % bounds for exp(xs)
    INTLAB_STDFCTS_EXP = getappdata(0,'INTLAB_STDFCTS_EXP');
    setround(-1)
    expinf = exp(xs)*(1-INTLAB_STDFCTS_EXP.EPS);
    setround(1)
    expsup = exp(xs)*(1+INTLAB_STDFCTS_EXP.EPS);
    expinf(xs==0) = 1;                    % gives accuracy for small x
    expsup(xs==0) = 1;

    % use sinh(xx) = sinh(xs+d) = sinh(xs)*cosh(d) + cosh(xs)*sinh(d)
    setround(rnd)
    corr = INTLAB_STDFCTS_SINH.EPS;
    if rnd==-1
      % coshxs <= cosh(xs)
      coshxs = 0.5*( expinf + 1./expsup );  % rounded downwards
      dd = d.*d;
      % 0 <= err(sinhd) <= sinh(d)*d^4/4! < 6e-8*d^3/6
      % 0 <= err(coshd) <= sinh(d)*d^5/5! < 4.8e-8*d^4/24
      y(index) = sinhxs + ...
          ( ( sinhxs .* ( -corr + (1-corr)*dd/2.*( 1 + dd/12 ) ) + ...
              coshxs.*d.*dd/6 ) + ...
            coshxs.*d ...
          );
    else
      % coshxs >= cosh(xs)
      coshxs = 0.5*( expsup + 1./expinf );  % correctly rounded according to rnd
      dd = d.*d;
      y(index) = sinhxs + ...
          ( ( sinhxs .* ( corr + (1+corr)*dd/2.*( 1 + dd/12.*( 1+4.8e-8 ) ) ) + ...
              coshxs.*d.*dd/6.*( 1+6e-8 ) ) + ...
            coshxs.*d ...
          );
    end
  end

  % medium input
  index = ( ~index ) & ( x<709 );         % 8 <= x < 709
  if any(index)
    xx = x(index);
    exprnd = exp_rnd(xx,rnd);
    setround(rnd)
    y(index) = 0.5 * ( exprnd + (-1)./exprnd );
  end

  % large input
  index = ( x>=709 );
  if any(index)
    INTLAB_STDFCTS_E = getappdata(0,'INTLAB_STDFCTS_E');
    exprnd = exp_rnd( x(index)-1 , rnd )/2;
    setround(rnd)
    if rnd==-1
      y(index) = exprnd .* INTLAB_STDFCTS_E.INF;
    else
      y(index) = exprnd .* INTLAB_STDFCTS_E.SUP;
    end
  end

repo: OpenMAS


text:
close all; clear all; clc

% include tclab.m for initialization
tclab;

% Write data to the digital pin
disp('LED on and off repeatedly for 5 seconds')
for i = 0:5
    led(1)
    pause(0.5);
    led(0)
    pause(0.5);
end

% Change the brightness from maximum to minimum using the
%  digital pin's PWM duty cycle.
disp('LED dimming from bright to off over 10 seconds')
for brightness = 1:-0.1:0
    led(brightness);
    pause(1.0);
end

disp('LED Test Complete')

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function vbmc_iterplot(vp,gp,optimState,stats,elbo)
%VBMC_ITERPLOT Plot current iteration of the VBMC algorithm.

D = vp.D;
iter = optimState.iter;
fontsize = 14;

if D == 1
    hold off;
    gplite_plot(gp);
    hold on;
    xlims = xlim;
    xx = linspace(xlims(1),xlims(2),1e3)';
    yy = vbmc_pdf(vp,xx,false,true);
    hold on;
    plot(xx,yy+elbo,':');
    drawnow;

else
    if ~isempty(vp)
        Xrnd = vbmc_rnd(vp,1e5,1,1);
    else
        Xrnd = gp.X;
    end
    X_train = gp.X;

    if iter == 1
        idx_new = true(size(X_train,1),1);
    else
        X_trainold = stats.gp(iter-1).X;
        idx_new = false(size(X_train,1),1);
        [~,idx_diff] = setdiff(X_train,X_trainold,'rows');
        idx_new(idx_diff) = true;
    end
    idx_old = ~idx_new;

    if ~isempty(vp.trinfo); X_train = warpvars_vbmc(X_train,'inv',vp.trinfo); end

    Pdelta = optimState.PUB_orig - optimState.PLB_orig;
    X_min = min(X_train,[],1) - Pdelta*0.1;
    X_max = max(X_train,[],1) + Pdelta*0.1;
    bounds = [max(min(optimState.PLB_orig,X_min),optimState.LB_orig); ...
        min(max(optimState.PUB_orig,X_max),optimState.UB_orig)];

    try
        for i = 1:D; names{i} = ['x_{' num2str(i) '}']; end
        [~,ax] = cornerplot(Xrnd,names,[],bounds);
        for i = 1:D-1
            for j = i+1:D
                axes(ax(j,i));  hold on;
                if any(idx_old)
                    scatter(X_train(idx_old,i),X_train(idx_old,j),'ok');
                end
                if any(idx_new)
                    scatter(X_train(idx_new,i),X_train(idx_new,j),'or','MarkerFaceColor','r');
                end
            end
        end

        h = axes(gcf,'Position',[0 0 1 1]);
        set(h,'Color','none','box','off','XTick',[],'YTick',[],'Units','normalized','Xcolor','none','Ycolor','none');
        text(0.9,0.9,['VBMC (iteration ' num2str(iter) ')'],'FontSize',fontsize,'HorizontalAlignment','right');

        drawnow;
    catch
        % pause
    end
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function Chi2 = chistart (D,L,a,ncands,factor)
%CHISTART: Computes the initial size of the search ellipsoid
%
% This routine computes or approximates the initial size of the search
% ellipsoid. If the requested number of candidates is not more than the
% dimension + 1, this is done by computing the squared distances of partially
% conditionally rounded float vectors to the float vector in the metric of the
% covariance matrix. Otherwise an approximation is used.
%
% Input arguments
%    L,D   : LtDL-decomposition of the variance-covariance matrix of
%            the float ambiguities (preferably decorrelated)
%    a     : float ambiguites (preferably decorrelated)
%    ncands: Requested number of candidates (default = 2)
%    factor: Multiplication factor for the volume of the resulting
%            search ellipsoid (default = 1.5)
%
% Output arguments:
%    Chi2  : Size of the search ellipsoid

% ----------------------------------------------------------------------
% File.....: chistart.m
% Date.....: 19-MAY-1999
% Modified.: 05-MAR-2001, by P. Joosten
% Author...: Peter Joosten
%            Mathematical Geodesy and Positioning
%            Delft University of Technology
% ----------------------------------------------------------------------

% ------------------
% --- Initialize ---
% ------------------

if nargin < 4; ncands = 2  ; end;
if nargin < 5; factor = 1.5; end;

n = max(size(a));

% ----------------------------------------------------------------------
% --- Computation depends on the number of candidates to be computed ---
% ----------------------------------------------------------------------

if ncands <= n+1;

    % --------------------------------------------------------
    % --- Computation based on the bootstrapping estimator ---
    % --------------------------------------------------------

    Chi = [];

    for k = n:-1:0;

        afloat = a;
        afixed = a;

        for i = n:-1:1;

            dw = 0;
            for j = n:-1:i;
                dw = dw + L(j,i) * (afloat(j) - afixed(j));
            end;

            afloat(i) = afloat(i) - dw;
            if (i ~= k);
                afixed(i) = round (afloat(i));
            else
                if isequal (afloat(i),afixed(i));
                    afixed(i) = round(afixed(i) + 1);
                else
                    afixed(i) = round (afloat(i) + sign (afloat(i) - afixed(i)));
                end;
            end;

        end;

        Chi = [Chi (a-afixed)' * (L'*diag(D)*L)^(-1) * (a-afixed)];

    end;

    % ---------------------------------------------------------------
    % --- Sort the results, and return the appropriate number     ---
    % --- Add an "eps", to make sure there is no boundary problem ---
    % ---------------------------------------------------------------

    Chi  = sort(Chi);
    Chi2 = Chi(ncands) + 1d-6;

else

    % -----------------------------------------------------
    % An approximation for the squared norm is computed ---
    % -----------------------------------------------------

    Linv = (L)^(-1);
    Dinv = 1./D;

    Vn   = (2/n) * (pi ^ (n/2) / gamma(n/2));
    Chi2 = factor * (ncands / sqrt((prod(1 ./ Dinv)) * Vn)) ^ (2/n);

end;

% ----------------------------------------------------------------------
% End of routine: chistart
% ----------------------------------------------------------------------

author: goGPS-Project


text:
% evalreprerror ... computes and plots the final error statistics
%
% cam = evalreprerror(cam,config)
% cam, config ... see the main GOCAL script
%
% $Id: evalreprerror.m,v 2.0 2003/06/19 12:07:03 svoboda Exp $

function cam = evalreprerror(cam,config)

Octave = exist('OCTAVE_VERSION', 'builtin') ~= 0;

disp('2D reprojection error')
disp(sprintf('All points: mean  %2.2f pixels, std is %2.2f',mean([cam.err2d]), std([cam.err2d])'));
% disp(sprintf('Inliers:    mean  %2.2f pixels, std is %2.2f',mean([cam.inerr2d]), std([cam.inerr2d])'));
if mean([cam.err2d])>1.5 | std([cam.err2d])>1.5
	disp('***************************************************')
	disp('W A R N I N G: the reprojection error is relatively high !')
end

%%%
% evaluate the reprojection error for each camera separately to detect possible problems
%%%
for i=1:size(config.cal.cams2use,2),
	cam(i).mean2Derr = mean(cam(i).err2d);
	cam(i).std2Derr  = std(cam(i).err2d);
end
% sort the values and print them to the 2D graphs

if ~Octave,
figure(30),
clf
plot(config.cal.cams2use,[cam.mean2Derr],'bd'),
hold on, grid on,
plot(config.cal.cams2use,[cam.mean2Derr],'b-'),
plot(config.cal.cams2use,[cam.std2Derr],'rd')
plot(config.cal.cams2use,[cam.std2Derr],'r-')
xlabel('Id of the camera')
title('2D error: mean (blue), std (red)')
ylabel('pixels')

figure(31)
clf
bar(config.cal.cams2use,[cam.mean2Derr;cam.std2Derr]',1.5)
grid on
xlabel('Id of the camera')
title('2D error: mean (blue), std (red)')
ylabel('pixels')

figure(31),
eval(['print -depsc ', config.paths.data, 'reprerrors.eps'])

figure(4),
eval(['print -depsc ', config.paths.data, 'reconstructedsetup.eps'])
end

Ret = 1;
return

author: strawlab

repo: MultiCamSelfCal


text:
function [v,x] = max(S1F,varargin)

x = linspace(0,2*pi);
y = S1F.eval(x);
[v,ind] = max(y);

x = x(ind);

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [PC, powers] = powercentersPD(T, E, wts)
% function [PC, powers] = powercentersPD(T, E, wts)
%
% T: triangulation
% E: set of points
% wts: weights of points in E
%
% The output array PC contains the power centers of the triangles in T.
% That is, row i of PC contains the point that is equidistant from each
% vertex of the triangle specified by row i of T, with respect to the power
% of each vertex. The output array powers contains the power of each power
% center.

[m, n] = size(T);
PC = zeros(m,n-1);
powers = zeros(m,1);

for i=1:m
    tr = E(T(i,:),:);
    wt = wts(T(i,:));
    p = tr(1,:);
    Rp = repmat(p, n-1, 1);
    Pts = tr(2:n,:);
    Ac = 2*(Pts - Rp);

    Rw1 = repmat(wt(1), n-1, 1);
    Wts = wt(2:n);
    Sp1 = repmat(sum(p.^2), n-1, 1);
    SPts = sum(Pts.^2, 2);
    Bc = Rw1 - Wts - Sp1 + SPts;

    pc = Ac \ Bc;

    PC(i,:) = pc;
    powers(i,1) = norm(pc - p')^2 - wt(1);
end

author: optimaltransport


text:
function [val] = filetype_check_extension(filename, ext)

% FILETYPE_CHECK_EXTENSION helper function to determine the file type
% by performing a case insensitive string comparison of the extension.

% Copyright (C) 2003-2012 Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are for remembering the type on subsequent calls with the same input arguments
persistent previous_argin previous_argout

current_argin = {filename, ext};
if isequal(current_argin, previous_argin)
  % don't do the detection again, but return the previous value from cache
  val = previous_argout;
  return
end

if iscell(filename)
  % compare the extension of multiple files
  val = zeros(size(filename));
  for i=1:length(filename)
    val(i) = filetype_check_extension(filename{i}, ext);
  end
elseif iscell(ext)
  val = zeros(size(ext));
  for i=1:length(ext)
    val(i) = filetype_check_extension(filename, ext{i});
  end
else
  % compare the extension of a single file
  if numel(filename)<numel(ext)
    val = false;
  else
    val = strcmpi(filename((end-length(ext)+1):end), ext);
  end
end

% remember the current input and output arguments, so that they can be
% reused on a subsequent call in case the same input argument is given
current_argout  = val;
previous_argin  = current_argin;
previous_argout = current_argout;

return % main()

author: spm


text:
function s = spm_existfile(filename)
% Check if a file exists on disk - a compiled routine
% FORMAT S = SPM_EXISTFILE(FILENAME)
% FILENAME - filename (can also be a relative or full pathname to a file)
% S        - logical scalar, true if the file exists and false otherwise
%_______________________________________________________________________
%
% This compiled routine is equivalent to:
% >> s = exist(filename,'file') == 2;
% and was written for speed purposes. The differences in behaviour are:
%  * spm_existfile returns true for directory names
%  * spm_existfile does not look in MATLAB's search path
%  * spm_existfile returns false for an existing file that does not have
%    read permission
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_existfile.m 3100 2009-05-06 19:00:39Z guillaume $

%-This is merely the help file for the compiled routine
%error('spm_existfile.c not compiled - see Makefile')
persistent runonce
if isempty(runonce)
    warning('spm_existfile is not compiled for your platform.');
    runonce = 1;
end

s = exist(filename,'file') > 0;

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function DVH_out = dvh_Rectum_shift_1cm(planC,structNum,doseNum)
%function DVH_out = dvh_Rectum_shift_1cm(planC,structNum,doseNum)
%
%APA, 11/09/2010
%
% 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/>.

%Number of Fractions and Trials
numFractions = 30;
numTrials    = 100;

% Systematic Shifts in cm (Inter-treatment)
XdispSys_Std = 0;
YdispSys_Std = 0;
ZdispSys_Std = 0;

% Random Shifts in cm (inter-fraction)
XdispRnd_Std = 1;
YdispRnd_Std = 1;
ZdispRnd_Std = 1;

indexS = planC{end};

%Distances are in cm
XpdfMean = 0;

YpdfMean = 0;

ZpdfMean = 0;

%Get scan associated with doseNum
scanNum = getAssociatedScan(planC{indexS.dose}(doseNum).assocScanUID, planC);

if isempty(scanNum) %Assume dose is associated with this scan
    scanNum = getStructureAssociatedScan(structNum,planC);
end

%Get reference transformation matrix for doseNum
if ~isempty(scanNum) && isempty(planC{indexS.dose}(doseNum).transM)
    referenceTransM = planC{indexS.scan}(scanNum).transM;
else
    referenceTransM = planC{indexS.dose}(doseNum).transM;
end

%Store the reference doseArray
referenceDoseArray = planC{indexS.dose}(doseNum).doseArray;

%Try and get a binWidth from stateS.  If it doesnt exist, get it from
%the CERROptions file (allows this function to be called outside CERR)
global stateS;
if ~isempty(stateS) && isfield(stateS, 'optS') && isfield(stateS.optS, 'DVHBinWidth') && ~isempty(stateS.optS.DVHBinWidth)
    binWidth = stateS.optS.DVHBinWidth;
else
    optS = CERROptions;
    binWidth = optS.DVHBinWidth;
end

%Compute DVH at given dose
[dosesCurrentV, volsV] = getDVH(structNum, doseNum, planC);

%Divide doseArray into fractions
planC{indexS.dose}(doseNum).doseArray = planC{indexS.dose}(doseNum).doseArray / numFractions;

%Get systematic errors
deltaX_systematic = randn(1, numTrials) * XdispSys_Std + XpdfMean;
deltaY_systematic = randn(1, numTrials) * YdispSys_Std + YpdfMean;
deltaZ_systematic = randn(1, numTrials) * ZdispSys_Std + ZpdfMean;

hWait = waitbar(0,'Computing plan robustness...');

if isempty(referenceTransM)
    referenceTransM_matrix = eye(4);
else
    referenceTransM_matrix = referenceTransM;
end

%Obtain DVH calculation points in blocks
optS    = planC{indexS.CERROptions};

ROIImageSize = [planC{indexS.scan}(scanNum).scanInfo(1).sizeOfDimension1  planC{indexS.scan}(scanNum).scanInfo(1).sizeOfDimension2];

deltaY = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;

%Get raster segments for structure.
[segmentsM, planC, isError] = getRasterSegments(structNum, planC);

if isempty(segmentsM)
    isError = 1;
end
numSegs = size(segmentsM,1);

%Relative sampling of ROI voxels in this place, compared to CT spacing.
%Set when rasterSegments are generated (usually on import).
sampleRate = optS.ROISampleRate;

%Sample the rows
indFullV =  1 : numSegs;
if sampleRate ~= 1
    rV = 1 : length(indFullV);
    rV([rem(rV+sampleRate-1,sampleRate)~=0]) = [];
    indFullV = rV;
end

%Block process to avoid swamping on large structures
if isfield(optS, 'DVHBlockSize') & ~isempty(optS.DVHBlockSize)
    DVHBlockSize = optS.DVHBlockSize;
else
    DVHBlockSize = 5000;
end

blocks = ceil(length(indFullV)/DVHBlockSize);

start = 1;

for b = 1 : blocks

  %Build the interpolation points matrix

  dummy = zeros(1,DVHBlockSize * ROIImageSize(1));
  x1V = dummy;
  y1V = dummy;
  z1V = dummy;
  volsSectionV =  dummy;

  if start+DVHBlockSize > length(indFullV)
    stop = length(indFullV);
  else
    stop = start + DVHBlockSize - 1;
  end

  indV = indFullV(start:stop);

  mark = 1;
  for i = indV

    tmpV = segmentsM(i,1:10);
    delta = tmpV(5) * sampleRate;
    xV = tmpV(3): delta : tmpV(4);
    len = length(xV);
    rangeV = ones(1,len);
    yV = tmpV(2) * rangeV;
    zV = tmpV(1) * rangeV;
    sliceThickness = tmpV(10);
    %v = delta^2 * sliceThickness;
    v = delta * (deltaY*sampleRate) * sliceThickness;
    x1V(mark : mark + len - 1) = xV;
    y1V(mark : mark + len - 1) = yV;
    z1V(mark : mark + len - 1) = zV;
    volsSectionV(mark : mark + len - 1) = v;
    mark = mark + len;

  end

  %cut unused matrix elements
  x1V = x1V(1:mark-1);
  y1V = y1V(1:mark-1);
  z1V = z1V(1:mark-1);
  volsSectionV = volsSectionV(1:mark-1);

  %Get transformation matrices for both dose and structure.
  transMDose    = getTransM('dose', doseNum, planC);
  transMStruct  = getTransM('struct', structNum, planC);

  %Forward transform the structure's coordinates.
  if ~isempty(transMStruct)
      [x1V, y1V, z1V] = applyTransM(transMStruct, x1V, y1V, z1V);
  end

  dvhCalcPtsC{b} = [x1V', y1V', z1V'];

  %Interpolate.
%   [dosesSectionV] = getDoseAt(doseNum, x1V, y1V, z1V, planC);

%   dosesV = [dosesV, dosesSectionV];
%   volsV  = [volsV, volsSectionV];

  start = stop + 1;

end

%Loop over to generate DVH
for iTrial = 1:numTrials

    deltaX = deltaX_systematic(iTrial) + randn(1, numFractions) * XdispRnd_Std + XpdfMean;
    deltaY = deltaY_systematic(iTrial) + randn(1, numFractions) * YdispRnd_Std + YpdfMean;
    deltaZ = deltaZ_systematic(iTrial) + randn(1, numFractions) * ZdispRnd_Std + ZpdfMean;

    tmpDoseV = zeros(1,length(dosesCurrentV));

    for iFraction = 1:numFractions

        waitbar(((iTrial-1)*numFractions + iFraction)/(numTrials*numFractions),hWait)

        transM = referenceTransM_matrix;
        transM(1:3,4) = transM(1:3,4) + [deltaX(iFraction); deltaY(iFraction); deltaZ(iFraction)];

        %Apply the new transM to dose
        planC{indexS.dose}(doseNum).transM = transM;

        %Get doses and volumes of points in structure.
        %[dosesV, volsV] = getDVH(structNum, doseNum, planC);

        volsV  = [];
        dosesV = [];

        for b = 1 : blocks
            x1V = dvhCalcPtsC{b}(:,1);
            y1V = dvhCalcPtsC{b}(:,2);
            z1V = dvhCalcPtsC{b}(:,3);

            %Back transform the coordinates into the doses' coordinate system.
            [x1V, y1V, z1V] = applyTransM(inv(transM), x1V, y1V, z1V);

            %Interpolate.
            [dosesSectionV] = getDoseAt(doseNum, x1V, y1V, z1V, planC);

            dosesV = [dosesV, dosesSectionV];
            volsV  = [volsV, volsSectionV];
        end

        tmpDoseV = tmpDoseV + dosesV;

    end

    DVHm(:,iTrial) = single(tmpDoseV(:));

    %Compute Mean and Std for binned dose
    %for iTrial=1:numTrials
        doseTrialV = DVHm(:,iTrial);
        [doseBinsV, volsHistV] = doseHist(doseTrialV, volsV, binWidth);
        doseBinsTmpV = doseBinsV;
        volsHistTmpV = volsHistV;
        if iTrial > 1
            AddToCurrent = ~ismember(doseBins{1},doseBinsV);
            AddToPrevious = ~ismember(doseBinsV,doseBins{1});
            doseBinsV = [doseBinsV doseBins{1}(AddToCurrent)];
            volsHistV = [volsHistV volsHist{1}(AddToCurrent)];
            for jTrial = 1:iTrial-1
                doseBins{jTrial} = [doseBins{jTrial} doseBinsTmpV(AddToPrevious)];
                volsHist{jTrial} = [volsHist{jTrial} volsHistTmpV(AddToPrevious)];
                [doseBins{jTrial},indSort] = sort(doseBins{jTrial});
                volsHist{jTrial} = volsHist{jTrial}(indSort);
            end
        end
        doseBins{iTrial} = doseBinsV;
        volsHist{iTrial} = volsHistV;
        [doseBins{iTrial},indSort] = sort(doseBins{iTrial});
        volsHist{iTrial} = volsHist{iTrial}(indSort);
    %end

end

close(hWait)

%Compute Std of binned volumes
blockSize = 100;
for iBlock = 1:ceil(length(volsHist{1})/blockSize)
    if iBlock == ceil(length(volsHist{1})/blockSize)
        indicesV = (iBlock-1)*blockSize+1:length(volsHist{1});
    else
        indicesV = (iBlock-1)*blockSize+1:iBlock*blockSize;
    end
    volsHistStdM = [];
    for iTrial = 1:length(volsHist)
        volsHistStdM(:,iTrial) = volsHist{iTrial}(indicesV);
    end
    volsHistStdV(indicesV) = std(volsHistStdM,0,2)';
end

%Reassign reference transM and doseArray to doseNum
planC{indexS.dose}(doseNum).transM = referenceTransM;
planC{indexS.dose}(doseNum).doseArray = referenceDoseArray;

%Compute Mean and Std Dev in blocks
numVoxels = length(volsV);
DVHBlockSize = 5000;
blocks = ceil(numVoxels/DVHBlockSize);
start = 1;

for b = 1 : blocks

  if start+DVHBlockSize > numVoxels
    stop = numVoxels;
  else
    stop = start + DVHBlockSize - 1;
  end

  DVH_block = DVHm(start:stop,:);

  meanDoseV(start:stop) = mean(DVH_block,2);
  stdDoseV(start:stop)  = std(DVH_block,0,2);

  start = stop + 1;

end

%Histogram for the expectation of dose over all trials
[doseBinsV, volsHistV] = doseHist(meanDoseV, volsV, binWidth);

DVH_out = [doseBinsV(:)'; volsHistV(:)'];

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%%********************************************************
%% qops:  Fu = qops(pblk,w,f,options,u);
%%
%% options = 1, Fu(i) = <wi,fi>
%%         = 2, Fu(i) = 2*wi(1)*fi(1)-<wi,fi>
%%         = 3, Fui   = w(i)*fi
%%         = 4, Fui   = w(i)*fi, Fui(1) = -Fui(1).
%% options = 5, Fu = w [ f'*u ; ub + fb*alp ], where
%%              alp = (f'*u + u0)/(1+f0);
%% options = 6, compute Finv*u.
%%
%% Note F = w [f0  fb'; fb  I+ fb*fb'/(1+f0) ], where
%%      f0*f0 - fb*fb' = 1.
%%      Finv = (1/w) [f0  -fb'; -fb  I+ fb*fb'/(1+f0) ].
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%********************************************************

function  Fu = qops(pblk,w,f,options,u)

if (options >= 1 && options <= 4)
    Fu = mexqops(pblk{2},w,f,options);
elseif (options == 5)
    s = 1 + [0 cumsum(pblk{2})];
    idx1 = s(1:length(pblk{2}));
    inprod = mexqops(pblk{2},f,u,1);
    tmp = (u(idx1)+inprod)./(1+f(idx1));
    Fu  = u + mexqops(pblk{2},tmp,f,3);
    Fu(idx1) = inprod;
    Fu = mexqops(pblk{2},w,Fu,3);
elseif (options == 6)
    s = 1 + [0 cumsum(pblk{2})];
    idx1 = s(1:length(pblk{2}));
    gamprod = mexqops(pblk{2},f,u,2);
    tmp = (u(idx1)+gamprod)./(1+f(idx1));
    Fu  = u - mexqops(pblk{2},tmp,f,3);
    Fu(idx1) = gamprod;
    Fu = mexqops(pblk{2},1./w,Fu,3);
end
%%********************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [E_z, E_zz] = estep_compute_Z_distr(P, S_bar, V, RR, Tr, sigma_sq)
%[E_z, E_zz] = estep_compute_Z_distr(P, S_bar, V, RR, Tr, sigma_sq)

% Computes the distribution over Z given the current parameter estimates (see Eq 17-18)

K = size(V,1)/3;
[T, J] = size(P);
T = T/2;

Pc = P - Tr(:)*ones(1,J);

M_t = zeros(2*J, K);
P_hat_t = zeros(2*J, 1);

E_z = zeros(K, T);
E_zz = zeros(T*K, K);

invSigmaSq_p = eye(2*J)./sigma_sq;
for t=1:T,
   R_t = [RR(t,:); RR(t+T,:)];

   for kk = 1:K,
      M_t(1:J, kk) = (R_t(1,:)*V(1+(kk-1)*3:kk*3, :))';
      M_t(J+1:end, kk) = (R_t(2,:)*V(1+(kk-1)*3:kk*3, :))';
   end
   P_hat_t(1:J) = (R_t(1,:)*S_bar)';
   P_hat_t(J+1:end) = (R_t(2,:)*S_bar)';

   %beta_t = M_t' * inv(M_t*M_t' + sigma_sq*eye(2*J));                      % (Eq 16)
   % Can be computed much more efficiently using the matrix inversion lemma:
   AA = M_t./sigma_sq;
   beta_t = M_t'*(invSigmaSq_p - AA*inv(eye(K) + M_t'*M_t./sigma_sq)*AA');   % (Eq 16)

   E_z(:, t) = beta_t*([Pc(t, :) Pc(t+T, :)]' - P_hat_t);                    % (Eq 17)
   E_zz((t-1)*K+1:t*K,:) = eye(K) - beta_t*M_t + E_z(:,t)*E_z(:,t)';         % (Eq 18)
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function view = addROIrectMontage(view,sgn)
%
% view = addROIrectMontage(view,[sgn])
%
% If sgn~=0, adds user-specified rectangle to selected ROI in
% current slice. If sgn==0, removes the rectangle from the ROI.
%
% djh, 7/98
% ras, 9/04 -- updated for montage view
% If you change this function make parallel changes in:
%   all addROI*.m functions

% error if no current ROI
if view.selectedROI == 0
  myErrorDlg('No current ROI');
  return
end

if ~exist('sgn','var')
  sgn = 1;
end

% Get current ROI coords
curCoords = getCurROIcoords(view);

% Save prevSelpts for undo
view.prevCoords = curCoords;

% Get region from user.
% (rgn is short for region)
rgn = round(ginput(2));

% Make sure 2nd value is larger than the 1st.
for i=1:2
  if rgn(2,i) < rgn(1,i)
    rgn(:,i)=flipud(rgn(:,i));
  end
end

% Note: ginput hands them back in x, y order (1st col is x and
% 2nd col is y).  But we use them in the opposite order (y then
% x), so flip 'em.
rgn = fliplr(rgn);

% Get rgnSize
rgnSize = rgn(2,:) - rgn(1,:) + 1;

% Compute new coordinates
indices = 1:prod(rgnSize);
tmpCoords = indices2Coords(indices,rgnSize);
montageCoords = zeros(2,length(indices));
montageCoords(1,:) = tmpCoords(1,:) + rgn(1,1);
montageCoords(2,:) = tmpCoords(2,:) + rgn(1,2);
newCoords = montage2Coords(view,montageCoords);

% Do an (inverse) rotation if necessary
if (strcmp(view.viewType,'Flat'))
    % The '1' here signifies inverse rotation
	% (from screen to zero (internal) representation)
    newCoords = rotateCoords(view,newCoords,1);
end

% Convert coords to canonical frame of reference
newCoords = curOri2CanOri(view,newCoords);

% Merge/remove coordinates
if sgn
  coords = mergeCoords(curCoords,newCoords);
else
  coords = removeCoords(newCoords,curCoords);
end
view.ROIs(view.selectedROI).coords = coords;
view.ROIs(view.selectedROI).modified = datestr(now);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function m = p03_m ( )

%*****************************************************************************80
%
%% P03_M returns the number of equations M for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer M, the number of equations.
%
  m = 5;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [bResult] = gui_IsStringPositiveFloat(sInput, sMessage)

bResult = true;
Tmp = str2double(sInput);
if (isnan(Tmp) | (Tmp <= 0))
  errordlg([sMessage ' must be a positive float value']);
  bResult = false;
end

author: CelsoReyes


text:
function amplitudes = GetSpikeAmplitudes(units,varargin)

%GetSpikeAmplitudes - Get spike amplitudes.
%
%  USAGE
%
%    amplitudes = GetSpikeAmplitudes(units,<options>)
%
%    units          list of [electrode group,cluster] pairs
%                   special conventions:
%                     cluster = -1   all clusters except artefacts and MUA
%                     cluster = -2   all clusters except artefacts
%                     cluster = -3   all clusters
%                   (artefacts are assumed to be in cluster 0, and MUA in 1)
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'restrict'    list of time intervals to read from file
%    =========================================================================
%
%  SEE
%
%    See also GetSpikeTimes, GetSpikeWaveforms.

% Copyright (C) 2004-2013 by Michaël Zugaro
%
% This program 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.

global DATA;
if isempty(DATA),
	error('No session defined (did you forget to call SetCurrentSession? Type ''help Data'' for details).');
end

% Default values
intervals = [];

if nargin < 1,
	error('Incorrect number of parameters (type ''help GetSpikeAmplitudes'' for details).');
end

if ~isimatrix(units,'>=-3') || size(units,2) ~= 2,
	error('Incorrect unit list (type ''help GetSpikeAmplitudes'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+1) ' is not a property (type ''help GetSpikeAmplitudes'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'restrict',
		intervals = varargin{i+1};
		if ~isdmatrix(intervals) || size(intervals,2) ~= 2,
			error('Incorrect value for property ''restrict'' (type ''help GetSpikeAmplitudes'' for details).');
		end
		otherwise,
		error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help GetSpikeAmplitudes'' for details).']);
	end
end

groups = unique(units(:,1));

% Check for incompatible list items
for i = 1:length(groups),
	group = groups(i);
	clusters = units(units(:,1)==group,2);
	if	(ismember(-1,clusters) && ~all(clusters==-1)) || (ismember(-2,clusters) && ~all(clusters==-2)) || (ismember(-3,clusters) && ~all(clusters==-3)),
		error('Incompatible list items (type ''help GetSpikeAmplitudes'' for details).');
	end
end

amplitudes = [];
for i = 1:length(groups),

	% Clusters for this group
	group = groups(i);
	clusters = units(units(:,1)==group,2);

	% Load all waveforms for this group
	filename = [DATA.session.path '/' DATA.session.basename '.spk.' int2str(group)];
	nChannels = length(DATA.spikeGroups.groups{group});
	nSamplesPerWaveform = DATA.spikeGroups.nSamples(group);
	peak = DATA.spikeGroups.peakSamples(group);
	a = LoadSpikeAmplitudes(filename,nChannels,nSamplesPerWaveform,peak,DATA.rates.wideband);

	% Select appropriate clusters
	if ismember(-1,clusters),
		keep = a(:,3) ~= 0 & a(:,3) ~= 1;
	elseif ismember(-2,clusters),
		keep = a(:,3) ~= 0;
	elseif ismember(-3,clusters),
		keep = logical(ones(size(a,1),1));
	else
		keep = ismember(a(:,3),clusters);
	end
	a = a(keep,:);

	% Select timeframes
	if ~isempty(intervals),
		a = Restrict(a,intervals);
	end

	% Add to growing list
	amplitudes = [amplitudes;a];

end

% Sort list by time
amplitudes = sortrows(amplitudes);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%
% MATLAB code that performs Homomorphic filtering, Using Butterworth
% High Pass Filter for performing filtering.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
close all
clear all
d=10;
order=2;
im=double(imread('tun.jpg'));
subplot(121)
imshow(im./255);
[r c]=size(im);
homofil(im,d,r,c,order);

author: Sable

repo: mcbench-benchmarks


text:
function r = ne(p,q)
% this function compares the inequality condition of x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if (isa(p, 'FaceVariable')&&isa(q, 'FaceVariable'))
    r=p;
    r.xvalue = p.xvalue~=q.xvalue;
    r.yvalue = p.yvalue~=q.yvalue;
    r.zvalue = p.zvalue~=q.zvalue;
elseif isa(p, 'FaceVariable')
    r=p;
    r.xvalue = p.xvalue~=q;
    r.yvalue = p.yvalue~=q;
    r.zvalue = p.zvalue~=q;
else
    r=q;
    r.xvalue = p~=q.xvalue;
    r.yvalue = p~=q.yvalue;
    r.zvalue = p~=q.zvalue;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function b = dgesl ( a, lda, n, ipvt, b, job )

%*****************************************************************************80
%
%% DGESL solves a real general linear system A * X = B.
%
%  Discussion:
%
%    DGESL can solve either of the systems A * X = B or A' * X = B.
%
%    The system matrix must have been factored by DGECO or DGEFA.
%
%    A division by zero will occur if the input factor contains a
%    zero on the diagonal.  Technically this indicates singularity
%    but it is often caused by improper arguments or improper
%    setting of LDA.  It will not occur if the subroutines are
%    called correctly and if DGECO has set 0.0 < RCOND
%    or DGEFA has set INFO == 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real A(LDA,N), the output from DGECO or DGEFA.
%
%    Input, integer LDA, the leading dimension of A.
%
%    Input, integer N, the order of the matrix A.
%
%    Input, integer IPVT(N), the pivot vector from DGECO or DGEFA.
%
%    Input, real B(N), the right hand side vector.
%
%    Input, integer JOB.
%    0, solve A * X = B;
%    nonzero, solve A' * X = B.
%
%    Output, real B(N), the solution vector.
%

%
%  Solve A * X = B.
%
  if ( job == 0 )

    for k = 1 : n-1

      l = ipvt(k);
      t = b(l);

      if ( l ~= k )
        b(l) = b(k);
        b(k) = t;
      end

      b(k+1:n) = daxpy ( n-k, t, a(k+1:n,k)', 1, b(k+1:n), 1 );

    end

    for k = n : -1 : 1
      b(k) = b(k) / a(k,k);
      t = -b(k);
      b(1:k-1) = daxpy ( k-1, t, a(1:k-1,k)', 1, b(1:k-1), 1 );
    end

  else
%
%  Solve A' * X = B.
%
    for k = 1 : n
      t = b(1:k-1) * a(1:k-1,k);
      b(k) = ( b(k) - t ) / a(k,k);
    end

    for k = n-1 : -1 : 1

      b(k) = b(k) + b(k+1:n) * a(k+1:n,k);
      l = ipvt(k);

      if ( l ~= k )
        [ b(l), b(k) ] = r8_swap ( b(l), b(k) );
      end

    end

  end

  return
end


text:
classdef nninterp < nntest
  methods (Test)

    function basicShrink(test)
      zoom = 1 ;
      shrink = 3 ;
      batchSize = 10 ;
      x = test.randn([5 5 3 batchSize]) ;
      y = vl_nninterp(x, shrink, zoom) ;

      % check derivatives with numerical approximation
      dzdy = test.randn(size(y)) ;
      dzdx = vl_nninterp(x, shrink, zoom, dzdy) ;
      test.der(@(x) vl_nninterp(x, shrink, zoom), x, dzdy, dzdx, 1e-3*test.range) ;
    end

    function basicShrinkZoom(test)
      zoom = 4 ;
      shrink = 3 ;
      batchSize = 10 ;
      padBeg = 2 ;
      padEnd = 1 ;
      pad = {'padBeg', padBeg, 'padEnd', padEnd} ;
      x = test.randn([5 5 3 batchSize]) ;
      y = vl_nninterp(x, shrink, zoom, pad{:}) ;

      % check derivatives with numerical approximation
      dzdy = test.randn(size(y)) ;
      dzdx = vl_nninterp(x, shrink, zoom, dzdy, pad{:}) ;
      test.der(@(x) vl_nninterp(x, shrink, zoom, pad{:}), ...
                               x, dzdy, dzdx, 1e-3*test.range) ;
    end

  end
end

repo: MFT


text:
function insidepoly_sglengine

error('Please compile MEX files by runing insidepoly_install');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Copyright (C) Daphne Koller, Stanford University, 2012

function F = NormalizeFactorValues( F )

  for i=1:length(F)
    ThisFactor = F(i);
    ThisFactor.val = ThisFactor.val / sum(ThisFactor.val);
    F(i) = ThisFactor;
  end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
clear

load('menpo_68_pts_flip');

pdmLoc = ['../../models/pdm/pdm_68_aligned_wild.mat'];
load(pdmLoc);
% Before plugging into NRSFM, align the points in 2D
num_pts = size(all_pts,1)/2;

num_lmks = numel(M) / 3;
m = reshape(M, num_lmks, 3)';
width_model = max(m(1,:)) - min(m(1,:));
height_model = max(m(2,:)) - min(m(2,:));

for i=1:size(all_pts,1)/2

    shape2D = cat(2, all_pts(i,:)', all_pts(i+num_pts,:)');

    M_n = M;

    if(sum(shape2D(:)==-1) > 0)

        hidden = true;
        % which indices to remove
        inds_to_rem = shape2D(:,1) == -1 | shape2D(:,2) == -1;

        shape2D = shape2D(~inds_to_rem,:);

        inds_to_rem = repmat(inds_to_rem, 3, 1);

        M_n = M(~inds_to_rem);

    end

    % To deal with really extreme cases of roll
    M2D = cat(2, M_n(1:end/3), M_n(end/3+1:2*end/3));
    [ A, t, error, alignedShape, s ] = AlignShapesWithScale(M2D, shape2D);
    R = A/s;

    % Transform the shape
    shape2D(:,1) = shape2D(:,1) - t(1);
    shape2D(:,2) = shape2D(:,2) - t(2);

    shape2D = (R' * shape2D')/s;
    shape2D = shape2D';

    all_pts(i,all_pts(i,:)~=-1) = shape2D(:,1);
    all_pts(i+num_pts,all_pts(i+num_pts,:)~=-1) = shape2D(:,2);
end

%%
% to_rem = randperm(round(0.1*num_pts)); % Remove 10% to break some symmetry
% all_pts([to_rem, to_rem+num_pts],:) = [];
num_pts = size(all_pts,1)/2;

left_ids = all_pts(1:num_pts,10) == -1;
right_ids = all_pts(1:num_pts,8) == -1;
frontal = true(num_pts,1);
frontal(left_ids | right_ids) = false;

%%
xs = all_pts(1:num_pts,:);
ys = all_pts(num_pts+1:end,:);

% Randperm the data, as a test

% xs_f = xs(frontal,:);
% ys_f = ys(frontal,:);
% scatter(xs_f(:), -ys_f(:));

%% Perform NRSFM by Torresani
addpath('../nrsfm-em');
% (T is the number of frames, J is the number of points)

J = size(all_pts,2);
T = size(all_pts,1)/2;

use_lds = 0; % not modeling a linear dynamic system here
max_em_iter = 200;
tol = 0.001;
K = 30; % number of deformation shapes

MD = all_pts(1:end/2,:)==-1;

[P3, S_hat, V, RO, Tr, Z] = em_sfm(all_pts, MD, K, use_lds, tol, max_em_iter);

save('Torr_menpo', 'P3', 'S_hat', 'V', 'RO', 'Tr', 'Z');

%%
% xs = P3(1:num_pts,:);
% ys = P3(num_pts+1:2*num_pts,:);
% zs = P3(2*num_pts+1:end,:);
% %
% xs_f = xs(left_ids,:);
% ys_f = ys(left_ids,:);
% zs_f = zs(left_ids,:);
% scatter3(xs_f(:), ys_f(:), zs_f(:));


text:
function [x,varargout] = blasst(x,lineFrequencies,frequencyRanges,samplingRate,varargin)
% blasst(): EEGLAB helper function for OCW line noise removal.
% Takes as input an array of signals x, along with relevant parameters, and
% performs BLASST filtering at specified frequencies. For each specified
% frequency, blasst() iteratively calls blasst_internal() and then uses
% blasst_test() to test for convergence based on the distributions of
% convolution coefficients in the target and surrounding frequency bands.
%
% INPUT:
% x                 an [n,N] array of n signals of length N.
% lineFrequencies   an array of target frequencies (not normalized).
% frequenyRanges    an array of target frequency ranges. Must be the same
%                   size as lineFrequencies.
% samplingRate      the sampling rate of the signal.
% varargin          optional 'key',value pairs:
%   'key'
%           [default] purpose
%   'Scale'
%           [2^(log2(samplingRate)+2)] Manually set the scale that
%           indicates the spread of Gabor atoms.
%   'ContinuousEpochs'
%           [0] If x is epoched, but epochs are temporally adjacent,
%           setting to 1 will flatten x for processing. Otherwise,
%           blasst is run on individual epochs.
%   'Verbose'
%           [1] When on, progress is printed on command line.
%   'Resolution'
%           [2] May be an integer value >= 1, sets 'resolution' in blasst.
%           Specificies density of Gabor atoms. May also be an array of
%           integer values of size(lineFrequencies).
%   'MaxIterations'
%           [50] Maximum number of external iterations of blasst run on
%           each frequency. May be either a scalar integer or array of
%           integers of size(lineFrequencies).
%   'ManualOffset'
%           [log2(scaleBases)+1] A scalar value that offsets the
%           arrangement of Gabor atoms at each iteration of blasst.
%   'Channels'
%           [1:size(x,1)] An array of integers indexing channels to
%           be computed. Allows manually selection of channels for
%           processing.
%
% OUTPUT:
% x                 the processed, or ``cleaned'' signal.
% varargout{1}      the aggregate of the target signal feature removed, an
%                   array of size(x).
%
% DEPENDENCIES:
% blasst_internal()         primary line noise removal algorithm.
% blasst_test()        convergence test for iterative blasst algorithm.
%
% EXAMPLE:
%  Suppose we wish to remove line noise frequency at 60 and 120 Hz, and the
%  noise is mostly stationary at 120 Hz but non-stationary varying by about
%  2 Hz, around 60 Hz. Then we might call the method as:
% >> x = blasst(x,[60,120],[2,.25],<sampling rate>);
%
%  If we want to use more densely packed Gabor atoms, we could call:
% >> x = blasst(x,[60,120],[2,.25],<sampling rate>,'Resolution',4);
%
% AUTHOR: Kenneth Ball, 2015.
%
% IF YOU FIND BLASST USEFUL IN YOUR WORK, PLEASE CITE:
%
% Ball, K. R., Hairston, W. D., Franaszczuk, P. J., Robbins, K. A.,
% BLASST: Band Limited Atomic Sampling with Spectral Tuning with
% Applications to Utility Line Noise Filtering, [Under Review].
%
% Copyright 2015 Kenneth Ball
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

% Set defaults:
flattenData = 0;
verbose = 1;
resolution = 2;
maxIterations = 50; % Default is high so that the method generally will generally converge.
channels = 1:size(x,1);

% Adjust for optional inputs
if (nargin> 1)
    if (nargin> 4 && rem(nargin,2) == 1)
        if length(varargin) == 1
            varargin = varargin{1};
        else
            fprintf('blasst(): Optional key and value pairs do not match.')
            return
        end
    end

    for ii = 1:2:length(varargin)
        key = varargin{ii};
        val = varargin{ii+1};
        switch key
            case 'SamplingRate'
                samplingRate = val;
            case 'Scale'
                scale = val;
            case 'ContinuousEpochs'
                flattenData = val;
            case 'Verbose'
                verbose = val;
            case 'Resolution'
                resolution  = val;
            case 'MaxIterations'
                maxIterations = val;
            case 'ManualOffset'
                manualOffset = val;
            case 'Channels'
                channels = val;
        end
    end
end

if length(frequencyRanges) == 1
    frequencyRanges = ones(size(lineFrequencies))*frequencyRanges;
elseif length(frequencyRanges) ~= length(lineFrequencies)
    error('Number of specified noise frequencies does not match number of specified range values.');
end

if length(resolution) == 1
    resolution = ones(size(lineFrequencies))*resolution;
elseif length(resolution) ~= length(lineFrequencies)
    error('Number of specified noise resolutions does not match number of specified range values.');
end

if length(maxIterations) == 1
    maxIterations = ones(size(lineFrequencies))*maxIterations;
elseif length(maxIterations) ~= length(lineFrequencies)
    error('Number of specified max iterations does not match number of specified range values.');
end

if ~exist('samplingRate','var')
    error('No sampling rate specified.');
elseif isempty(samplingRate) || samplingRate == 0
    error('Null or zero sampling rate specified.');
end

if ~exist('scaleBases','var')
    scale = 2^(log2(samplingRate)+2);
end

if ~exist('manualOffset','var')
    manualOffset = log2(scale)+1;
end

if flattenData
    dataSizeTemp = size(x);
    x = reshape(x,size(x,1),size(x,2)*size(x,3));
end

% Initialize holders for fitted noise and cleaned signals.
y = zeros(size(x));
yTemp = zeros(size(x));
xTemp = zeros(size(x));

countTrack = 1;

for ii = 1:length(lineFrequencies)

    if verbose
        fprintf(1,'Frequency: %d Hz\n',lineFrequencies(ii));
    end
    for jj = channels
        if verbose
            fprintf(1,'Computing Channel:                ');
        end
        %         BDist = Inf;
        [BDist,~,~] = blasst_test(reshape(x(jj,:,:),1,size(x,2)*size(x,3)),lineFrequencies(ii),frequencyRanges(ii),samplingRate,scale);
        for mm = 1:maxIterations(ii)
            if verbose
                fprintf(1,'\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b% 3i   Pass: % 3i',jj,mm);
            end
            % Run blasst_internal for each epoch of data in the EEG struct.
            for kk = 1:size(x,3)
                [yTemp(jj,:,kk),xTemp(jj,:,kk)] = blasst_internal(x(jj,:,kk),lineFrequencies(ii),frequencyRanges(ii),samplingRate,scale,resolution(ii),(mm-1)*manualOffset);
                %                 y(jj,:,kk) = y(jj,:,kk) + tempY;
                %                 EEG.data(jj,:,kk) = tempX;
            end

            % Compute Bhatt. distance for testing convergence. Overwrite
            % BDist1, then compare to BDist (from the last pass). If BDist1
            % exceeds BDist, we presume we have passed the minimum distance
            % between the distributions, and we should halt the algorithm
            % for this channel and frequency, retaining only previous
            % passes.

            [BDist1,maxFlag,~] = blasst_test(reshape(xTemp(jj,:,:),1,size(xTemp,2)*size(xTemp,3)),lineFrequencies(ii),frequencyRanges(ii),samplingRate,scale);
            countTrack = countTrack+1;
            if BDist1 >= BDist && ~maxFlag % && mm > 1
                if verbose
                    fprintf(1,'\nBreak at pass % 2i\n',mm-1);
                end
                break
            else
                BDist = BDist1;
                y(jj,:,:) = y(jj,:,:) + yTemp(jj,:,:);
                x(jj,:,:) = xTemp(jj,:,:);
            end
        end
        if verbose
            fprintf(1,'\n');
        end
    end
    if verbose
        fprintf(1,'\n');
    end

end

varargout{1} = y; % This is the aggregate of all line noise that was removed. That is, y+EEG.data is the original dataset.

if flattenData
    x = reshape(x,dataSizeTemp(1),dataSizeTemp(2),dataSizeTemp(3));
end

end

author: VisLab

repo: EEG-Clean-Tools


text:
function image_db = load_image_dataset(varargin)
%
% This file is part of the code that implements the following paper:
% Title      : "LocNet: Improving Localization Accuracy for Object Detection"
% Authors    : Spyros Gidaris, Nikos Komodakis
% Institution: Universite Paris Est, Ecole des Ponts ParisTech
% ArXiv link : http://arxiv.org/abs/1511.07763
% code       : https://github.com/gidariss/LocNet
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2016 Spyros Gidaris
%
% Title     : "LocNet: Improving Localization Accuracy for Object Detection"
% ArXiv link: http://arxiv.org/abs/1511.07763
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

% ************************** OPTIONS *************************************
ip = inputParser;
ip.addParamValue('dataset', 'pascal');
ip.addParamValue('image_set', {})
ip.addParamValue('voc_year',  {})
ip.addParamValue('proposals_method', {});
ip.addParamValue('use_flips', false, @islogical);

ip.parse(varargin{:});
opts = ip.Results;

voc_year         = opts.voc_year;
image_set        = opts.image_set;
proposals_method = opts.proposals_method;
use_flips        = opts.use_flips;

voc_path  = [pwd, '/datasets/VOC%s/'];
coco_path = [pwd, '/datasets/MSCOCO/'];

if ~iscell(proposals_method), proposals_method = {proposals_method}; end
if ~iscell(voc_year), voc_year = {voc_year}; end
if ~iscell(image_set), image_set = {image_set}; end

assert(ischar(voc_year{1}));

num_sets = length(image_set);
assert(length(image_set) == length(voc_year));
image_db.all_regions = {};
for i = 1:num_sets
    % add the paths to the images and the grount truth bounding boxes in
    % the image_db structure
    switch opts.dataset
        case 'pascal'
            image_db_this = load_pascal_dataset(...
                voc_year{i}, image_set{i}, voc_path);
        case 'mscoco'
            image_db_this = load_mscoco_dataset(...
                voc_year{i}, image_set{i}, coco_path );
    end
    % add the bounding box proposals in the image_db structure
    if ~isempty(proposals_method)
        assert(ischar(proposals_method{1}));
        image_db.all_regions = [image_db.all_regions(:); ...
            load_box_proposals(image_db_this, proposals_method)];
    end
    image_db_all(i) = image_db_this;
end
image_db_all = image_db_all(:);
image_db.all_bbox_gt = vertcat(image_db_all.all_bbox_gt);
image_db.image_paths = vertcat(image_db_all.image_paths);
image_db.image_sizes = single(vertcat(image_db_all.image_sizes));

proposals_suffix = '';
if ~isempty(proposals_method)
    for p = 1:length(proposals_method)
        proposals_suffix = [proposals_suffix, proposals_method{p}, '_'];
    end
    proposals_suffix = proposals_suffix(1:end-1);
end

switch opts.dataset
    case 'pascal', image_set_name_all = 'voc_';
    case 'mscoco', image_set_name_all = 'coco_';
end

for i = 1:num_sets
    image_set_name_all = [image_set_name_all, voc_year{i}, '_' image_set{i},'_'];
end
image_set_name_all = image_set_name_all(1:end-1);

image_db.image_set_name = image_set_name_all;
image_db.proposals_suffix = proposals_suffix;

% add flipped versions of the images
if use_flips, image_db = add_flip_data(image_db); end

switch opts.dataset
    case 'pascal'
        voc_path_year = sprintf(voc_path, '2007');
        VOCopts  = initVOCOpts(voc_path_year,'2007');
        classes = VOCopts.classes;
    case 'mscoco'
        save_file_categories = fullfile(coco_path,'annotations','categoties.mat');
        ld = load(save_file_categories,'category_names','cat_ids','catIdsToIndex');
        classes = ld.category_names;
end

image_db.classes = classes;

end

function image_db = load_pascal_dataset(voc_year, image_set, voc_path)
image_db = struct;
voc_path_dataset = sprintf(voc_path, voc_year);
assert(exist(voc_path_dataset,'dir')>0);
voc_path_devkit = fullfile(voc_path_dataset,'VOCdevkit');
assert(exist(voc_path_devkit,'dir')>0);

[image_db.image_paths, image_db.image_set_name] = get_image_paths_from_voc( voc_path_devkit, image_set, voc_year );
image_db.all_bbox_gt = get_grount_truth_bboxes_from_voc(voc_path_devkit, image_set, voc_year, true, voc_path_dataset);
image_db.image_sizes = get_img_size(image_db.image_paths);
end

function image_db = load_mscoco_dataset(coco_year, image_set, data_dir )
image_db = struct;

suffix = sprintf('%s%s',image_set,coco_year);
image_db.image_set_name = sprintf('coco_%s_%s',coco_year,image_set);
save_file_bbox_gt = fullfile(data_dir,'annotations',sprintf('all_bbox_gt_%s.mat',suffix));

save_file_image_names = fullfile(data_dir,'images',sprintf('all_image_names_%s.mat',suffix));
ld = load(save_file_image_names, 'image_names', 'image_sizes');
image_names          = ld.image_names;
image_db.image_paths = strcat(sprintf('%s/images/%s/', data_dir, suffix), image_names);
image_db.image_sizes = ld.image_sizes;

if (exist(save_file_bbox_gt,'file') > 0)
    ld = load(save_file_bbox_gt, 'all_bbox_gt');
    image_db.all_bbox_gt = ld.all_bbox_gt;
else
    image_db.all_bbox_gt = cell(length(image_db.image_paths),1);
    for i = 1:length(image_db.all_bbox_gt)
        image_db.all_bbox_gt{i} = zeros(0,6,'single');
    end
end

end

function image_sizes = get_img_size(image_paths)
num_imgs    = numel(image_paths);
image_sizes = zeros(num_imgs,2);

for img_idx = 1:num_imgs
    im_info                = imfinfo(image_paths{img_idx});
    image_sizes(img_idx,1) = im_info.Height;
    image_sizes(img_idx,2) = im_info.Width;
end
end

function image_db = add_flip_data(image_db)
num_imgs          = numel(image_db.image_paths);
image_paths_new   = cell(num_imgs,1);

there_are_regions = isfield(image_db, 'all_regions') && ~isempty(image_db.all_regions);
there_are_bbox_gt  = isfield(image_db, 'all_bbox_gt')  && ~isempty(image_db.all_bbox_gt);
there_are_feature_paths = isfield(image_db, 'feature_paths')  && ~isempty(image_db.feature_paths);

if there_are_regions,   all_regions_new = image_db.all_regions; end
if there_are_bbox_gt,   all_bbox_gt_new = image_db.all_bbox_gt; end
if there_are_feature_paths, feature_paths_new = cell(num_imgs,1); end

for img_idx = 1:num_imgs
    [img_dir, img_name, img_ext]  = fileparts(image_db.image_paths{img_idx});
    image_paths_new{img_idx} = [img_dir, filesep, img_name, '_flip', img_ext];

    if there_are_feature_paths
        [feat_dir, feat_name, feat_ext]  = fileparts(image_db.feature_paths{img_idx});
        feature_paths_new{img_idx} = [feat_dir, filesep, feat_name, '_flip', feat_ext];
    end

    img_weight = image_db.image_sizes(img_idx,2);
    if there_are_regions
        all_regions_new{img_idx}(:,[1,3]) = img_weight + 1 - all_regions_new{img_idx}(:,[3,1]);
    end
    if there_are_bbox_gt
        all_bbox_gt_new{img_idx}(:,[1,3]) = img_weight + 1 - all_bbox_gt_new{img_idx}(:,[3,1]);
    end
end

image_db.image_paths = cat(1, image_db.image_paths, image_paths_new);
image_db.image_sizes = cat(1, image_db.image_sizes,  image_db.image_sizes);

if there_are_bbox_gt, image_db.all_bbox_gt = cat(1, image_db.all_bbox_gt, all_bbox_gt_new); end
if there_are_regions, image_db.all_regions = cat(1, image_db.all_regions, all_regions_new); end
if there_are_feature_paths, image_db.feature_paths = cat(1,image_db.feature_paths, feature_paths_new); end
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function mesh = prepare_mesh_headshape(cfg)

% PREPARE_MESH_HEADSHAPE
%
% Configuration options:
%   cfg.headshape   = a filename containing headshape, a Nx3 matrix with surface
%                     points, or a structure with a single or multiple boundaries
%   cfg.smooth      = a scalar indicating the number of non-shrinking
%                     smoothing iterations (default = no smoothing)
%
% See also PREPARE_MESH_MANUAL, PREPARE_MESH_SEGMENTATION

% Copyrights (C) 2009, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% get the specific options
cfg.headshape    = ft_getopt(cfg, 'headshape');
cfg.smooth       = ft_getopt(cfg, 'smooth');   % no default

if isa(cfg, 'config')
  % convert the config-object back into a normal structure
  cfg = struct(cfg);
end

if isa(cfg.headshape, 'config')
  % convert the nested config-object back into a normal structure
  cfg.headshape = struct(cfg.headshape);
end

% get the surface describing the head shape
if isstruct(cfg.headshape) && numel(cfg.headshape)>1
  % this applies for multilayer BEM models and concentric sphere models
  headshape = [];
  for i=1:numel(cfg.headshape)
    [headshape(i).pos, headshape(i).tri] = headsurface([], [], 'headshape', cfg.headshape(i));
  end
else
  [headshape.pos, headshape.tri] = headsurface([], [], 'headshape', cfg.headshape);
end

if numel(headshape)>1 && numel(cfg.numvertices)==1
  % use the same number of vertices for each of the head shape surfaces
  cfg.numvertices = repmat(cfg.numvertices, size(headshape));
end

if ~isempty(cfg.numvertices) && ~isequal(cfg.numvertices, arrayfun(@(x) size(x.pos, 1), headshape))
  for i=1:numel(headshape)
    tri1 = headshape(i).tri;
    pos1 = headshape(i).pos;
    % The number of vertices is multiplied by 3 in order to have more
    % points on the original mesh than on the sphere mesh (see below).
    % The rationale for this is that every projection point on the sphere
    % has three corresponding points on the mesh
    if (cfg.numvertices(i)>size(pos1,1))
      [tri1, pos1] = refinepatch(headshape(i).tri, headshape(i).pos, 3*cfg.numvertices(i));
    else
      [tri1, pos1] = reducepatch(headshape(i).tri, headshape(i).pos, 3*cfg.numvertices(i));
    end

    % remove double vertices
    [pos1, tri1] = remove_double_vertices(pos1, tri1);

    % replace the probably unevenly distributed triangulation with a regular one
    % and retriangulate it to the desired accuracy
    [pos2, tri2] = mysphere(cfg.numvertices(i)); % this is a regular triangulation
    [pos1, tri1] = retriangulate(pos1, tri1, pos2, tri2, 2);
    [pos1, tri1] = fairsurface(pos1, tri1, 1); % this helps redistribute the superimposed points

    % remove double vertices
    [headshape(i).pos, headshape(i).tri] = remove_double_vertices(pos1, tri1);
    fprintf('returning %d vertices, %d triangles\n', size(headshape(i).pos,1), size(headshape(i).tri,1));
  end
end

% smooth the mesh
if ~isempty(cfg.smooth)
  for i=1:numel(headshape)
    [headshape(i).pos,headshape(i).tri] = fairsurface(headshape(i).pos, headshape(i).tri, cfg.smooth);
  end
end

% the output should only describe one or multiple boundaries and should not
% include any other fields
mesh = keepfields(headshape, {'pos', 'tri'});

function [tri1, pos1] = refinepatch(tri, pos, numvertices)
fprintf('the original mesh has %d vertices, the requested number of vertices is %d\n',size(pos,1),numvertices/3);
fprintf('trying to refine the mesh...\n');
[pos1, tri1] = refine(pos, tri, 'updown', numvertices);

function [pos, tri] = mysphere(N)
% This is a copy of MSPHERE without the confusing output message
% Returns a triangulated sphere with approximately M vertices
% that are nicely distributed over the sphere. The vertices are aligned
% along equally spaced horizontal contours according to an algorithm of
% Dave Russel.
%
% Use as
%  [pos, tri] = msphere(M)
%
% See also SPHERE, NSPHERE, ICOSAHEDRON, REFINE
% Copyright (C) 1994, Dave Rusin

storeM    = [];
storelen  = [];
increaseM = 0;
while (1)

  % put a single vertex at the top
  phi = [0];
  th  = [0];

  M = round((pi/4)*sqrt(N)) + increaseM;
  for k=1:M
    newphi = (k/M)*pi;
    Q = round(2*M*sin(newphi));
    for j=1:Q
      phi(end+1) = newphi;
      th(end+1)  = (j/Q)*2*pi;
      % in case of even number of contours
      if mod(M,2) && k>(M/2)
        th(end) = th(end) + pi/Q;
      end
    end
  end

  % put a single vertex at the bottom
  phi(end+1) = [pi];
  th(end+1)  = [0];

  % store this vertex packing
  storeM(end+1).th  = th;
  storeM(end  ).phi = phi;
  storelen(end+1) = length(phi);
  if storelen(end)>N
    break;
  else
    increaseM = increaseM+1;
    % fprintf('increasing M by %d\n', increaseM);
  end
end

% take the vertex packing that most closely matches the requirement
[m, i] = min(abs(storelen-N));
th  = storeM(i).th;
phi = storeM(i).phi;

% convert from spherical to cartehsian coordinates
[x, y, z] = sph2cart(th, pi/2-phi, 1);
pos = [x' y' z'];
tri = convhulln(pos);

function [pos1, tri1] = fairsurface(pos, tri, N)

% FAIRSURFACE modify the mesh in order to reduce overlong edges, and
% smooth out "rough" areas. This is a non-shrinking smoothing algorithm.
% The procedure uses an elastic model : At each vertex, the neighbouring
% triangles and vertices connected directly are used. Each edge is
% considered elastic and can be lengthened or shortened, depending
% on their length. Displacement are done in 3D, so that holes and
% bumps are attenuated.
%
% Use as
%   [pos, tri] = fairsurface(pos, tri, N);
% where N is the number of smoothing iterations.
%
% This implements:
%   G.Taubin, A signal processing approach to fair surface design, 1995

% This function corresponds to spm_eeg_inv_ElastM
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
%                    Christophe Phillips & Jeremie Mattout
% spm_eeg_inv_ElastM.m 1437 2008-04-17 10:34:39Z christophe
%
% $Id$

ts = [];
ts.XYZmm = pos';
ts.tri   = tri';
ts.nr(1) = size(pos,1);
ts.nr(2) = size(tri,1);

% Connection vertex-to-vertex
%--------------------------------------------------------------------------
M_con = sparse([ts.tri(1,:)';ts.tri(1,:)';ts.tri(2,:)';ts.tri(3,:)';ts.tri(2,:)';ts.tri(3,:)'], ...
  [ts.tri(2,:)';ts.tri(3,:)';ts.tri(1,:)';ts.tri(1,:)';ts.tri(3,:)';ts.tri(2,:)'], ...
  ones(ts.nr(2)*6,1),ts.nr(1),ts.nr(1));

kpb   = .1;                       % Cutt-off frequency (default: .1)
lam   = .5; mu = lam/(lam*kpb-1); % Parameters for elasticity. (default: .5)
XYZmm = ts.XYZmm;

% smoothing iterations
%--------------------------------------------------------------------------
for j=1:N

  XYZmm_o = zeros(3,ts.nr(1)) ;
  XYZmm_o2 = zeros(3,ts.nr(1)) ;

  for i=1:ts.nr(1)
    ln = find(M_con(:,i));
    d_i = sqrt(sum((XYZmm(:,ln)-XYZmm(:,i)*ones(1,length(ln))).^2));
    if sum(d_i)==0
      w_i = zeros(size(d_i));
    else
      w_i = d_i/sum(d_i);
    end
    XYZmm_o(:,i) = XYZmm(:,i) + ...
      lam * sum((XYZmm(:,ln)-XYZmm(:,i)*ones(1,length(ln))).*(ones(3,1)*w_i),2);
  end

  for i=1:ts.nr(1)
    ln = find(M_con(:,i));
    d_i = sqrt(sum((XYZmm(:,ln)-XYZmm(:,i)*ones(1,length(ln))).^2));
    if sum(d_i)==0
      w_i = zeros(size(d_i));
    else
      w_i = d_i/sum(d_i);
    end
    XYZmm_o2(:,i) = XYZmm_o(:,i) + ...
      mu * sum((XYZmm_o(:,ln)-XYZmm_o(:,i)*ones(1,length(ln))).*(ones(3,1)*w_i),2);
  end

  XYZmm = XYZmm_o2;

end

% collect output results
%--------------------------------------------------------------------------

pos1 = XYZmm';
tri1 = tri;

if 0
  % this is some test/demo code
  mesh = [];
  [mesh.pos, mesh.tri] = mesh_sphere(162);

  scale = 1+0.3*randn(size(pos,1),1);
  mesh.pos = mesh.pos .* [scale scale scale];

  figure
  ft_plot_mesh(mesh)

  [mesh.pos, mesh.tri] = fairsurface(mesh.pos, mesh.tri, 10);

  figure
  ft_plot_mesh(mesh)
end


text:
%
% Performs contour (circular) integral with
% center x_0,y_0 and radius r, using discrete riemann approach.
% The angs parameter specifies region of the circle
% considering clockwise 0-2pi notation.
% Result not normalized.
% Author: J. Rouces
%
function sum = ContourIntegralCircular(imagen,y_0,x_0,r,angs)

% EXHAUSTIVE EXTENSIVE ALGORITHM
% rc = r^2;
% sum = 0;
% for x = max(1,x_0-r):min(size(imagen,2),x_0+r)
%     for y = max(1,y_0-r):min(size(imagen,1),y_0+r)
%         if abs((x-x_0)^2+(y-y_0)^2-rc)<2
%             sum = sum + imagen(y,x);
%         end
%     end
% end
% sum = sum/r;

% LIGHT ALGORITHM
sum = 0;
for ang = angs
    y = round(y_0-cos(ang)*r);
    x = round(x_0+sin(ang)*r);
    if y<1
        y = 1;
    elseif y>size(imagen,1)
        y = size(imagen,1);
    end
    if x<1
        x = 1;
    elseif x>size(imagen,2)
        x = size(imagen,2);
    end
    sum = sum + imagen(y,x);
end

author: Qingbao

repo: iris

save_path: github-repos/MATLAB/Qingbao-iris


text:
function x = grid2 ( i1, i2, ndim, nstep, x1, x2 )

%*****************************************************************************80
%
%% GRID2 computes grid points between X1 and X2 in N dimensions.
%
%  Discussion:
%
%    However, X1 need not be the first point computed, nor X2 the last.
%    The user must specify the steps on which X1 and X2 are passed
%    through.  These steps may even be outside the range of 1 through NSTEP.
%
%    We assume that a set of equally spaced points have
%    been drawn on the line through X1 and X2, and that
%    they have been numbered, with X1 labeled I1 and X2
%    labeled I2.  I1 or I2 may be between 1 and NSTEP,
%    in which case X1 or X2 will actually be returned in the
%    X array, but there is no requirement that I1 or I2
%    satisfy this condition.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I1, I2.  I1 specifies the step on which
%    X1 would be computed, and similarly for I2.
%    I1 and I2 must be distinct.
%
%    Input, integer NDIM, the dimension of the points X1 and X2.
%
%    Input, integer NSTEP, this is the number of equally
%    spaced points that are to be generated.
%    NSTEP should be at least 1.
%
%    Input, real X1(NDIM), X2(NDIM), the points that define
%    the line along which the equally spaced points are generated, and
%    which may or may not be included in the set of computed points.
%
%    Output, real X(NDIM,NSTEP), the set of equally spaced points.
%    Each column of X represents one point.
%    If 1 <= I1 <= NSTEP, then X(*,I1) = X1, and similarly for I2.
%
  if ( i1 == i2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GRID2 - Fatal error!\n' );
    fprintf ( 1, '  I1 = I2, leading to zero denominator.\n' );
    fprintf ( 1, '  I1 = %d\n', i1 );
    fprintf ( 1, '  I2 = %d\n', i2 );
    error ( 'GRID2 - Fatal error!' );
  end

  for i = 1 : nstep
    for j = 1 : ndim
      x(j,i) = ( ( i2 - i      ) * x1(j)   ...
               + (      i - i1 ) * x2(j) ) ...
               / ( i2     - i1 );
    end
  end

  return
end


text:
function s=dev(t);
%
%s=dev(t)
%
%This function calculates the deviatoric part s(ij) of a tensor t(ij)
%But, for the purpose of computations we use the vector representation
%of a tensor t(ij), so that:
%t - is a vector representation of a tensor t(ij)
%    with the components in the following order:
%
%t=[txx tyy tzz txy txz tzx tyx tyz txz]
%or
%t=[txx tyy tzz txy txz tyz]
%or
%t=[txx tyy tzz]
%
%Dimension of t can be [m n]
%where: m - represents the number of tensors to calculate
%       n = 9, 6 or 3 (see above)
%
%   Copyright (c) 2001 by Aleksander Karolczuk.
%   $Revision: 1.0 $  $Date: 2001/02/01

error(nargchk(1,1,nargin))

[m n]=size(t);

if n==3
   I=[1 1 1];
elseif n==6
   I=[1 1 1 0 0 0];
elseif n==9
   I=[1 1 1 0 0 0 0 0 0];
else
   error('Improper matrix dimension')
end

s=t-((1/3)*(t*I'))*I;