text:
function [x, c, funVal, ValueL]=tree_mcLogisticR(A, y, z, opts)
%
%%
% Function mcLogisticR:
%      Logistic Loss for Multi-class (task) Learning
%             with the tree structured group Lasso Regularization
%
%% Problem
%
%  min  - sum_{il} weight_{il} log( p_{il} ) + z * sum_i sum_j w_j ||x^i_{G_j}||
%
%  p_{il}= 1 / (1+ exp(-y_i (x_i' * a_i + c_l) ) ) denotes the probability
%  weight_{il} is the weight for the i-th sample in the l-th classifier
%                is a m x k matrix
%  c_l is the intercept for the l-th classfier, and is a 1xk vector
%  x_i denotes the i-th column of x
%  x^j denotes the j-th row of x
%  a_i' denotes the i-th row of A
%
%  y_i (either 1 or -1) is the response
%
%  In this implementation, we assume weight_{il}=1/(mk)
%
%  x^i denotes the i-th row of x
%
%  G_j's are nodes with tree structure
%
%  The tree structured group information is contained in
%  opts.ind, which is a 3 x nodes matrix, where nodes denotes the number of
%  nodes of the tree.
%
%  opts.ind(1,:) contains the starting index
%  opts.ind(2,:) contains the ending index
%  opts.ind(3,:) contains the corresponding weight (w_j)
%
%  Note:
%  1) If each element of x^j is a leaf node of the tree and the weight for
%  this leaf node are the same, we provide an alternative "efficient" input
%  for this kind of node, by creating a "super node" with
%  opts.ind(1,1)=-1; opts.ind(2,1)=-1; and opts.ind(3,1)=the common weight.
%
%  2) If the features are well ordered in that, the features of the left
%  tree is always less than those of the right tree, opts.ind(1,:) and
%  opts.ind(2,:) contain the "real" starting and ending indices. That is to
%  say, x^j( opts.ind(1,j):opts.ind(2,j) ) denotes x^j_{G_j}. In this case,
%  the entries in opts.ind(1:2,:) are within 1 and k (the number of tasks).
%
%
%  If the features are not well ordered, please use the input opts.G for
%  specifying the index so that
%   x^j( opts.G ( opts.ind(1,j):opts.ind(2,j) ) ) denotes x^j_{G_j}.
%  In this case, the entries of opts.G are within 1 and n, and the entries of
%  opts.ind(1:2,:) are within 1 and length(opts.G).
%
% The following example shows how G and ind works:
%
% G={ {1, 2}, {4, 5}, {3, 6}, {7, 8},
%     {1, 2, 3, 6}, {4, 5, 7, 8},
%     {1, 2, 3, 4, 5, 6, 7, 8} }.
%
% ind={ [1, 2, 100]', [3, 4, 100]', [5, 6, 100]', [7, 8, 100]',
%       [9, 12, 100]', [13, 16, 100]', [17, 24, 100]' },
%
% where each node has a weight of 100.
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size m x k)
%  z -        Regularization parameter (z >=0)
%  opts-      optional inputs (default value: opts=[])
%
%% Output parameters:
%  x-         The obtained weight of size n x k
%  c-         The obtained intercept if size 1 x k
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on October 3, 2010.
%
%% Related papers
%
% [1] Jun Liu and Jieping Ye, Moreau-Yosida Regularization for
%     Grouped Tree Structure Learning, NIPS 2010
%
%% Related functions:
%
%  sll_opts
%
%%

%% Verify and initialize the parameters
%%
if (nargin <4)
    error('\n Inputs: A, y, z, and opts (.ind) should be specified!\n');
end

[m,n]=size(A);

k=size(y,2);

if (size(y,1) ~=m)
    error('\n Check the length of y!\n');
end

if (z<0)
    error('\n z should be nonnegative!\n');
end

opts=sll_opts(opts); % run sll_opts to set default values (flags)

% restart the program for better efficiency
%  this is a newly added function
if (~isfield(opts,'rStartNum'))
    opts.rStartNum=opts.maxIter;
else
    if (opts.rStartNum<=0)
        opts.rStartNum=opts.maxIter;
    end
end

%% Detailed initialization
%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Group and others

% Initialize ind
if (~isfield(opts,'ind'))
    error('\n In tree_mcLeastR, the field .ind should be specified');
else
    ind=opts.ind;

    if (size(ind,1)~=3)
        error('\n Check opts.ind');
    end
end

GFlag=0;
% if GFlag=1, we shall apply general_altra
if (isfield(opts,'G'))
    GFlag=1;

    G=opts.G;
    if (max(G) >k || max(G) <1)
        error('\n The input G is incorrect. It should be within %d and %d',1,k);
    end
end

%% Starting point initialization

% % % compute AT y
% % if (opts.nFlag==0)
% %     ATy =A'* y;
% % elseif (opts.nFlag==1)
% %     ATy= A'* y -  mu' * sum(y, 1);    ATy=ATy./repmat(nu, 1, k);
% % else
% %     invNu=y./repmat(nu, 1, k);        ATy=A'*invNu-mu' * sum(invNu, 1);
% % end

p_flag=(y==1);                  % the indices of the postive samples
m1=sum(p_flag,1);               % the total number of the positive samples
m2=m-m1;                        % the total number of the negative samples
% m1 and m2 are 1 x k vectors

% process the regularization parameter
if (opts.rFlag==0)
    lambda=z;
else % z here is the scaling factor lying in [0,1]
%     if (z<0 || z>1)
%         error('\n opts.rFlag=1, and z should be in [0,1]');
%     end

    % we compute ATb for computing lambda_max, when the input z is a ratio

    computedFlag=0;
    if (isfield(opts,'lambdaMax'))
        if (opts.lambdaMax~=-1)
            lambda=z*opts.lambdaMax;
            computedFlag=1;
        end
    end

    if (~computedFlag)

        b=zeros(m, k);                  % b is of size m x k
        b= b + p_flag.*repmat(m2,m,1);
        b= b -(~p_flag).*repmat(m1,m,1);
        b=b/m;

        % compute AT b
        if (opts.nFlag==0)
            ATb =A'* b;
        elseif (opts.nFlag==1)
            ATb= A'* b -  mu' * sum(b, 1);    ATb=ATb./repmat(nu, 1, k);
        else
            invNu=b./repmat(nu, 1, k);        ATb=A'*invNu-mu' * sum(invNu, 1);
        end

        if (GFlag==0)
            lambda_max=findLambdaMax_mt(ATb, n, k, ind, size(ind,2));
        else
            lambda_max=general_findLambdaMax_mt(ATb, n, k, G, ind, size(ind,2));
        end

        lambda_max= lambda_max / (m*k);

        % As .rFlag=1, we set lambda as a ratio of lambda_max
        lambda=z*lambda_max;
    end
end

% The following is for computing lambdaMax
% we use opts.lambdaMax=-1 to show that we need the computation.
%
% One can use this for setting up opts.lambdaMax
if (isfield(opts,'lambdaMax'))

    if (opts.lambdaMax==-1)

        b=zeros(m, k);                  % b is of size m x k
        b= b + p_flag.*repmat(m2,m,1);
        b= b -(~p_flag).*repmat(m1,m,1);
        b=b/m;

        % compute AT b
        if (opts.nFlag==0)
            ATb =A'* b;
        elseif (opts.nFlag==1)
            ATb= A'* b -  mu' * sum(b, 1);    ATb=ATb./repmat(nu, 1, k);
        else
            invNu=b./repmat(nu, 1, k);        ATb=A'*invNu-mu' * sum(invNu, 1);
        end

        % compute lambda_max
        if (GFlag==0)
            lambda_max=findLambdaMax_mt(ATb, n, k, ind, size(ind,2));
        else
            lambda_max=general_findLambdaMax_mt(ATb, n, k, G, ind, size(ind,2));
        end

        x=lambda_max;
        funVal=lambda_max;
        ValueL=lambda_max;

        return;
    end
end

% initialize a starting point
if opts.init==2
    x=zeros(n,k); c=log(m1./m2);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if ( size(x,1)~=n || size(x,2)~=k )
            error('\n Check the input .x0');
        end
    else
        x=zeros(n,k);
    end

    if isfield(opts,'c0')
        c=opts.c0;

        if ( length(c)~=k )
            error('\n Check the input .c0');
        end
    else
        c=log(m1./m2);
    end
end

% compute Ax= A * x
if (opts.nFlag==0)
    Ax=A* x;
elseif (opts.nFlag==1)
    invNu=x./repmat(nu, 1, k);  mu_invNu=mu * invNu;
    Ax=A*invNu -repmat(mu_invNu, m, 1);
else
    Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./repmat(nu, 1, k);
end

%% The main program

if (opts.mFlag==0 && opts.lFlag==0)
    % The Armijo Goldstein line search schemes + accelearted gradient descent

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1/(m*k); % the intial guess of the Lipschitz continuous gradient

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,k);
    cp=c;         ccp=zeros(1,k);

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;  sc=c + beta* ccp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute the gradient (g) at s
        As=Ax + beta* (Ax-Axp);

        % aa= - diag(y) * (A * s + sc)
        aa=- y.*(As+ repmat(sc,m,1));

        % fun_s is the logistic loss at the search point
        bb=max(aa,0);
        fun_s= sum(sum ( log( exp(-bb) +  exp(aa-bb) ) + bb ) ) / (m * k);

        % compute prob=[p_1;p_2;...;p_m]
        prob=1./( 1+ exp(aa) );

        % b= - diag(y) * (1 - prob)
        b= -y.*(1-prob) / (m * k);

        gc=sum(b); % the gradient of c

        % compute g= AT b, the gradient of x
        if (opts.nFlag==0)
            g=A'*b;
        elseif (opts.nFlag==1)
            g= A'* b -  mu' * sum(b, 1);  g=g./repmat(nu, 1, k);
        else
            invNu=b./repmat(nu, 1, k);        g=A'*invNu-mu' * sum(invNu, 1);
        end

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;
        cp=c;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the L1/Lq-norm regularized projection
            v=s-g/L; c= sc- gc/L;

            % tree overlapping group Lasso projection
            ind_work(1:2,:)=ind(1:2,:);
            ind_work(3,:)=ind(3,:) * (lambda / L);

            if (GFlag==0)
                x=altra_mt(v, n, k, ind_work, size(ind_work,2));
            else
                x=general_altra_mt(v, n, k, G, ind_work, size(ind_work,2));
            end

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute Ax= A * x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./repmat(nu, 1, k);  mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./repmat(nu, 1, k);
            end

            % aa= - diag(y) * (A * x + c)
            aa=- y.*(Ax+ repmat(c,m,1));

            % fun_x is the logistic loss at the new approximate solution
            bb=max(aa,0);
            fun_x= sum(sum ( log( exp(-bb) +  exp(aa-bb) ) + bb ) ) / (m * k);

            r_sum=(norm(v,'fro')^2 + norm(c-sc,2)^2) / 2;
            l_sum=fun_x - fun_s - sum(sum(v.* g)) - (c-sc)* gc';

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is fun_x <= fun_s + <v, g> + <c ,gc>
            %                           + L/2 * (<v,v> + <c-sc,c-sc> )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        ValueL(iterStep)=L;

        xxp=x-xp;  ccp=c-cp;

        funVal(iterStep)=fun_x;

        for i=1:n
            xRow=x(i,:);

            if (GFlag==0)
                tree_norm=treeNorm(xRow, k, ind, size(ind,2));
            else
                tree_norm=general_treeNorm(xRow, k, G, ind, size(ind,2));
            end

            funVal(iterStep)=funVal(iterStep)+ lambda*tree_norm;
        end

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=norm(xp,'fro');    norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
else
    error('\n The function does not support opts.mFlag neq 0 & opts.lFlag neq 0!');
end

author: thomaskuestner


text:
function write(hdr,data,filename);

fid=fopen(filename,'w','native');
precision=get_precision(hdr);

hdr.dim.vox_offset=352;
hdr.hist.magic='n+1';
write_hdr(fid,hdr);

fwrite(fid,ones(1,hdr.dim.vox_offset-ftell(fid)),'uint8');
data=data(:);
fwrite(fid,data(:),precision);
fclose(fid);

%%%%
function write_hdr(fid,hdr)

fwrite(fid,hdr.key.sizeof_hdr(1),'int32');

pad = zeros(1, 10-length(hdr.key.data_type));
hdr.key.data_type =[hdr.key.data_type,char(pad)];
fwrite(fid,hdr.key.data_type(1:10),'uchar');

pad = zeros(1, 18-length(hdr.key.db_name));
hdr.key.db_name = [hdr.key.db_name,char(pad)];
fwrite(fid,hdr.key.db_name(1:18),'uchar');

fwrite(fid,hdr.key.extents(1),       'int32');
fwrite(fid,hdr.key.session_error(1), 'int16');
fwrite(fid,hdr.key.regular(1),       'uchar');
fwrite(fid,hdr.key.dim_info(1),      'uchar');

fwrite(fid,hdr.dim.dim(1:8),        'int16');
fwrite(fid,hdr.dim.intent_p1(1),  'float32');
fwrite(fid,hdr.dim.intent_p2(1),  'float32');
fwrite(fid,hdr.dim.intent_p3(1),  'float32');
fwrite(fid,hdr.dim.intent_code(1),  'int16');
fwrite(fid,hdr.dim.datatype(1),     'int16');
fwrite(fid,hdr.dim.bitpix(1),       'int16');
fwrite(fid,hdr.dim.slice_start(1),  'int16');
fwrite(fid,hdr.dim.pixdim(1:8),   'float32');
fwrite(fid,hdr.dim.vox_offset(1), 'float32');
fwrite(fid,hdr.dim.scl_slope(1),  'float32');
fwrite(fid,hdr.dim.scl_inter(1),  'float32');
fwrite(fid,hdr.dim.slice_end(1),    'int16');
fwrite(fid,hdr.dim.slice_code(1),   'uchar');
fwrite(fid,hdr.dim.xyzt_units(1),   'uchar');
fwrite(fid,hdr.dim.cal_max(1),    'float32');
fwrite(fid,hdr.dim.cal_min(1),    'float32');
fwrite(fid,hdr.dim.slice_duration(1), 'float32');
fwrite(fid,hdr.dim.toffset(1),    'float32');
fwrite(fid,hdr.dim.glmax(1),        'int32');
fwrite(fid,hdr.dim.glmin(1),        'int32');

pad=zeros(1,80-length(hdr.hist.descrip));
hdr.hist.descrip=[hdr.hist.descrip,char(pad)];
fwrite(fid,hdr.hist.descrip(1:80),'uchar');

pad=zeros(1,24-length(hdr.hist.aux_file));
hdr.hist.aux_file=[hdr.hist.aux_file,char(pad)];
fwrite(fid,hdr.hist.aux_file(1:24),'uchar');

fwrite(fid,hdr.hist.qform_code,    'int16');
fwrite(fid,hdr.hist.sform_code,    'int16');
fwrite(fid,hdr.hist.quatern_b,   'float32');
fwrite(fid,hdr.hist.quatern_c,   'float32');
fwrite(fid,hdr.hist.quatern_d,   'float32');
fwrite(fid,hdr.hist.qoffset_x,   'float32');
fwrite(fid,hdr.hist.qoffset_y,   'float32');
fwrite(fid,hdr.hist.qoffset_z,   'float32');
fwrite(fid,hdr.hist.srow_x(1:4), 'float32');
fwrite(fid,hdr.hist.srow_y(1:4), 'float32');
fwrite(fid,hdr.hist.srow_z(1:4), 'float32');

pad=zeros(1,16-length(hdr.hist.intent_name));
hdr.hist.intent_name=[hdr.hist.intent_name ,char(pad)];
fwrite(fid,hdr.hist.intent_name(1:16),'uchar');

pad=zeros(1,4-length(hdr.hist.magic));
hdr.hist.magic=[hdr.hist.magic,char(pad)];
fwrite(fid,hdr.hist.magic(1:4),'uchar');

return;

%%%%
function precision=get_precision(hdr)

switch hdr.dim.datatype
    case   1,
        precision = 'ubit1';
    case   2,
        precision = 'uint8';
    case   4,
        precision = 'int16';
    case   8,
        precision = 'int32';
    case  16,
        precision = 'float32';
    case  32,
        precision = 'float32';
    case  64,
        precision = 'float64';
    case 128,
        precision = 'uint8';
    case 256
        precision = 'int8';
    case 511
        precision = 'float32';
    case 512
        precision = 'uint16';
    case 768
        precision = 'uint32';
    case 1024
        precision = 'int64';
    case 1280
        precision = 'uint64';
    case 1792,
        precision = 'float64';
    otherwise
        error('Unknown data precision.\n');
end

return;

author: yetianmed

repo: subcortex

save_path: github-repos/MATLAB/yetianmed-subcortex


text:
function v=headingSpeed2Vel(plhPoint,speed,geoEastOfNorth,angUpFromLevel)
%%HEADINGSPEED2VEL Given a target location in latitude and logitude a
%    target speed, the heading East of North and any angle it is making up
%    from level flight, obtain a 3D velocity vector. This function can be
%    useful for parameterizing
%
%INPUTS: plhPoint The 2XN or 3XN location matrix of the points in geodetic
%              latitude  and longitude in radians at which the headings are
%              taken. The point can be [latitude;longitude] or
%              [latitude;longitude;height]. The height component is ignored
%              if included because it does not change the result.
%        speed An NX1 or 1XN array of the speed of the target associated
%              with each point. If a single scalar is passed, then the
%              speed is assumed to be the same for all targets.
% geoEastOfNorth An NX1 or 1XN array of N geographic headings in radians
%              clockwise from North that should be turned into ECEF unit
%              vectors. A geographic heading is a direction in the local
%              tangent plane of an East-North-Up coordinate system as
%              defined on a specific reference ellipsoid. If all headings
%              are the same (but the angles up from level vary), then this
%              can be a scalar value. If this parameter is omitted or an
%              empty matrix is passed, then the default of 0 is used.
% angUpFromLevel An NX1 or 1XN array of N angles of the trajectory above
%              the local tangent plane to the reference ellipsoid. If all
%              elevations are the same (but geographic headings might vary),
%              then this can be a scalar value. If this parameter is
%              omitted or an empty matrix is passed, then the default of 0
%              is used.
%
%OUTPUTS: v A 3XN matrix of velocity vectors. The units are the same as the
%           units of the speed input.
%
%This function just calls geogHeading2uVec and multiplies the result by
%speed.
%
%May 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(angUpFromLevel))
    angUpFromLevel=0;
end

if(nargin<3||isempty(geoEastOfNorth))
    geoEastOfNorth=0;
end

v=bsxfun(@times,speed(:).',geogHeading2uVec(plhPoint,geoEastOfNorth,angUpFromLevel));

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% Local Regression and Likelihood, Figure 6.2.
% Author: Catherine Loader
%
% Local smooth of CO2 dataset, main trend.

load co2;
fit = locfit(year+month/12,co2,'alpha',0.5,'deg',1);
figure('Name','fig6_2: CO2 Dataset: Local smoothing, main trend' );
lfplot(fit);
xlabel('Date');
ylabel('CO2')

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function value = p07_i4 ( action, name, value )

%*****************************************************************************80
%
%% P07_I4 sets or gets I4 parameters for problem 07.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION,
%    'D' to set a parameter to its default value.
%    'S' to set a parameter,
%    'G' to get a parameter,
%    'I' to increment a parameter.
%
%    Input, string NAME, the name of the variable.
%    '#' is the number of calls to the integrand routine.
%
%    Input/output, integer VALUE.
%    * If ACTION = 'I', then VALUE is an input quantity, and is the
%      new value to be added to NAME.
%    * If ACTION = 'S', then VALUE is an input quantity, and is the
%      new value to be assigned to NAME.
%    * If ACTION = 'G', then VALUE is an output quantity, and is the
%      current value of NAME.
%
  persistent calls;

  if ( action == 'D' | action == 'd' )

    if ( name == '#' | name == '*' )
      calls = 0;
    end

  elseif ( action == 'G' | action == 'g' )

    if ( name == '#' )
      value = calls;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P07_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P07_I4 - Fatal error!' );
    end

  elseif ( action == 'I' | action == 'i' )

    if ( name == '#' )
      calls = calls + value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P07_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P07_I4 - Fatal error!' );
    end

  elseif ( action == 'S' | action == 's' )

    if ( name == '#' )
      calls = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P07_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P07_I4 - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P07_I4 - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action = "%s".\n', action );
    error ( 'P07_I4 - Fatal error!' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Score = CalWHV(PopObj,RefPoint,Weight)
% Calculate the exact weighted hypervolume value of the population

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    [N,M] = size(PopObj);
    PopObj(any(PopObj>repmat(RefPoint,N,1),2),:) = [];
    if isempty(PopObj)
        Score = 0;
    else
        pl = sortrows(PopObj);
        S  = {1,pl};
        for k = 1 : M-1
            S_ = {};
            for i = 1 : size(S,1)
                Stemp = Slice(cell2mat(S(i,2)),k,RefPoint);
                for j = 1 : size(Stemp,1)
                    temp(1) = {cell2mat(Stemp(j,1))*cell2mat(S(i,1))};
                    temp(2) = Stemp(j,2);
                    S_      = Add(temp,S_);
                end
            end
            S = S_;
        end
        Score = 0;
        for i = 1 : size(S,1)
            p     = Head(cell2mat(S(i,2)));
            Score = Score + cell2mat(S(i,1))*abs(p(M)-RefPoint(M))*Weight(i);
        end
    end
end

function S = Slice(pl,k,RefPoint)
    p  = Head(pl);
    pl = Tail(pl);
    ql = [];
    S  = {};
    while ~isempty(pl)
        ql  = Insert(p,k+1,ql);
        p_  = Head(pl);
        cell_(1,1) = {abs(p(k)-p_(k))};
        cell_(1,2) = {ql};
        S   = Add(cell_,S);
        p   = p_;
        pl  = Tail(pl);
    end
    ql = Insert(p,k+1,ql);
    cell_(1,1) = {abs(p(k)-RefPoint(k))};
    cell_(1,2) = {ql};
    S  = Add(cell_,S);
end

function ql = Insert(p,k,pl)
    flag1 = 0;
    flag2 = 0;
    ql    = [];
    hp    = Head(pl);
    while ~isempty(pl) && hp(k) < p(k)
        ql = [ql;hp];
        pl = Tail(pl);
        hp = Head(pl);
    end
    ql = [ql;p];
    m  = length(p);
    while ~isempty(pl)
        q = Head(pl);
        for i = k : m
            if p(i) < q(i)
                flag1 = 1;
            else
                if p(i) > q(i)
                    flag2 = 1;
                end
            end
        end
        if ~(flag1 == 1 && flag2 == 0)
            ql = [ql;Head(pl)];
        end
        pl = Tail(pl);
    end
end

function p = Head(pl)
    if isempty(pl)
        p = [];
    else
        p = pl(1,:);
    end
end

function ql = Tail(pl)
    if size(pl,1) < 2
        ql = [];
    else
        ql = pl(2:end,:);
    end
end

function S_ = Add(cell_,S)
    n = size(S,1);
    m = 0;
    for k = 1 : n
        if isequal(cell_(1,2),S(k,2))
            S(k,1) = {cell2mat(S(k,1))+cell2mat(cell_(1,1))};
            m = 1;
            break;
        end
    end
    if m == 0
        S(n+1,:) = cell_(1,:);
    end
    S_ = S;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Max manipulability index ALONG A LINE.
% Use stomp like to optimize along a surface/line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function experiment2B_K5_N30
close all;
global parameters

%STOMP PARAMETERS
%number of particles
K = 5;

%repeat experiment number of times
parameters.n_repeat = 3;
parameters.experiment_name = 'experiment2B_K5_N30_time.mat';
parameters.experiment_number = 2;
parameters.animate = 0;

%repeat the experiment E times
random_manips=[];
Gout = [];
ela = [];
for i=1:parameters.n_repeat
    close all
    tic;
    [pk, final_manip] = path_planning_SCO_4DOF(K);
    Gout{i}=pk;
    random_manips = [random_manips; final_manip];
    save(parameters.experiment_name)
    elapsed = toc;
    ela = [ela elapsed];
end
fprintf('TIC TOC mean: %g\n', mean(ela));
fprintf('TIC TOC cov: %g\n', cov(ela));
fprintf('TIC TOC 2sigma: %g\n', 2*sqrt(cov(ela)));

'ended'
parameters.experiment_name

repo: ARTE


text:
function [ml] = likelihood(hyp, inf, mean, cov, lik, input, target)
% Calculates the negative log marginal likelihood.
%
%% Syntax
%  ml = likelihood(hyp, inf, mean, cov, lik, input, target)
%
%% Description
% Function for calculating the negative log marginal likelihood of given
% hyperparameters, covariance function and data. If more sets of
% hyperparameters are given it outputs a vector of calculated negative log
% marginal likelihoods.
%
% Based on the work of C.E.Rasmussen.
%
% Input:
% * hyp    ... the hyperparameter struct(s), row vector of structs
% * inf    ... the function specifying the inference method
% * mean   ... the prior mean function
% * cov    ... the prior covariance function
% * lik    ... the likelihood function
% * input  ... the input part of the training data,  NxD matrix
% * target ... the output part of the training data (ie. target), Nx1 vector
%
% Output:
% * ml     ... the negative log marginal likelihood(s), vector
%
% See Also:
% gpx, minimize, covFunction, gp_initial
%
% Examples:
% gp_initial.m
%%

if nargin < 7 % input validation
  error('Too few parameters are given.'); %
end

ml=zeros(size(hyp)); % allocate

for i = 1 : size(hyp,2)
  try

    mlt = gp(hyp(i), inf, mean, cov, lik, input, target);
  catch
    mlt = +Inf;
  end

  % stability
  if (~isfinite(mlt) || isnan(mlt) || ~isreal(mlt))
    mlt = +Inf;
  end
  ml(i) = mlt;
end

author: Dynamic-Systems-and-GP

repo: GPdyn


text:
function h = clickA3DPoint(pointCloud)
%CLICKA3DPOINT
%   H = CLICKA3DPOINT(POINTCLOUD) shows a 3D point cloud and lets the user
%   select points by clicking on them. The selected point is highlighted
%   and its index in the point cloud will is printed on the screen.
%   POINTCLOUD should be a 3*N matrix, represending N 3D points.
%   Handle to the figure is returned.
%
%   other functions required:
%       CALLBACKCLICK3DPOINT  mouse click callback function
%       ROWNORM returns norms of each row of a matrix
%
%   To test this function ...
%       pointCloud = rand(3,100)*100;
%       h = clickA3DPoint(pointCloud);
%
%       now rotate or move the point cloud and try it again.
%       (on the figure View menu, turn the Camera Toolbar on, ...)
%
%   To turn off the callback ...
%       set(h, 'WindowButtonDownFcn','');
%
%   by Babak Taati
%   http://rcvlab.ece.queensu.ca/~taatib
%   Robotics and Computer Vision Laboratory (RCVLab)
%   Queen's University
%   May 4, 2005
%   revised Oct 30, 2007
%   revised May 19, 2009

if nargin ~= 1
    error('Requires one input arguments.')
end

if size(pointCloud, 1)~=3
    error('Input point cloud must be a 3*N matrix.');
end

% show the point cloud
h = gcf;
plot3(pointCloud(1,:), pointCloud(2,:), pointCloud(3,:), 'c.');
cameratoolbar('Show'); % show the camera toolbar
hold on; % so we can highlight clicked points without clearing the figure

% set the callback, pass pointCloud to the callback function
set(h, 'WindowButtonDownFcn', {@callbackClickA3DPoint, pointCloud});

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [ x, y ] = sswap ( n, x, incx, y, incy )

%*****************************************************************************80
%
%% SSWAP interchanges two vectors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 May 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vectors.
%
%    Input, real X(*), one of the vectors to swap.
%
%    Input, integer INCX, the increment between successive entries of X.
%
%    Input, real Y(*), one of the vectors to swap.
%
%    Input, integer INCY, the increment between successive elements of Y.
%
%    Output, real X(*), the swapped vector.
%
%    Output, real Y(*), the swapped vector.
%
  if ( n <= 0 )

  elseif ( incx == 1 & incy == 1 )

    m = mod ( n, 3 );

    for i = 1 : m
      temp = x(i);
      x(i) = y(i);
      y(i) = temp;
    end

    for i = m+1 : 3 : n

      temp = x(i);
      x(i) = y(i);
      y(i) = temp;

      temp = x(i+1);
      x(i+1) = y(i+1);
      y(i+1) = temp;

      temp = x(i+2);
      x(i+2) = y(i+2);
      y(i+2) = temp;

    end

  else

    if ( 0 <= incx )
      ix = 1;
    else
      ix = ( - n + 1 ) * incx + 1;
    end

    if ( 0 <= incy )
      iy = 1;
    else
      iy = ( - n + 1 ) * incy + 1;
    end

    for i = 1 : n
      temp = x(ix);
      x(ix) = y(iy);
      y(iy) = temp;
      ix = ix + incx;
      iy = iy + incy;
    end

  end

  return
end


text:
classdef InverseGammaD
%%INVERSEGAMMAD Functions to handle the inverse gamma distribution.
%Implemented methods are: mean, var, PDF, CDF, rand, entropy
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)

function val=mean(k, beta)
%%MEAN Obtain the mean of the inverse gamma distribution for given shape
%      and scale parameters.
%
%INPUTS: k The shape parameter of the distribution; k>0.
%     beta The scale parameter of the distribution; beta>0.
%
%OUTPUTS: val The mean of the inverse gamma distribution.
%
%If k<1, then the distribution has no mean and a NaN is returned.
%
%The distribution arises when estimating radar cross sections in Ch.
%10.3.2 of [1].
%
%REFERENCES:
%[1] W. Koch, Tracking and Sensor Fusion: Methodological Framework and
%    Selected Applications. Berlin, Germany: Springer, 2014.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(k<1)
        val=NaN;
    else
        val=beta/(k-1);
    end
end

function val=var(k,beta)
%%VAR Obtain the variance of the inverse gamma distribution for given shape
%     and scale parameters.
%
%INPUTS: k The shape parameter of the distribution; k>0.
%     beta The scale parameter of the distribution; beta>0.
%
%OUTPUTS: val  The variance of the inverse gamma distribution.
%
%If k<2, then the distribution has no variance and a NaN is returned.
%
%The distribution arises when estimating radar cross sections in Ch.
%10.3.2 of [1].
%
%REFERENCES:
%[1] W. Koch, Tracking and Sensor Fusion: Methodological Framework and
%    Selected Applications. Berlin, Germany: Springer, 2014.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(k<2)
        val=NaN;
    else
        val=beta^2/((k-1)^2*(k-2));
    end
end

function val=PDF(x,k,beta)
%%PDF Evaluate the probability density function (PDF) of the inverse gamma
%     distribution at the desired points.
%
%INPUTS: x The point(s) at which the inverse gamma PDF is to be evaluated.
%          Note that x>=0.
%        k The shape parameter of the distribution; k>0.
%     beta The scale parameter of the distribution; beta>0.
%
%OUTPUTS: val The value(s) of the PDF of the inverse gamma distribution
%             with parameters k and beta evaluated at x.
%
%The inverse gamma distribution is the distribution of the inverse of a
%gamma random variable. Thus, with a little bit of algebra, it can be
%expressed in terms of the gamma PDF.
%
%The distribution arises when estimating radar cross sections in Ch.
%10.3.2 of [1].
%
%REFERENCES:
%[1] W. Koch, Tracking and Sensor Fusion: Methodological Framework and
%    Selected Applications. Berlin, Germany: Springer, 2014.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=(beta^k/gamma(k)).*x.^(-k-1).*exp(-beta./x);
end

function val=CDF(x,k,beta)
%%CDF Evaluate the cumulative distribution function (CDF) of the inverse
%     gamma distribution at desired points.
%
%INPUTS: x The point(s) at which the inverse gamma CDF is to be  evaluated.
%          Note that x>=0.
%        k The shape parameter of the distribution; k>0.
%     beta The scale parameter of the distribution; beta>0.
%
%OUTPUTS: prob The value(s) of the CDF of the inverse gamma distribution
%              with parameters k and beta evaluated at x.
%
%The inverse gamma distribution is the distribution of the inverse of a
%gamma random variable. Thus, with a little bit of algebra, it can be
%expressed in terms of the gamma PDF.
%
%The distribution arises when estimating radar cross sections in Ch.
%10.3.2 of [1].
%
%REFERENCES:
%[1] W. Koch, Tracking and Sensor Fusion: Methodological Framework and
%    Selected Applications. Berlin, Germany: Springer, 2014.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=1-GammaD.CDF(1./x,k,1/beta);
end

function vals=rand(N,k,beta)
%%RAND Generate inverse gamma distributed random variables with the given
%      parameters.
%
%INPUTS: N If N is a scalar, then rand returns an NXN matrix of random
%          variables. If N=[M,N1] is a two-element row  vector, then rand
%          returns an MXN1 matrix of random variables.
%        k The shape parameter of the distribution; k>0.
%     beta The scale parameter of the distribution; beta>0.
%
%OUTPUTS: vals A matrix whose dimensions are determined by N of the
%              generated inverse gamma random variables.
%
%This just takes the inverse of gamma distributed random variables.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    vals=1./randGamma(N,k,1/beta);
end

function entropyVal=entropy(k,beta)
%%ENTROPY Obtain the differential entropy of the inverse gamma
%         distribution given in nats. The differential entropy of a
%         continuous distribution is entropy=-int_x p(x)*log(p(x)) dx where
%         the integral is over all values of x. Units of nats mean that the
%         natural logarithm is used in the definition. Unlike the Shannon
%         entropy for discrete variables, the differential entropy of
%         continuous variables can be both positive and negative.
%
%INPUTS: k The shape parameter of the distribution; k>0.
%     beta The scale parameter of the distribution; beta>0.
%
%OUTPUTS: entropyVal The value of the differential entropy in nats.
%
%Differential entropy is defined in Chapter 8 of [1].
%
%REFERENCES:
%[1] T. M. Cover and J. A. Thomas, Elements of Information Theory, 2nd ed.
%    Hoboken, NJ: Wiley-Interscience, 2006.
%
%April 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

    entropyVal=k+log(beta)+gammaln(k)-(1+k)*psi(k);
end

end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
function G = rayxgridfZ(F, transla, val)
%------------------------------------------------------------------------------
% Multiplies gridfunction F with element of stencil given by its value
% val and its position transla.
% Excess area is filled by padding with zeroes.
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: June 23, 2000.
% (c) 1999-2002 Stichting CWI, Amsterdam.
%------------------------------------------------------------------------------
o=[0 0];
if ~all(size(o) == size(transla))
  error(' rayxgridfZ - unexpected dimensions of transla ')
end
G = val * moveLRZ(moveUDZ(F, -transla(1)), -transla(2));
%------------------------------------------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function [irf_record_allt]=tvbvarirf4(beta_gibbs,sigma_t_gibbs,It,Bu,IRFperiods,n,m,p,k,T,signresperiods,signrestable)

% now identify all the periods concerned with restrictions
% first expand the non-empty entries in signresperiods since they are only expressed in intervals: transform into list
% for instance, translate [1 4] into [1 2 3 4]; I don't think this can done without a loop
temp=cell2mat(signresperiods(~cellfun(@isempty,signresperiods)));
periods=[];
for ii=1:size(temp,1)
    periods=[periods temp(ii,1):temp(ii,2)];
end
% suppress duplicates and sort
periods=sort(unique(periods))';
% count the total number of restriction periods (required for IRF matrix)
nperiods=size(periods,1);

% Identify the restriction matrices
% create five cells, corresponding to the three possible restrictions:
% one cell for sign restrictions, three cells for magnitude restrictions, one cell for zero restrictions
Scell=cell(1,n);
Mcell=cell(1,n);
Mlcell=cell(1,n);
Mucell=cell(1,n);
Zcell=cell(1,n);

% loop over rows and columns of the period matrix
for ii=1:n
    for jj=1:n
        % if entry (ii,jj) of the period matrix is not empty...
        if ~isempty(signresperiods{ii,jj})
            % ... then there is a restriction over one (or several) periods
            % loop overt those periods
            for kk=signresperiods{ii,jj}(1,1):signresperiods{ii,jj}(1,2)
                % identify the position of the considered period within the list of all periods (required to build the matrix)
                position=find(periods==kk);
                % now create the restriction matrix: this will depend on the type of restriction
                % if it is a positive sign restriction...
                if strcmp(signrestable{ii,jj},'+')
                    % ... then input a 1 entry in the corresponding S matrix
                    Scell{1,jj}=[Scell{1,jj};zeros(1,n*nperiods)];
                    Scell{1,jj}(end,(position-1)*n+ii)=1;
                    % if it is a negative sign restriction...
                elseif strcmp(signrestable{ii,jj},'-')
                    % ... then input a -1 entry in the corresponding S matrix
                    Scell{1,jj}=[Scell{1,jj};zeros(1,n*nperiods)];
                    Scell{1,jj}(end,(position-1)*n+ii)=-1;
                    % if it is a zero restriction...
                elseif strcmp(signrestable{ii,jj},'0')
                    % ... then input a 1 entry in the corresponding Z matrix
                    Zcell{1,jj}=[Zcell{1,jj};zeros(1,n*nperiods)];
                    Zcell{1,jj}(end,(position-1)*n+ii)=1;
                    % else, a non-empty entry being neither a sign nor a zero restriction has to be a magnitude restriction
                else
                    % fill the corresponding M matrices:
                    % input a 1 in M
                    Mcell{1,jj}=[Mcell{1,jj};zeros(1,n*nperiods)];
                    Mcell{1,jj}(end,(position-1)*n+ii)=1;
                    % input the lower value of the interval in Ml
                    temp=str2num(signrestable{ii,jj});
                    Mlcell{1,jj}=[Mlcell{1,jj};temp(1,1)];
                    % input the upper value of the interval in Mu
                    Mucell{1,jj}=[Mucell{1,jj};temp(1,2)];
                end
            end
        end
    end
end

% now check what kind of restrictions apply among sign, zero and magnitude restrictions
% check for sign restrictions: if there are any, at least one entry in the cell Scell is non-empty
if sum(~cellfun(@isempty,Scell))~=0
    signres=1;
else
    signres=0;
end
% similarly check for zero restrictions
if sum(~cellfun(@isempty,Zcell))~=0
    zerores=1;
else
    zerores=0;
end
% and finally, check for magnitude restrictions
if sum(~cellfun(@isempty,Mcell))~=0
    magnres=1;
else
    magnres=0;
end

% create the cell aray that will store the values from the simulations
irf_record_allt=cell(n,n);
storage1=cell(It-Bu,T);

% loop over sample periods
for tt=1:T

    % loop over iterations
    parfor ii=1:It-Bu
        % initiate the variable 'success'; this variable will be used to check whether the restrictions are satisfied
        % if there are only zero restrictions, they will be satisfied by construction, and 'success' will simply be ignored
        success=0;
        % how the algorithm will be conducted will depend on the types of restrictions implemented

        % if there are only zero restrictions, the algorithm is simple as no checking is required: the conditions are satisfied by construction
        if zerores==1 && signres==0 && magnres==0
            % draw beta and sigma
            beta=beta_gibbs{tt,1}(:,ii);
            sigma=sigma_t_gibbs{tt,1}(:,:,ii);
            D=chol(bear.nspd(sigma),'lower');
            % obtain orthogonalised IRFs
            [~,ortirfmatrix]=bear.irfsim(beta,D,n,m,p,k,max(IRFperiods,max(periods)));
            % generate the stacked IRF matrix
            stackedirfmat=[];
            for jj=1:numel(periods)
                stackedirfmat=[stackedirfmat;ortirfmatrix(:,:,periods(jj,1)+1)];
            end
            % draw an entire random matrix Q satisfying the zero restrictions
            [Q]=bear.qzerores(n,Zcell,stackedirfmat);
            % there is no need to verify the restrictions: there are satisfied by construction

            % if there are sign/magnitude restrictions, possibly associated with zero restrictions
        else
            % the algorithm becomes a bit more complicated as conditions now need to be checked
            % to maintain efficiency, the algorithm proceeds recursively shock by shock, and stops as soon as a condition on the considered shock fails
            % repeat algorithm for the iteration as long as not all conditions are satisfied
            while success==0
                % switch 'success' to 1; it will be turned back to zero if at any time Q is detected as a candidate not satisfying the restrictions
                success=1;
                % draw randomly the vector of VAR coefficients: draw a random index
                index=floor(rand*(It-Bu))+1;
                % draw beta and sigma
                beta=beta_gibbs{tt,1}(:,index);
                sigma=sigma_t_gibbs{tt,1}(:,:,index);
                D=chol(bear.nspd(sigma),'lower');
                % obtain orthogonalised IRFs
                [~,ortirfmatrix]=bear.irfsim(beta,D,n,m,p,k,max(IRFperiods,max(periods)));
                % generate the stacked IRF matrix
                stackedirfmat=[];
                for jj=1:numel(periods)
                    stackedirfmat=[stackedirfmat;ortirfmatrix(:,:,periods(jj,1)+1)];
                end
                % initiate Qj
                Qj=[];
                % now start looping over the shocks and checking sequentially whether conditions on these shocks hold
                % stop as soon as one restriction fails
                mm=1;
                while success==1 && mm<=n
                    % build column j of the random matrix Q
                    [qj]=bear.qrandj(n,Zcell{1,mm},stackedirfmat,Qj);
                    % obtain the candidate column fj
                    fj=stackedirfmat*qj;
                    % check restrictions: first sign restrictions
                    [success, qj]=bear.checksignres(Scell{1,mm},qj,fj);
                    % if 'success' is still equal to 1, also check for magnitude restrictions
                    if success==1
                        [success]=bear.checkmagres(Mcell{1,mm},Mlcell{1,mm},Mucell{1,mm},fj);
                    end
                    % also, if 'success' is still equal to 1, update Qj by concatenating qj
                    if success==1
                        Qj=[Qj qj];
                    end
                    mm=mm+1;
                end
                % repeat this loop until a succesful draw is obtained
            end
            % with succesful Qj at hand, eventually set Q as Qj
            Q=Qj;
        end

        % store
        for jj=1:IRFperiods
            storage1{ii,tt}(:,:,jj)=ortirfmatrix(:,:,jj)*Q;
        end

    end

end

% reorganise storage
% loop over iterations
for tt=1:T
    for ii=1:It-Bu
        % loop over IRF periods
        for jj=1:IRFperiods
            % loop over variables
            for kk=1:n
                % loop over shocks
                for ll=1:n
                    irf_record_allt{kk,ll}(ii,jj,tt)=storage1{ii,tt}(kk,ll,jj);
                end
            end
        end
    end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function plotOcclusionBoundaries(bndinfo, blabels)
% im = drawOcclusionBoundaries(im, bndinfo, blabels)
% blabels(i) = (0, 1, or 2) for off, left, or right
%

if numel(blabels)==bndinfo.ne*2
    blabels = blabels(1:end/2) + 2*blabels(end/2+1:end);
end

hold on
imsize = bndinfo.imsize;

color = [0.01 0.01 0.01];

indices = bndinfo.edges.indices;

arrowdist = ceil(sqrt(imsize(1).^2 + imsize(2).^2)/10);

for k = 1:numel(indices)

    if blabels(k)>0

        ind = double(indices{k});
        [ey, ex] = ind2sub(bndinfo.imsize(1:2), ind);
        npix = numel(ind);

        narrows = ceil(npix/arrowdist);

        epos = ceil((1:narrows) / (narrows+1) * npix);

        for j = 1:numel(epos)

            [ay, ax] = ind2sub(imsize(1:2), ind(epos(j)));

            if blabels(k)==1
                [y1, x1] = ind2sub(imsize(1:2), ind(max(epos(j)-10,1)));
                [y2, x2] = ind2sub(imsize(1:2), ind(min(epos(j),npix)));
            else % blabels(k)==2;
                [y1, x1] = ind2sub(imsize(1:2), ind(min(epos(j)+10,npix)));
                [y2, x2] = ind2sub(imsize(1:2), ind(min(epos(j),npix)));
            end
            theta = atan2(y2-y1, x2-x1);
%             if blabels(k)==2
%                 theta = mod(theta+pi, 2*pi);
%             end

            asx = ax - cos(theta);
            asy = ay - sin(theta);

            %[ax,ay] = dsxy2figxy(gca, ax, ay);
            %[asx, asy] = dsxy2figxy(gca, asx, asy);
            ax = (ax-1)/(imsize(2)-1);   ay = 1-(ay-1)/(imsize(1)-1);
            asx = (asx-1)/(imsize(2)-1);   asy = 1-(asy-1)/(imsize(1)-1);

            plot(ex, ey, 'Color', 1-color, 'LineWidth', 3);
            plot(ex, ey, 'Color', color, 'LineWidth', 1);

            annotation('arrow', [asx ax], [asy ay], 'LineStyle', 'none', ...
                'HeadWidth', 17, 'HeadLength', 10);
        end
    end
end

% arrhsv = ones(5, 3);
% arrhsv(:, 1) = [100 0 135 169 42]/255;
%
% colors = max(hsv2rgb(arrhsv), 1/255); % avoid pixel intensity of 0
%
% [imh, imw] = size(bndinfo.wseg);
% edgeim = zeros([imh imw 3]);
%
% hasarrow = false(bndinfo.ne*2, 1);
% arrowpos = zeros(1000, 2);
% narrow = 0;
%
% for k = 1:numel(blabels)
%     if blabels(k) > 0
%
%         ku = mod(k-1, bndinfo.ne)+1;
%
%         cols = colors(blabels(k), :);
%         pix = bndinfo.edges.indices{ku};
%         npix = size(edgeim, 1)*size(edgeim, 2);
%         for b = 1:3
%             edgeim(pix + (b-1)*npix) = cols(b);
%         end
%
%         jcts = bndinfo.edges.junctions(ku, :);
%         if k~=ku
%             jcts = jcts([2 1]);
%         end
%
%         x = bndinfo.junctions.position(jcts, 1);
%         y = bndinfo.junctions.position(jcts, 2);
%
%         if ~any(sum((arrowpos(1:narrow, :) - repmat([x(2) y(2)], narrow, 1)).^2, 2) < 20^2) ...
%                 && (numel(pix)>10)
%
%             narrow = narrow + 1;
%             arrowpos(narrow, :) = [x(2) y(2)];
%
%             arrow = struct('x', x(2), 'y', y(2), ...
%                 'angle', atan2(-(y(2)-y(1)), x(2)-x(1)), ...
%                 'radius', 0, 'head_length', round(sqrt(npix)/50), ...
%                 'head_base_angle', 30);
%             arrow.angle = mod(arrow.angle/pi*180, 360);
%             edgeim = draw_arrow_image(edgeim, arrow, cols);
%         end
%
%     end
% end
%
% sz = ceil(sqrt(size(im, 1).^2 + size(im,2).^2) / 500);
%
% for b = 1:3
%     edgeim(:, :, b) = ordfilt2(edgeim(:, :, b), sz*sz, ones(sz, sz));
% end
%
% if strcmp(class(im), 'uint8')
%     edgeim = im2uint8(edgeim);
% end
%
% im(edgeim>0) = edgeim(edgeim>0);
%
% % hold off, clf, imagesc(im), axis image, hold on
% % drawnow;
%
% % dy = -1;
% % dx = -1;
% %
% % for k = 1:numel(blabels)
% %     if blabels(k) > 0
% %         ku = mod(k-1, bndinfo.ne)+1;
% %         jcts = bndinfo.edges.junctions(ku, :);
% %         if k~=ku
% %             jcts = jcts([2 1]);
% %         end
% %
% %         x = min(max((bndinfo.junctions.position(jcts, 1)+dx),1), imw);
% %         y = min(max((bndinfo.junctions.position(jcts, 2)+dy), 1), imh);
% %         [arrowx,arrowy] = dsxy2figxy(gca, x, y);
% %
% %         if mod(k, 5)==0
% %             annotation('arrow', arrowx, 1-arrowy, ...
% %                 'Color', colors(blabels(k), :), 'HeadStyle', 'vback2', 'LineStyle', 'none');
% %         end
% %
% %     end
% % end
%
%

author: Cloud-CV

repo: object-proposals


text:
function plotData(X, y)
%PLOTDATA Plots the data points X and y into a new figure
%   PLOTDATA(x,y) plots the data points with + for the positive examples
%   and o for the negative examples. X is assumed to be a Mx2 matrix.

% Create New Figure
figure; hold on;
% Find Indices of Positive and Negative Examples
pos = find(y==1); neg = find(y == 0);
% Plot Examples
plot(X(pos, 1), X(pos, 2), 'k+','LineWidth', 2, ...
'MarkerSize', 7);
plot(X(neg, 1), X(neg, 2), 'ko', 'MarkerFaceColor', 'y', ...
'MarkerSize', 7);
% ====================== YOUR CODE HERE ======================
% Instructions: Plot the positive and negative examples on a
%               2D plot, using the option 'k+' for the positive
%               examples and 'ko' for the negative examples.
%

% =========================================================================

hold off;

end

author: AvaisP


text:
%generalized_sylvester_direct2 is a function.
%    [Tzz, retcode] = generalized_sylvester_direct2(A, B, C, fv_lamba_bf, W)
%

author: jmaih


text:
% Liou-Steffen scheme with limiting: Eulerstep in Runge-Kutta time stepping
% Called by Liou_Steffen_MUSCL_scheme

% Theory in Sections 10.3 and 10.8 of:

% 	P. Wesseling: Principles of Computational Fluid Dynamics
% 	Springer, Heidelberg, 2000 ISBN 3-5453-0. XII, 642 pp.
% 	See http://ta.twi.tudelft.nl/nw/users/wesseling/cfdbook.html

% Function called: extrap

  [U1L, U1R] = extrap(rhostar,limtype);   	% Left and right
  [U2L, U2R] = extrap(mstar,limtype);		%extrapolated
  [U3L, U3R] = extrap(totenstar,limtype);	%states
  uL = U2L./U1L; uR = U2R./U1R;
  cL = sqrt(gamma*gam1*(U3L./U1L - 0.5*uL.^2));
  cR = sqrt(gamma*gam1*(U3R./U1R - 0.5*uR.^2));
  pL = U1L.*cL.^2/gamma; pR = U1R.*cR.^2/gamma;
  machL = uL./cL;  machR = uR./cR;
  rhocL = U1L.*cL; rhocR = U1R.*cR;
  enthalpyL = cL.^2/(gamma-1) + 0.5*uL.^2;
  enthalpyR = cR.^2/(gamma-1) + 0.5*uR.^2;

% Preallocations
  machplus = machL; machminus = machL; machhalf = machL;
  presplus = machL; presminus = machL;

  for j = 1:J-1
    if abs(machL(j)) > 1
      machplus(j) = 0.5*(machL(j) + abs(machL(j)));
      presplus(j) = 0.5*pL(j)*(1 + abs(machL(j))/machL(j));
    else
      machplus(j) = 0.25*(machL(j) + 1)^2;
      presplus(j) = 0.5*pL(j)*(1 + machL(j));
    end
    if abs(machR(j)) > 1
      machminus(j) = 0.5*(machR(j) - abs(machR(j)));
      presminus(j) = 0.5*pR(j)*(1 - abs(machR(j))/machR(j));
    else
      machminus(j) = - 0.25*(machR(j) - 1)^2;
      presminus(j) = 0.5*pR(j)*(1 - machR(j));
    end
  end

% Liou-Steffen fluxes
  machhalf = machplus + machminus;
  flux1a = 0.5*(machhalf + abs(machhalf)).*rhocL;
  flux1b = 0.5*(machhalf - abs(machhalf)).*rhocR;
  flux1  = flux1a + flux1b;
  flux2  = flux1a.*uL + flux1b.*uR + presplus + presminus;
  flux3  = flux1a.*enthalpyL + flux1b.*enthalpyR;

% Update of state variables
  rhostar(1) = rholeft; rhostar(J) = rhoright;
  mstar(1) = rholeft*uleft; mstar(J) = rhoright*uright;
  totenstar(1) = totenleft; totenstar(J) = totenright;
  for j = 2:J-1
    rhostar(j)   = rhoold(j)   - rkalpha*lambda*(flux1(j) - flux1(j-1));
    mstar(j)     = mold(j)     - rkalpha*lambda*(flux2(j) - flux2(j-1));
    totenstar(j) = totenold(j) - rkalpha*lambda*(flux3(j) - flux3(j-1));
  end

repo: Aero-matlab


text:
function [res] = tt_axpy2(log_a, sign_a, x, log_p, sign_p, y, eps, max_rank)
%Returns A*X+P*Y in the TT-format (stabilized version)
%   [RES]=TT_AXPY(A,X,P,Y,EPS,MAX_RANK) Stabilized a*X+p*Y where  log_a =
%   log(abs(a)), sign_a = sign(a),  log_p = log(abs(p)), sign_p = sign(p),
%   EPS, MAX_RANK - compression parameters

if ((nargin<8)||(isempty(max_rank)))
    max_rank=[];
end;

ax = tt_scal2(x, log_a, sign_a);
py = tt_scal2(y, log_p, sign_p);

res = tt_add(ax,py);
res = tt_compr2(res, eps, max_rank);

end

author: oseledets

repo: TT-Toolbox


text:
function brnd = betarand(a,b,r,c)
%  BETARAND Random matrices from beta distribution.
%     R = BETARAND(A,B) returns an array of random numbers chosen from the
%     beta distribution with parameters A and B.  The size of R is the common
%     size of A and B if both are arrays.  If either parameter is a scalar,
%     the size of R is the size of the other parameter.
%
%     R = BETARAND(A,B,M,N) returns an M-by-N matrix.

% Copyright (c) 1995-1997, 2005-2007 Kurt Hornik

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.
if (nargin > 1)
  if (~isscalar(a) || ~isscalar(b))
    if (size(a,1) ~= size(b,1))
      error ('betarnd: a and b must be of common size or scalar');
    end
  end
end

if (nargin == 4)
  if (~(isscalar(r) && (r > 0) && (r == round(r))))
    error ('betarnd: r must be a positive integer');
  end
  if (~(isscalar (c) && (c > 0) && (c == round (c))))
    error ('betarnd: c must be a positive integer');
  end
  sz = [r, c];

  if (any (size (a) ~= 1) && (length (size (a)) ~= length (sz) || any (size (a) ~= sz)))
    error ('betarnd: a and b must be scalar or of size [r,c]');
  end
elseif (nargin == 3)
  if (isscalar (r) && (r > 0))
    sz = [r, r];
  elseif (isvector(r) && all (r > 0))
    sz = r(:)';
  else
    error ('betarnd: r must be a positive integer or vector');
  end

  if (any (size (a) ~= 1) && (length (size (a)) ~= length (sz) || any (size (a) ~= sz)))
    error ('betarnd: a and b must be scalar or of size sz');
  end
elseif (nargin == 2)
  sz = size(a);
else
  error('must provide atleast 2 parameters')
end

if (isscalar(a) && isscalar(b))
  if (find ((a < 0) | isinf(a) | (b < 0) | isinf(b)))
    brnd = NaN * ones (sz);
  else
    r1 = gamrand(a,2.*a,sz(1),sz(2));
    brnd = r1 ./ (r1 + gamrand(b,2.*b,sz(1),sz(2)));
  end
else
  brnd = zeros (sz);

  k = find ((a < 0) | isinf(a) | (b < 0) | isinf(b));
  if (any (k))
    brnd(k) = NaN * ones (size (k));
  end

  k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf));
  if (any (k))
    r1 = gamrnd(a(k),1,size(k,1),size(k,2));
    brnd(k) = r1 ./ (r1 + gamrnd(b(k),1,size(k,1),size(k,2)));
  end
end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function recon_im = SENSE_recon(smap, reduced_fft, reduced_dim, np, regularizer, figs_on)
% function recon_im = SENSE_recon(smap, reduced_fft, reduced_dim, np,
%       regularizer, figs_on)
% in:
%   smap          [N M nc]                complex sensitivity maps
%   reduced_fft   [N/np M] or [N M/np]    undersampled fft
%   reduced_dim   1 or 2                  dimension of undersampling
%   np            degree of undersampling
%   regularizer   string, 'none' or 'tikhonov'
%   figs_on       boolean that shows intermediate images for debugging
% out:
%   recon_im      [N M]                   SENSE reconstructed image
% This implementation assumes that the fft will be reduced in such a way
% that the DC term is still included (i.e. even DFT coefficients retained)
% 2012-06-07 Mai Le, University of Michigan

% dimensions of original image
dims = [size(smap,1) size(smap,2)];

% number of coils
nc = size(smap,3);

% create aliased intermediate images
aliased_im = zeros(size(reduced_fft,1), size(reduced_fft,2), nc);
for ii = 1:nc
     aliased_im(:,:,ii) = ifft2(reduced_fft(:,:,ii));
end

% plot aliased images
if (figs_on)
    figure;
    for ii = 1:nc
        subplot(2,2,ii); imshow(abs(aliased_im(:,:,ii)),[]); colorbar;
    end
end

% calculate regularization parameter
% I found empirically that somewhere between 0.5% and 1% of the max SVD
% value works well for beta
C = Cdiff1(np);
[U, SIG, V] = svd(C'*C);
beta = 0.01*max(SIG(:));
beta = 0.005*max(SIG(:));

% reconstruct image by finding (regularized) LS solution for each set of pixels
recon_im = zeros(dims);
if (reduced_dim == 1)
    for ii = 1:dims(1)/np
        for jj = 1:dims(2)
            S = constructS(reduced_dim, smap, ii, jj, np);
            a = squeeze([aliased_im(ii,jj,1:nc)]);
            v = recon_pixels(S,a,regularizer, beta);
            recon_im = place_pixels(recon_im, v, reduced_dim, dims, ii, jj, np);
        end
    end
else
    for ii = 1:dims(1)
        for jj = 1:dims(2)/np
            S = constructS(reduced_dim, smap, ii, jj, np);
            a = squeeze([aliased_im(ii,jj,1:nc)]);
            v = recon_pixels(S,a,regularizer, beta);
            recon_im = place_pixels(recon_im, v, reduced_dim, dims, ii, jj, np);
        end
    end
end

% selects values from sensitivity map relevant for specific pixels
function S = constructS(reduced_dim, smap, ii, jj, np)
dims = size(smap);
indeces = (0:np-1)*dims(reduced_dim)/np;
if (reduced_dim == 1)
    S = squeeze([smap(ii+indeces,jj,:)]).';
else
    S = squeeze([smap(ii,jj+indeces,:)]).';
end

% places each of the reconstructed pixels in the reconstructed image
function recon_im = place_pixels(recon_im, v, reduced_dim, dims, ii, jj, np)
indeces = (0:np-1)*dims(reduced_dim)/np;
if (reduced_dim == 1)
    for kk = 1:np
        recon_im(ii+indeces(kk),jj) = v(kk);
    end
else
    for kk = 1:np
        recon_im(ii,jj+indeces(kk)) = v(kk);
    end
end

% reconstruct pixels with or without Tikhonov regularization
function pixels = recon_pixels(S,a,regularizer, beta)
switch regularizer
    case 'none'
        pixels = S\a;
    case 'tikhonov'
        M = min(size(S));
        pixels = inv(S'*S+(beta^2)*eye(M))*S'*a;
    otherwise
        pixels = S\a;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function demo_designOptimization ()
% // VBA toolbox //////////////////////////////////////////////////////////
%
% [posterior, out] = demo_designOptimization ()
% demo of off-line and online design optimisation
%
% This demo simulates a psychophysics paradigm similar to a signal
% detection task, whereby the detection probability is a sigmoidal function
% of the stimulus contrast (which is the design control variable).
% Our goal is to estimate  the inflexion point (detection threshold) and the
% sigmoid steepness (d prime) or the response function.
% In order to provide the most efficient estimate of these model parameters,
% we will show how to use offline (before the experiment) and online
% (during the experiment) design optimization to decide given trial-by-trial
% subjects' binary choice data (seen/unseen) which stimulus to show next.
%
% /////////////////////////////////////////////////////////////////////////

%% Global values
% =========================================================================

% number of trials for the respective simulations
N = 50;

% true parameter values for the simulations
phi = [- .5; 2.5]; % simulated parameters: [inflexion point, log-slope]

% range of potential simuli
uRange = linspace(- 1, 1, N);

% prepare display
VBA_figure('Name', 'demo_designOptimisation');

%% Define the model
% =========================================================================

% observation function
% -------------------------------------------------------------------------
function [gx, dgdx, dgdp] = g_psychometric (~, phi, u, ~)
    % tip: VBA_sigmoid returns derivatives wrt parameters in alphabetical
    % order, ie. center -> slope.
    [gx, ~, dgdp] = VBA_sigmoid (u,...
    'center', phi(1), ...
    'slope', exp (phi(2)) ...
    );
    dgdp(2,:) = dgdp(2,:) * exp (phi(2));
    dgdx = [];
end

% dimensions
% -------------------------------------------------------------------------
dim.n_phi = 2;
dim.n_t = 1;
dim.p = N;

% general options
% -------------------------------------------------------------------------
% binary observations
options.sources.type = 1;

% no display
options.DisplayWin = 0;
options.verbose = 0;

%% 1) No optimisation
% =========================================================================
% here, we'll use the most naive approach, a full swipe of all possible
% stimulus intensities

% experimental design
% -------------------------------------------------------------------------
u = uRange';

% simulate responses
% -------------------------------------------------------------------------
y = VBA_simulate (1,[],@g_psychometric,[],phi,u,[],[],options);

% estimate parameters
% -------------------------------------------------------------------------
posterior_naive = VBA_NLStateSpaceModel (y,u,[],@g_psychometric,dim,options);

% display results
% -------------------------------------------------------------------------
plot_design(1, 'no optimisation', u, y, posterior_naive);

%% 2) Offline optimisation
% =========================================================================
% Here, we will try to find a better design before running the experiment

% experimental design
% -------------------------------------------------------------------------
% number of designs to try
nAttempts = 1e4;

% initialization
fprintf('Offline optimisation: optimizing (  0%%)');
efficiency_offOpt = - Inf;
efficiencyDesign = nan(1, nAttempts);
keepDesign = [1];

% loop over designs
for attempt = 1 : nAttempts

    % draw random design
    u_attempt = uRange(randi (numel (uRange), 1, N))';

    % estimate efficiency
    efficiencyDesign(attempt) = VBA_designEfficiency([],@g_psychometric,dim,options,u_attempt,'parameters');

    % if better, store and display
    if efficiencyDesign(attempt) > efficiency_offOpt
        efficiency_offOpt = efficiencyDesign(attempt);
        u = sort(u_attempt);
        keepDesign(end+1) = attempt;
    end

    if efficiencyDesign(attempt) > efficiency_offOpt || mod(attempt, 50) == 0
        plot_design(2, 'offline optimisation', u, y, [], [], efficiencyDesign(keepDesign));
    end

    % progress bar
    fprintf('\b\b\b\b\b%3d%%)', round(100* attempt / nAttempts));
end

% simulate responses
% -------------------------------------------------------------------------
y = VBA_simulate (1,[],@g_psychometric,[],phi,u,[],[],options);

% estimate parameters
% -------------------------------------------------------------------------
posterior_offline = VBA_NLStateSpaceModel (y,u,[],@g_psychometric,dim,options);

% display results
% -------------------------------------------------------------------------
plot_design(2, 'offline optimisation', u, y, posterior_offline);

%% 3) Online optimisation
% =========================================================================
% Here, we will optimize the design during the experiment by taking into
% account trial-by-trial subject's responses to adaptively select the next
% best stimulus to present

% initialization
opt = options;
u = nan(N, 1);
efficiencyInput = nan(1, length (uRange));
efficiencyDesign = nan(1, length (uRange));

% run experiment
for t = 1 : N

    % extend design
    % ---------------------------------------------------------------------
    % start from current posterior belief
    try
        opt.priors = posterior_online;
    end

    % compute efficiency of potential stimuli
    dim.p = 1;
    for i = 1 : length (uRange)
        efficiencyInput(i) = VBA_designEfficiency([],@g_psychometric,dim,opt,uRange(i),'parameters');
    end

    % find best next stimulus to present
    [efficiencyDesign(t), idxMaxEff] = max (efficiencyInput);
    u(t) = uRange(idxMaxEff);

    % simulate 1 responses
    % ---------------------------------------------------------------------
    y(t) = VBA_simulate (1,[],@g_psychometric,[],phi,u(t),Inf,[],options);

    % estimate parameters given data acquired so far
    % ---------------------------------------------------------------------
    dim.p = t;
    posterior_online = VBA_NLStateSpaceModel(y(1:t),u(1:t),[],@g_psychometric,dim,options);

    % display
    % ---------------------------------------------------------------------
    plot_design(3, 'online optimisation', u, y, posterior_online, efficiencyInput, efficiencyDesign);

end

% display
% ---------------------------------------------------------------------
plot_design(3, 'online optimisation', u, y, posterior_online);

%% show results
% =========================================================================

fprintf('\nSimulation results:\n');

disp (table ( ...
    phi, ...
    posterior_naive.muPhi, ...
    posterior_offline.muPhi, ...
    posterior_online.muPhi, ...
    'RowNames', {'center','slope'}, ...
    'VariableNames',{'true','naive','offline','online'}));

%% ########################################################################
%  display subfunction
%  ########################################################################

function plot_design(idx, titleTxt, u, y, posterior, uEfficiency, dEfficiency)

    % jitter for data display
    persistent jitter;
    if isempty(jitter)
        jitter = 0.1 * (rand(N,1)-0.5);
    end

    % experimental design
    subplot(4,3,idx)

    if nargin > 5 && ~ isempty (uEfficiency)% if efficiency given

        [ax,h1,h2] = plotyy(u,10,uRange,uEfficiency,@myHistogram,@myPlot);
         xlim([uRange(1)-0.05, uRange(end)+0.05]);
         set(get(ax(1), 'YLabel'), 'String', 'freq. of presentation')
         set(get(ax(2), 'YLabel'), 'String', 'efficiency')
         set(ax(1),'YLim', [0 0.4],'YTick',0:.2:.4)
         xlabel('stimulus intensity')
         box off

    else

     % show stimuli density
     histogram(u, 10, 'EdgeColor','none','FaceColor',[.3 .3 .4],'Normalization','probability');
     xlim([uRange(1)-0.05, uRange(end)+0.05])
     ylim([0 0.4])
     xlabel('stimulus intensity')
     ylabel('freq. of presentation')
     box off
    end
     % show type of optimisation
     VBA_title(gca,titleTxt);

    if nargin > 6 % if efficiency given
        subplot(4,3,3+idx)
        plot(dEfficiency);
        xlim([1 numel(dEfficiency)])
        ylim([1.1*min(dEfficiency) 0])
        switch idx
            case 2
                xlabel('selected design');
            case 3
                xlabel('trial');
        end

        ylabel('efficiency');
        box off
    end

     % show results if any
     if ~isempty(posterior)

        % + observations

        subplot(4,3,6+idx)
        % predictions
        opt_plot = options;
        opt_plot.priors = posterior;
        dim_opt = dim;
        dim_opt.p = numel(uRange);
        muy = VBA_getLaplace(uRange',[],@g_psychometric,dim_opt,opt_plot);
        plot(uRange,muy,'r','LineWidth',2);
        % true model
        hold on
        plot(uRange,g_psychometric([],phi,uRange),'Color',[0 .8 0],'LineWidth',2);
        % data
        plot(u,y+jitter(1:numel(y)),'.k');

        % options
        ylim([-0.1 1.1])
        xlim([uRange(1)-0.05, uRange(end)+0.05])
        xlabel('stimulus intensity')
        ylabel('prob. of detection')
        hold off
        box off
        if idx == 1
            text(.2,.6,'true model','Color',[0 .8 0]);
            text(.2,.4,'observations','Color','k');
            text(.2,.2,'predicted','Color','r');
        end

        % + parameters

        subplot(4,3,9+idx);

        % posterior estimates
        plotUncertainTimeSeries(posterior.muPhi,sqrt(diag(posterior.SigmaPhi)),[],gca);
        % true values
        hold on
        plot(phi,'o','MarkerFaceColor',[0 .8 0], 'MarkerEdgeColor',[0 .8 0]);
        % options
        set(gca,'XTickLabel',{'center','slope'})
        ylim([-2 4])
        xlabel('parameter')
        ylabel('posterior estimate')
        hold off
        box off
    end

    drawnow
end

    function h = myHistogram (x,y)
        h = histogram(x,y,'EdgeColor','none','FaceColor',[.3 .3 .4],'Normalization','probability');
    end

    function h = myPlot (x,y)
        h = plot(x,y);
        hold on
        [mE, iE] = max (y);
        plot(uRange(iE),mE,'o','MarkerFaceColor',[.3 .3 .3],'MarkerEdgeColor',[.3 .3 .3])
        hold off
    end
end

author: MBB-team

repo: VBA-toolbox


text:
function Test_Mahal(file)
if nargin<1
    file='sample10k.csv';
end
UmapUtil.Initialize;
SAMPLE_FILE = 'sample10k.csv';
file = UmapUtil.RelocateExamples(file);
data= File.ReadCsv(file);

cov = nancov(data);
disp(['The covariance matrix for the data from ' SAMPLE_FILE ' is:']);
disp(cov);

invCov = inv(cov);
disp(['The inverse of the covariance matrix for the data from ' SAMPLE_FILE ' is:']);
disp(invCov);

X1 = knnsearch(data, data, 'K', 15, 'Distance', 'mahalanobis');
X2 = knnsearch(data, data, 'K', 15, 'Distance', 'mahalanobis', 'Cov', cov);

if isequal(X1, X2)
    disp('MATLAB''s knnsearch.m confirms that we correctly computed the covariance matrix!');
else
    disp('Hmm... did we calculate the covariance matrix incorrectly?');
end
X3 = KnnFind.Approximate(data, 15, 'mahalanobis', invCov, false, 3);
X4 = KnnFind.Approximate(data, 15, 'mahalanobis', [], false, 3);

acc = KnnFind.AssessApproximation(X1, X3);
disp(['nn_descent found nearest neighbors with ' num2str(100*acc) ' percent accuracy!']);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%{
Authors:
Jonatas Lopes de Paiva
Claudio Fabiano Motta Toledo
Helio Pedrini

%}

function pop = restartPop(pop, sizePop, denoisedImages, img, beta, lambda)

popTemp(sizePop).fitness = 0;
popTemp(sizePop).cromo = [];
popTemp(sizePop).id = '';
popTemp(1) = pop(1);

mut1 = denoisedImages(1).img;
mut2 = denoisedImages(2).img;
mut3 = denoisedImages(3).img;

fit1 = calcFitness(mut1, img, beta, lambda);
fit2 = calcFitness(mut2, img, beta, lambda);
fit3 = calcFitness(mut3, img, beta, lambda);

popTemp(2).fitness = fit1;
popTemp(2).cromo = mut1;
popTemp(2).id = char(java.util.UUID.randomUUID.toString);

popTemp(3).fitness = fit2;
popTemp(3).cromo = mut2;
popTemp(3).id = char(java.util.UUID.randomUUID.toString);

popTemp(4).fitness = fit3;
popTemp(4).cromo = mut3;
popTemp(4).id = char(java.util.UUID.randomUUID.toString);

for i = 5:sizePop
    tmp = denoisedImages(randi(3)).img;
    popTemp(i).cromo = mutation(tmp);
    popTemp(i).fitness = calcFitness(popTemp(i).cromo, img, beta, lambda);
    popTemp(i).id = char(java.util.UUID.randomUUID.toString);
end

pop = arrangePop(popTemp);

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function dy = ode_system(t, y, m, b, k)

%*****************************************************************************80
%
%% ODE_SYSTEM evaluates the right hand side of the ODE.
%
%  Discussion:
%
%    The second order ODE:
%
%      m * x'' + b * x' + k * x = 0
%
%    is transformed into a pair of first order ODE's
%    using the variables:
%
%
%      y(1) = x,
%      y(2) = x'
%
%    so that
%
%      y'(1) = y(2)
%      y'(2) = - ( k / m ) y(1) - ( b / m ) y(2)
%
%  Modified:
%
%    01 March 2010
%
%  Parameters:
%
%    Input, real T, the current time.
%
%    Input, real Y(2), the current solution.
%
%    Input, real M, B, K, the mass, damping and stiffness constants.
%
%    Output, real DY(2,1), the right hand sides of the ODE's,
%    returned as a column vector.
%
  dy(1,1) = y(2);
  dy(2,1) = - ( k / m ) * y(1) - ( b / m ) * y(2);

  return
end


text:
function [ value, ifault ] = digamma ( x )

%*****************************************************************************80
%
%% DIGAMMA calculates DIGAMMA ( X ) = d ( LOG ( GAMMA ( X ) ) ) / dX
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 June 2013
%
%  Author:
%
%    Original FORTRAN77 version by Jose Bernardo.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jose Bernardo,
%    Algorithm AS 103:
%    Psi ( Digamma ) Function,
%    Applied Statistics,
%    Volume 25, Number 3, 1976, pages 315-317.
%
%  Parameters:
%
%    Input, real X, the argument of the digamma function.
%    0 < X.
%
%    Output, real DIGAMMA, the value of the digamma function at X.
%
%    Output, integer IFAULT, error flag.
%    0, no error.
%    1, X <= 0.
%

%
%  Check the input.
%
  if ( x <= 0.0 )
    value = 0.0;
    ifault = 1;
    return
  end
%
%  Initialize.
%
  ifault = 0;
  value = 0.0;
%
%  Use approximation for small argument.
%
  if ( x <= 0.00001 )
    euler_mascheroni = 0.57721566490153286060;
    value = - euler_mascheroni - 1.0 / x;
    return
  end
%
%  Reduce to DIGAMA(X + N).
%
  while ( x < 8.5 )
    value = value - 1.0 / x;
    x = x + 1.0;
  end
%
%  Use Stirling's (actually de Moivre's) expansion.
%
  r = 1.0 / x;
  value = value + log ( x ) - 0.5 * r;
  r = r * r;
  value = value ...
    - r * ( 1.0 / 12.0 ...
    - r * ( 1.0 / 120.0 ...
    - r *   1.0 / 252.0 ) );

  return
end


text:
function shape = load_ply(filename)

[tri,vert] = ply_read(filename, 'tri');
shape = {};
shape.VERT = vert';
shape.TRIV = tri';

shape.n = size(shape.VERT,1);
shape.m = size(shape.TRIV,1);

end

author: OshriHalimi


text:
% Author: Guanghan Ning, Dec 2016
clear all;
addpath('../utils_eval_mpii/');
addpath('../eval_MPII/');

% Load annotations: validation
load('../eval_MPII/mpii_Tompson_validation', 'annolist');
load('../eval_MPII/mpii_Tompson_validation', 'RELEASE_img_index');
load('../eval_MPII/mpii_Tompson_validation', 'RELEASE_person_index');
annolist_validate = annolist;
num_validate_imgs = size(annolist_validate, 2);

% Load annotations: train + validation
load('../eval_MPII/mpii_human_pose_v1_u12_1','RELEASE');
annolist_all = RELEASE.annolist;

% Load scales extracted from image height
scale_in_cpp = load('../scale_in_cpp_MPII_validate_bigger_1.4.mat');
scale_in_cpp = scale_in_cpp.output_scale;

%%
img_folder_path = '/home/ngh/dev/POSE/POSE-dev/dataset/MPI/images/';%'../../../dataset/MPI/images/';

for scale_id = 1
    for iter_id = 1176058
        mat_folder_path = strcat('../dataset_mpii/results_all_validate_cropped_bigger_1.4/', num2str(iter_id), '_', num2str(scale_id-1), '/quantitative/');

        for img_id = 1:num_validate_imgs

            img_index = RELEASE_img_index(img_id);
            img_name = annolist_all(img_index).image.name;
            img_path = strcat(img_folder_path, '/', img_name);
            img_org_full = imread(img_path);

            % Display predicted results on the original image
            figure;
            imshow(img_org_full);
            hold on;

            % Get rect id. Note the difference between rect_ct and rect_id
            rect_id = RELEASE_person_index(img_id);

            % Find scale and mid-points in order to convert prediction
            x_mid = annolist_validate(img_id).annorect.objpos.x;
            y_mid = annolist_validate(img_id).annorect.objpos.y;
            mid_point = [x_mid, y_mid];  %scale = annolist_test(img_id).annorect.scale;

            % Load our algorithm predictions for this rect
            [folder_path, base_name, extension] = fileparts(img_name);
            [x_preds_raw, y_preds_raw] = load_pred_result(mat_folder_path, base_name, rect_id);

            for joint_id = 1:14
                % Convert predictions for un-cropped image scale for evaluation
                %id = convert_joint_order(joint_id);
                id = convert_joint_order_to_gt_order(joint_id);

                pred_point_raw = [x_preds_raw(joint_id), y_preds_raw(joint_id)];
                [x_pred, y_pred]= convert_trans_and_scale_validation(pred_point_raw, scale_in_cpp(img_id), mid_point);

                plot(x_pred, y_pred,'o');
            end
            hold off;
            pause;
        end
    end
end

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
function objectiveAbbr = printObjective(model)
% Prints out the Stoichiometric Coefficients for each
% Metabolite, with the name of the objective
%
% USAGE:
%
%    objectiveAbbr = printObjective(model)
%
% INPUT:
%    model:            COBRA model structure
%
% OUTPUT:
%    objectiveAbbr:    Objective reaction abbreviation
%
% .. Authors:
%       - Ronan Fleming 22/10/2008
%       - Thomas Pfau 15/12/2015 - Made the function compatible with sparse S matrices
%       - Laurent Heirendt March 2017 - Compatibility with large models and conversion to table

objRxnInd = find(model.c ~= 0);
objectiveAbbr = model.rxns(objRxnInd);
T = cell(length(objRxnInd), 1);

if isempty(objRxnInd)
    warning('There is no objective!');
else
    objMetVect = {};
    objRxnVect = {};
    objCoeffVect = {};
    for k = 1:length(objRxnInd)
        objMetInd = find(model.S(:, objRxnInd(k)));
        objMetVect{k} = model.mets(objMetInd);
        rxnName = model.rxns(objRxnInd(k));
        objCoeffVect{k} = full(model.S(objMetInd, objRxnInd(k)));

        objRxnVect = {};
        objRxnIDVect = {};
        % fill the list with the reaction name and ID
        for m = 1:length(objMetInd)
            objRxnVect{m} = char(rxnName);
            objRxnIDVect{m} = objRxnInd(k);
        end

        % save the table for reaction k
        T{k} = table(objCoeffVect{k}, categorical(objMetVect{k}), objMetInd, categorical(objRxnVect'), cell2mat(objRxnIDVect'), ...
                     'VariableNames', {'Coefficient', 'Metabolite', 'metID', 'Reaction', 'RxnID'});

        % concatenate the tables
        if k == 1
            summaryT = T{1};
        else
            summaryT = vertcat(summaryT, T{k});
        end
    end

    % display a summary
    display(summaryT);
end

author: opencobra

repo: cobratoolbox


text:
function [ labels, valid_ids, filenames  ] = extract_FERA2011_labels( FERA2011_dir, recs, aus )
%EXTRACT_SEMAINE_LABELS Summary of this function goes here
%   Detailed explanation goes here

    num_files = numel(recs);

    % speech invalidates lower face AUs
    labels = cell(num_files, 1);
    valid_ids = cell(num_files, 1);
    filenames = cell(num_files, 1);

    file_id = 1;

    for i=1:numel(recs)

        file = [FERA2011_dir, '/', recs{i}, '/', recs{i}, '-au.dat'];

        [~, filename,~] = fileparts(file);
        filenames{file_id} = filename;

        delimiter = {'   '};
        formatSpec = '%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%[^\n\r]';

        %% Open the text file.
        fileID = fopen(file,'r');

        %% Read columns of data according to the format.
        data = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string',  'ReturnOnError', false);
        data = [data{1:end-1}];

        %% Close the text file.
        fclose(fileID);

        %data = dlmread(file, ' '); %import annotations for one video file

        speech = data(:,end);

        labels{file_id} = data(:, aus);

        % Finding the invalid regions
        if(aus(1) >= 10)
            valid = speech == 0;
        else
            valid = true(size(speech,1), 1);
        end

        % all indices in SEMAINE are valid
        valid_ids{file_id} = valid;

        file_id = file_id + 1;
    end

    labels = labels(1:file_id-1);
    valid_ids = valid_ids(1:file_id-1);
    filenames = filenames(1:file_id-1);

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function dtiFiberGroupMniTransform(data_directory, fgfilename)

%Input is a file with a fiber group
%can be either .mat or .pdb file
%It is assumed that is is located in dtidir/fibers
%Output is the same fiber structure tranformed into MNI space.

%ER 04/2008
[pathstr, name, ext, versn] = fileparts(fgfilename);

    dt6File=fullfile(data_directory, 'dt6.mat');
    [sn, def]=dtiComputeDtToMNItransform(dt6File);

    filename=fullfile(data_directory, 'fibers', fgfilename);

    if strcmp(ext, '.pdb')
    [fg,filename] = mtrImportFibers(filename, eye(4));
    elseif strcmp(ext, '.mat')
        load(filename);
    else
display('Only .mat and .pdb types are supported for fgfilename');
    end

    fg = dtiXformFiberCoords(fg, def); %fg_sn = dtiXformFiberCoords(fg, def);

    coordinateSpace='MNI'; versionNum=1;
    newfilename=[filename(1:end-4) 'MNI'];
    save(newfilename, 'fg', 'coordinateSpace', 'versionNum');

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [x,e_conn,index_b] = twod_mesh(x_l, x_r, y_l, y_r, ...
                                        etype, ...
                                        n_nodesx, n_nodesy, ...
                                        fig_num )
%-----------------------------------------------------------------------
%  twod_mesh.m - Generate a rectangular mesh with a prescribed density.
%                This routine returns nodal coordinates, element
%                connectivity, and the nodal indices of boundary nodes.
%
%  Copyright (c) 2001, Jeff Borggaard, Virginia Tech
%  Version: 1.0
%
%  Usage:    [x,e_conn,index_b] = twod_mesh(x_l,x_r,y_l,y_r,etype,...
%                                           n_nodex,n_nodesy,fig_num)
%
%  Variables:     (x_l,y_l)
%                        coordinates of lower left corner
%                 (x_r,y_r)
%                        coordinates of upper right corner
%
%                 etype
%                        element type
%                          'linear', 'quadratic', 'cubic'
%
%                 n_nodesx, n_nodesy
%                        number of nodes in x and y directions
%                        (must be compatible with element type)
%
%                 fig_num  (optional)
%                        specifies the figure number for mesh plot
%    Outputs:
%                 x
%                        node coordinates of mesh
%                 e_conn
%                        element connectivity of mesh
%                 index_b
%                        node numbers of boundary nodes
%-----------------------------------------------------------------------

  if ( strcmp( etype, 'linear' ) )
    %  Generate node coordinates
    dx = (x_r-x_l)/(n_nodesx-1);
    dy = (y_r-y_l)/(n_nodesy-1);
    for i=1:n_nodesx
      for j=1:n_nodesy
        k = i+(j-1)*n_nodesx;
        x(k,1) = x_l + dx*(i-1);
        x(k,2) = y_l + dy*(j-1);
      end
    end

    %  Generate element connectivity
    e_conn = zeros(n_elements, 3);
    ie     = 0;
    for j=1:n_nodesy-1
      for i=1:n_nodesx-1
        ie = ie + 1;
        k = i + (j-1)*n_nodesx;
        e_conn(ie,:) = [k, k+1+n_nodesx, k+n_nodesx];
        ie = ie + 1;
        e_conn(ie,:) = [k, k+1, k+1+n_nodesx];
      end
    end

  elseif ( strcmp( etype, 'quadratic') )
    %  Generate a mesh with quadratic elements
    n_nodes  = n_nodesx*n_nodesy;
    n_elements = 2*(n_nodesx-1)*(n_nodesy-1)/4;

    dx = (x_r-x_l)/(n_nodesx-1);
    dy = (y_r-y_l)/(n_nodesy-1);
    for i=1:n_nodesx
      for j=1:n_nodesy
        k = i+(j-1)*n_nodesx;
        x(k,1) = x_l + dx*(i-1);
        x(k,2) = y_l + dy*(j-1);
      end
    end

    %  Generate element connectivity
    e_conn = zeros(n_elements, 6);
    ie     = 0;
    for j=1:2:n_nodesy-1
      for i=1:2:n_nodesx-1
        ie = ie + 1;
        k = i + (j-1)*n_nodesx;
        e_conn(ie,:) = [k, k+2+2*n_nodesx, k+2*n_nodesx, ...
                        k + 1 + n_nodesx, k+1+2*n_nodesx, k+n_nodesx ];
        ie = ie + 1;
        e_conn(ie,:) = [k, k+2, k+2+2*n_nodesx, ...
                        k+1, k+2+n_nodesx, k+1+n_nodesx ];
      end
    end

  elseif ( strcmp( etype, 'cubic' ) )
    error('twod_mesh: cubic elements are not currently implemented\n')
  else
    error('twod_mesh: etype is not a valid string\n')
  end

  if ( nargout==3 )
    index_b = [ 2:n_nodesx-1, ...
               1:n_nodesx:(n_nodesy-1)*n_nodesx, ...
               n_nodesx:n_nodesx:n_nodesy*n_nodesx, ...
               (n_nodesy-1)*n_nodesx+1:n_nodes-1 ];
  end

  if ( nargin==8 )
    twod_plotm2(fig_num,x,e_conn,'+')
  end

author: johannesgerer

repo: jburkardt-m