text:
% K-nearest Neighbor (9/12/2020)

function Acc = jknn(feat,label,opts)
% Default of k-value
k = 5;

if isfield(opts,'k'), k = opts.k; end
if isfield(opts,'Model'), Model = opts.Model; end

% Define training & validation sets
trainIdx = Model.training;    testIdx = Model.test;
xtrain   = feat(trainIdx,:);  ytrain  = label(trainIdx);
xvalid   = feat(testIdx,:);   yvalid  = label(testIdx);
% Training model
My_Model = fitcknn(xtrain,ytrain,'NumNeighbors',k);
% Prediction
pred     = predict(My_Model,xvalid);
% Accuracy
Acc      = sum(pred == yvalid) / length(yvalid);

fprintf('\n Accuracy: %g %%',100 * Acc);
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
function tf = ispolycw_N(x, y)
%ISPOLYCW_N True if polygon vertices are in clockwise order
%
%   TF = ISPOLYCW(X, Y) returns true if the polygonal contour vertices
%   represented by X and Y are ordered in the clockwise direction.  X and Y
%   are numeric vectors with the same number of elements.
%
%   Alternatively, X and Y can contain multiple contours, either in
%   NaN-separated vector form or in cell array form.  In that case,
%   ISPOLYCW returns a logical array containing one true or false value
%   per contour.
%
%   ISPOLYCW always returns true for polygonal contours containing two or
%   fewer vertices.
%
%   Vertex ordering is not well defined for self-intersecting polygonal
%   contours.  For such contours, ISPOLYCW returns a result based on the
%   order or vertices immediately before and after the left-most of the
%   lowest vertices.  In other words, of the vertices with the lowest Y
%   value, find the vertex with the lowest X value.  For a few special
%   cases of self-intersecting contours, the vertex ordering cannot be
%   determined using only the left-most of the lowest vertices; for these
%   cases, ISPOLYCW uses a signed area test to determine the ordering.
%
%   Class Support
%   -------------
%   X and Y may be any numeric class.
%
%   Example
%   -------
%   Orientation of a square:
%
%       x = [0 1 1 0 0];
%       y = [0 0 1 1 0];
%       ispolycw(x, y)                     % Returns 0
%       ispolycw(fliplr(x), fliplr(y))     % Returns 1
%
%   See also POLY2CW, POLY2CCW, POLYBOOL.

% Copyright 2004-2009 The MathWorks, Inc.
% $Revision: 1.1.4.3 $  $Date: 2009/08/11 15:44:28 $

if isempty(x)
   tf = true;
   return;
end

if ~iscell(x)
   checkxy_N(x, y, mfilename, 'X', 'Y', 1, 2)
   is_row = (size(x,1) == 1);
   [x, y] = polysplit_N(x, y);
   if is_row
      x = x';
      y = y';
   end
end

tf = false(size(x));
for k = 1:numel(x)
   tf(k) = isContourClockwise(x{k}, y{k});
end

%----------------------------------------------------------------------
function tf = isContourClockwise(x, y)

if numel(x) <= 1
    tf = true;
    return;
end

is_closed = (x(1) == x(end)) && (y(1) == y(end));
if is_closed
    x(end) = [];
    y(end) = [];
end

[x, y] = removeDuplicates(x, y);
num_vertices = numel(x);
if num_vertices <= 2
    tf = true;
    return;
end

idx = findExtremeVertices(x, y);

if numel(idx) > 1
    % The same extreme vertex appears multiple, nonsuccessive times in
    % the vertex list.  Use signed area test.
    tf = signedArea(x, y) <= 0;
    return;
end

% Find the three vertices we are interested in: the left-most of the
% lowest vertices, as well as the ones immediately before and after it.
p = mod((idx - 1) + [-1, 0, 1], num_vertices) + 1;
xx = x(p);
yy = y(p);

if ~isfloat(xx)
    xx = double(xx);
end
if ~isfloat(yy)
    yy = double(yy);
end

ux = xx(2) - xx(1);
uy = yy(2) - yy(1);

vx = xx(3) - xx(2);
vy = yy(3) - yy(2);

a = ux*vy;
b = uy*vx;
if a == b
    % The left-most lowest vertex is the end-point of a kind of linear
    % "spur."  The contour doubles back on itself, such as in this case:
    % x = [0 1 1 0 0 -1 0];
    % y = [0 0 1 1 0 -1 0];
    % The left-most lowest vertex is (-1,-1), but we since this vertex
    % is the end-point of a spur, we can't tell the direction from it.
    % Use the signed polygon test.
    tf = signedArea(x, y) <= 0;
else
    tf = a < b;
end

%----------------------------------------------------------------------
function [xout, yout] = removeDuplicates(x, y)
num_vertices = numel(x);
k1 = [2:num_vertices 1];
k2 = 1:num_vertices;
dups = (x(k1) == x(k2)) & (y(k1) == y(k2));
xout = x;
yout = y;
xout(dups) = [];
yout(dups) = [];

%----------------------------------------------------------------------
function idx = findExtremeVertices(x, y)
% Return the indices of all the left-most lowest vertices in (x,y).

% Find the vertices with the minimum y.
idx = find(y == min(y));

x_subset = x(idx);
idx2 = (x_subset == min(x_subset));

idx = idx(idx2);

%----------------------------------------------------------------------
function a = signedArea(x, y)
% a = signedArea(x,y) returns twice the signed area of the polygonal
% contour represented by vectors x and y.  Assumes (x,y) is NOT closed.

% Reference:
% http://geometryalgorithms.com/Archive/algorithm_0101/algorithm_0101.htm

x = x - mean(x);
n = numel(x);
if n <= 2
    a = 0;
else
    i = [2:n 1];
    j = [3:n 1 2];
    k = (1:n);
    a = sum(x(i) .* (y(j) - y(k)));
end

author: zouchuhang

repo: LayoutNet


text:
function [aw,ae,ap,su]=abc

global dx n2 ta tb
warning off
b=n2*(dx^2);

i=1;
aw(i)=0;
ae(i)=1;
su(i)=b*ta+2*tb;
sp(i)=-b-2;
ap(i)=aw(i)+ae(i)-sp(i);

i=2;
aw(i)=1;
ae(i)=1;
su(i)=b*ta;
sp(i)=-b;
ap(i)=aw(i)+ae(i)-sp(i);

i=3;
aw(i)=1;
ae(i)=0;
su(i)=b*ta;
sp(i)=-b;
ap(i)=aw(i)+ae(i)-sp(i);

author: Sable

repo: mcbench-benchmarks


text:
#!/usr/bin/env octave
%% Machine Learning Online Class
%  Exercise 7 | Principle Component Analysis and K-Means Clustering
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     pca.m
%     projectData.m
%     recoverData.m
%     computeCentroids.m
%     findClosestCentroids.m
%     kMeansInitCentroids.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% ================= Part 1: Find Closest Centroids ====================
%  To help you implement K-Means, we have divided the learning algorithm
%  into two functions -- findClosestCentroids and computeCentroids. In this
%  part, you shoudl complete the code in the findClosestCentroids function.
%
fprintf('Finding closest centroids.\n\n');

% Load an example dataset that we will be using
load('ex7data2.mat');

% Select an initial set of centroids
K = 3; % 3 Centroids
initial_centroids = [3 3; 6 2; 8 5];

% Find the closest centroids for the examples using the
% initial_centroids
idx = findClosestCentroids(X, initial_centroids);

fprintf('Closest centroids for the first 3 examples: \n')
fprintf(' %d', idx(1:3));
fprintf('\n(the closest centroids should be 1, 3, 2 respectively)\n');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ===================== Part 2: Compute Means =========================
%  After implementing the closest centroids function, you should now
%  complete the computeCentroids function.
%
fprintf('\nComputing centroids means.\n\n');

%  Compute means based on the closest centroids found in the previous part.
centroids = computeCentroids(X, idx, K);

fprintf('Centroids computed after initial finding of closest centroids: \n')
fprintf(' %f %f \n' , centroids');
fprintf('\n(the centroids should be\n');
fprintf('   [ 2.428301 3.157924 ]\n');
fprintf('   [ 5.813503 2.633656 ]\n');
fprintf('   [ 7.119387 3.616684 ]\n\n');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =================== Part 3: K-Means Clustering ======================
%  After you have completed the two functions computeCentroids and
%  findClosestCentroids, you have all the necessary pieces to run the
%  kMeans algorithm. In this part, you will run the K-Means algorithm on
%  the example dataset we have provided.
%
fprintf('\nRunning K-Means clustering on example dataset.\n\n');

% Load an example dataset
load('ex7data2.mat');

% Settings for running K-Means
K = 3;
max_iters = 10;

% For consistency, here we set centroids to specific values
% but in practice you want to generate them automatically, such as by
% settings them to be random examples (as can be seen in
% kMeansInitCentroids).
initial_centroids = [3 3; 6 2; 8 5];

% Run K-Means algorithm. The 'true' at the end tells our function to plot
% the progress of K-Means
[centroids, idx] = runkMeans(X, initial_centroids, max_iters, true);
fprintf('\nK-Means Done.\n\n');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ============= Part 4: K-Means Clustering on Pixels ===============
%  In this exercise, you will use K-Means to compress an image. To do this,
%  you will first run K-Means on the colors of the pixels in the image and
%  then you will map each pixel on to it's closest centroid.
%
%  You should now complete the code in kMeansInitCentroids.m
%

fprintf('\nRunning K-Means clustering on pixels from an image.\n\n');

%  Load an image of a bird
A = double(imread('bird_small.png'));

% If imread does not work for you, you can try instead
%   load ('bird_small.mat');

A = A / 255; % Divide by 255 so that all values are in the range 0 - 1

% Size of the image
img_size = size(A);

% Reshape the image into an Nx3 matrix where N = number of pixels.
% Each row will contain the Red, Green and Blue pixel values
% This gives us our dataset matrix X that we will use K-Means on.
X = reshape(A, img_size(1) * img_size(2), 3);

% Run your K-Means algorithm on this data
% You should try different values of K and max_iters here
K = 16;
max_iters = 10;

% When using K-Means, it is important the initialize the centroids
% randomly.
% You should complete the code in kMeansInitCentroids.m before proceeding
initial_centroids = kMeansInitCentroids(X, K);

% Run K-Means
[centroids, idx] = runkMeans(X, initial_centroids, max_iters);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================= Part 5: Image Compression ======================
%  In this part of the exercise, you will use the clusters of K-Means to
%  compress an image. To do this, we first find the closest clusters for
%  each example. After that, we

fprintf('\nApplying K-Means to compress an image.\n\n');

% Find closest cluster members
idx = findClosestCentroids(X, centroids);

% Essentially, now we have represented the image X as in terms of the
% indices in idx.

% We can now recover the image from the indices (idx) by mapping each pixel
% (specified by it's index in idx) to the centroid value
X_recovered = centroids(idx,:);

% Reshape the recovered image into proper dimensions
X_recovered = reshape(X_recovered, img_size(1), img_size(2), 3);

% Display the original image
subplot(1, 2, 1);
imagesc(A);
title('Original');

% Display compressed image side by side
subplot(1, 2, 2);
imagesc(X_recovered)
title(sprintf('Compressed, with %d colors.', K));

fprintf('Program paused. Press enter to continue.\n');
pause;

author: SaveTheRbtz

repo: ml-class


text:
function [ dbCA, cones, dir, Q, P, dualized ] = eliminate( prob, destructive, can_dual )
if nargin < 3, can_dual = nargout >= 6; end
if nargin < 2, destructive = false; end

% For the problem
%
%    minimize c' * x + d
%    s.t.     y : A * x + b == 0
%             x \in K
%
% The Lagrangian is
%
%  L(x,y,z) = c' * x + d - y' * ( A x + b ) - z' * x
%
%                            [ - d  b' 0 ]   [ 1 ]
%            = - [ 1, x' ] * [ - c  A' I ] * [ y ]
%                                            [ z ]
%
% This function provides a smaller [ d, b' ; c, A' ] with no more nonzeros
% that solves an equivalent problem. The original x and y can be recovered
% from the reduced xx and yy by Q*[1;xx] and P*[1;-yy], respectively.

[ dbCA, cones, dir, Q, P ] = extract( prob, destructive );
dualized = false;
if size( dbCA, 1 ) == 1,
    return;
end

%
% Negate the objective so that the transformation matrices P and Q are
% properly formed.
%

dbCA(:,1) = -dbCA(:,1);
if ~issparse( dbCA ),
    dbCA = sparse( dbCA );
end

for pass = 1 : 2,

    if pass == 1 || dualized,
        n_tot = 0;
        nn  = size(dbCA,1);
        rsv = sparse( 1, 1, 1, nn, 1 );
        nng = sparse( nn, 1 );
        for k = 1 : length( cones ),
            temp = cones(k).indices;
            n_tot = n_tot + numel(temp);
            temp = sparse( temp, 1, 1, nn, 1 );
            rsv = rsv + temp;
            if isequal( cones(k).type, 'nonnegative' ),
                nng = nng + temp;
            elseif can_dual && strncmp( cones(k).type, 'i_', 2 ),
                can_dual = false;
            end
        end
        rsv   = full( rsv );
        nng   = full( nng );
        ndxs  = ( 1 : nn )';
        nold  = nn;
    end

    cc   = dbCA( :, 1 );
    rcnt = sum( dbCA ~= 0, 2 );

    % In the first pass, we don't eliminate columns which have inequality
    % structure to them, so that we can make the best decision as to
    % whether or not to convert the problem to dual standard form. Exempted
    % from this are columns with trivial xi = xj constraints, where xi is free.
    if pass == 1,
        trivs = sum( dbCA(rsv==0,:) ~= 0, 1 ) == 1 & sum( dbCA(rsv~=0,:) ~= 0, 1 ) - ( dbCA( 1, : ) ~= 0 ) == 1;
        ineqs = full(any(dbCA(rsv~=0&rcnt==1,:),1)) & full(~trivs);
        ineqs = +ineqs;
    else % if dualized,
        ineqs = zeros(1,size(dbCA,2));
    end
    ineqs(1) = 1;

    while true,

        success = false;

        %
        % STEP 1: Look for free or nonnegative variables that do not appear
        % in any constraints. Unconstrained variables that also appear in
        % the objective are unbounded, as are nonnegative variables that
        % appear there with the right sign. Otherwise their values might
        % as well be zero. If we have multiple unbounded variables, keep
        % all but one so that the solver can still see this happen.
        %
        % Eliminated for now. Frankly, my suspicion is that this happens
        % very infrequently, and this code seems to have been the source
        % of bugs in the past.
        %

        if 0,
            rows = ( rcnt == ( cc ~= 0 ) ) & ( ~rsv | nng );
            nnzr = nnz( rows );
            if nnzr > 0,
                csgn = 1 - 2 * dualized;
                celm = csgn * cc( rows, 1 );
                celm( nng(rows) & celm < 0 ) = 0;
                nnzc = nnz( celm );
                if nnzc > 1 || nnzr > nnzc,
                    success = true;
                    if nnzc,
                        cnrm = norm( celm );
                        ndxq = find( rows );
                        ndxq = ndxq( celm ~= 0 );
                        ndxq = ndxq( 1 );
                        Q( :, ndxq ) = Q( :, rows ) * ( celm / cnrm );
                        dbCA( ndxq, 1 ) = csgn * cnrm; %#ok
                        rows( ndxq ) = 0;
                    end
                    rowX = ~rows;
                    dbCA = dbCA( rowX, : );
                    rsv  = rsv ( rowX, : );
                    nng  = nng ( rowX, : );
                    ndxs = ndxs( rowX, : );
                    Q    =    Q( :, rowX );
                end
            end
        end

        %
        % STEP 2: Look for columns which differ only by a constant factor.
        % These correspond to redundant equality constraints. These occur
        % often enough as as consequence of our tranformation method, and
        % they cause problems in solvers, so we must eliminate them. Of
        % course, if there are more complex linear dependencies in the
        % equality constraints, we can't do anything about that.
        %

        [ xR, dbCA ] = cvx_bcompress( dbCA, 'full', 1 );
        if size( xR, 1 ) ~= size( xR, 2 ),
            success = true;
            P       = P * cvx_invert_structure( xR );
            ineqs   = ( xR * ineqs(:) )' ~= 0;
            ineqs   = +ineqs;
        end

        while true,

            %
            % STEP 3: Look for variables that we can eliminate without
            % increasing fill-in. This means looking for rows or columns
            % with only 1, 2, or (in some cases) 3 nonzeros.
            %

            [ rows, cols ] = cvx_eliminate_mex( dbCA, 1, rsv, ineqs );
            if ~any( rows ), break; end
            success = true;
            rows = rows ~= 0;
            cols = cols ~= 0;
            rowX = ~rows;
            colX = ~cols;

            %
            % [ x1^T x2^T ] [ C1 A11 A12 ] [ 1  ]
            %               [ C2 A21 A22 ] [ y1 ] = 0
            %                              [ y2 ]
            %
            % [ x1^T x2^T ] = x1^T [ I -A12*A22i ]
            %
            % [ G Y1^T Y2^T ] = [ G Y1^T ] [ I  0  -C2'*A22i'  ]
            %                              [ 0  I  -A21'*A22i' ]

            %
            A11  = dbCA( rowX, colX );
            A12  = dbCA( rowX, cols );
            A21  = dbCA( rows, colX );
            A22  = dbCA( rows, cols );
            if ( size( A22, 1 ) ~= size( A22, 2 ) || nnz( A22 ) ~= size( A22, 1 ) ),
                error( 'There seems to be an error in the CVX presolver routine.\nPlease report this to the authors; and if possible, include the\ncvx model and data that gave you this error.', 1 ); %#ok
            end
            [ ii, jj, vv ] = find( A22 );
            A22i  = sparse( jj, ii, 1.0 ./ vv );
            temp  = - A22i * A21;
            P     = P( :, colX ) + P( :, cols ) * temp;
            temp  = - A12 * A22i;
            Q     = Q( :, rowX ) + Q( :, rows ) * temp';
            dbCA  = A11 + temp * A21;
            rsv   =   rsv( rowX, : );
            nng   =   nng( rowX, : );
            ndxs  =  ndxs( rowX, : );
            ineqs = ineqs( :, colX );

        end

        if ~success,
            break;
        end

        cc   = dbCA( :, 1 );
        rcnt = sum( dbCA ~= 0, 2 );

    end

    if pass == 2 || isempty(cones) || ~can_dual,
        break;
    end

    %
    % Check to see if dualization will result in smaller problem
    %
    ineqs(1) = 0; rsv(1) = 0;
    n_save = nnz(cvx_eliminate_mex(dbCA,1,rsv,zeros(size(ineqs))));
    % n_save = nnz(sum(dbCA(:,ineqs~=0)~=0,1)==1+(dbCA(1,ineqs~=0)~=0));
    n_ineq = nnz(any(dbCA(rsv&rcnt==(cc~=0)+1,:)));
    rsv(1) = 1; ineqs(1) = 1; %#ok
    [n1,m1] = size(dbCA);
    m_pri = m1 - n_save - 1;
    n_pri = n1 - n_save - 1;
    m_dua = n1 - n_ineq - 1;
    n_dua = nnz( rsv ) + m1 - n_ineq - 1;
    if ( ( m_pri > n_pri ) || ( m_pri * m_pri * n_pri > m_dua * m_dua * n_dua ) ) && ( m_dua <= n_dua ),
        ndxs = full(sparse(ndxs,1,1:n1));
        PP = cell(2,length(cones));
        n_cur = m1;
        for k = 1 : length(cones),
            temp = cones(k).indices;
            [nn,nv] = size(temp);
            temp = reshape(ndxs(temp),size(temp));
            switch cones(k).type,
                case 'semidefinite',
                    nt = 0.5*(sqrt(8*nn+1)-1);
                    SS = 'symmetric';
                case 'hermitian-semidefinite',
                    nt = sqrt(nn);
                    SS = 'hermitian';
                case 'exponential',
                    SS = sparse(inv([0,-1,0;-1,0,0;0,0,exp(1)]));
                    SS = cvx_replicate_structure(SS,nv);
                otherwise,
                    SS = [];
            end
            PP{k} = sparse(1:numel(temp),max(temp,1),temp~=0,numel(temp),n1);
            if ~isempty(SS),
                if ischar(SS),
                    SS = cvx_create_structure([nt,nt,nv],SS);
                    SS = SS * SS';
                end
                PP{k} = SS * PP{k};
            end
            cones(k).indices = reshape(n_cur+1:n_cur+nn*nv,nn,nv);
            n_cur = cones(k).indices(end);
        end
        dbCA  = vertcat(dbCA',PP{:});
        dir   = -dir;
        tmp   = Q; Q = P; P = tmp;
        nold  = size(dbCA,1);
        Q(:,nold) = 0;
        dualized = true;
    end

end

%
% Return the objective back to normal.
%

if dualized,
    P = -P;
    P(:,1) = -P(:,1);
else
    dbCA(:,1) = -dbCA(:,1);
end

%
% Move the cone indices to their new locations
%

ndxs = full( sparse( ndxs, 1, 1 : length( ndxs ), nold, 1 ) );
tt = zeros(1,length(cones));
for k = 1 : length( cones ),
    temp = ndxs(cones(k).indices);
    if all(temp),
        temp = reshape( temp, size(cones(k).indices) );
    else
        temp = nonzeros(temp);
        temp = reshape( temp, 1, length(temp) );
    end
    tt(k) = isempty(temp);
    cones(k).indices = temp;
end
if any(tt),
    cones(tt~=0) = [];
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% Assume a target positioned at x = 1, travelling with speed v = 0.1
state = [1;0.1;2;0];    % Assume state with four dimensions [x_pos, x_vel, y_pos, y_vel]
measurement = [1.2;1.8;3];% Measurement state [x_pos,y_pos]

% Create an instance of a 1D Constant Velocity model
obs = LinearGaussianX('NumMeasDims', 3,'NumStateDims',4,...
                      'MeasurementErrVariance',50,'Mapping',[1 3 2]);

% View the transition matrix and process covariance matrices
H = obs.feval();
R = obs.covar();

% Predict the target's position and velocity after the interval has passed
projectedState = obs.feval(state);

% Do the same as above, but this time add process noise to the prediction
projectedState2 = obs.feval(state,true);

% Generate 50 random noise samples from the dynamic model
noise = obs.random(50);

% Check how likely the predictions we made are
lik = obs.pdf(projectedState,state);
lik2 = obs.pdf(projectedState2,state); % HINT: newState2 should be less likely

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function b = r83_np_ml ( n, a_lu, x, job )

%*****************************************************************************80
%
%% R83_NP_ML computes A * x or x * A, where A has been factored by R83_NP_FA.
%
%  Discussion:
%
%    The R83 storage format is used for a tridiagonal matrix.
%    The superdiagonal is stored in entries (1,2:N), the diagonal in
%    entries (2,1:N), and the subdiagonal in (3,1:N-1).  Thus, the
%    original matrix is "collapsed" vertically into the array.
%
%  Example:
%
%    Here is how a R83 matrix of order 5 would be stored:
%
%       *  A12 A23 A34 A45
%      A11 A22 A33 A44 A55
%      A21 A32 A43 A54  *
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be at least 2.
%
%    Input, real A_LU(3,N), the LU factors from R83_FA.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Input, integer JOB, specifies the product to find.
%    0, compute A * x.
%    nonzero, compute A' * x.
%
%    Output, real B(N), the product A*x or A'*x.
%
  b(1:n) = x(1:n);

  if ( job == 0 )
%
%  Compute X := U * X
%
    for i = 1 : n

      b(i) = a_lu(2,i) * b(i);

      if ( i < n )
        b(i) = b(i) + a_lu(1,i+1) * b(i+1);
      end

    end
%
%  Compute X: = L * X.
%
    for i = n : -1 : 2
      b(i) = b(i) + a_lu(3,i-1) * b(i-1);
    end

  else
%
%  Compute X: = L' * X.
%
    for i = 1 : n-1
      b(i) = b(i) + a_lu(3,i) * b(i+1);
    end
%
%  Compute X: = U' * X.
%
    for i = n : -1 : 2
      b(i) = a_lu(2,i) * b(i);
      b(i) = b(i) + a_lu(1,i) * b(i-1);
    end
    b(1) = a_lu(2,1) * b(1);

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef MimCombinedImageDatabase < CoreBaseClass
    % MimCombinedImageDatabase. Part of the internal framework of the TD MIM Toolkit.
    %
    %     You should not use this class within your own code. It is intended to
    %     be used internally within the framework of the TD MIM Toolkit.
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = private)
        ImageDatabase
        MatNatDatabase
    end

    events
        DatabaseHasChanged
    end

    methods
        function obj = MimCombinedImageDatabase(image_database, matnat_database)
            obj.ImageDatabase = image_database;
            obj.MatNatDatabase = matnat_database;

            obj.AddEventListener(image_database, 'DatabaseHasChanged', @obj.ImageDatabaseChanged);
        end

        function [project_names, project_ids] = GetListOfProjects(obj)
            [project_names_1, project_ids_1] = obj.ImageDatabase.GetListOfProjects;
            [project_names_2, project_ids_2] = obj.MatNatDatabase.GetListOfProjects;
            project_names = [project_names_1, project_names_2];
            project_ids = [project_ids_1, project_ids_2];
        end

        function datasets = GetAllSeriesForThisPatient(obj, project_id, patient_id, group_patients)
            datasets_1 = obj.ImageDatabase.GetAllSeriesForThisPatient(project_id, patient_id, group_patients);
            datasets_2 = obj.MatNatDatabase.GetAllSeriesForThisPatient(project_id, patient_id);
            datasets = horzcat(datasets_1, datasets_2);
        end

        function [names, ids, short_visible_names, patient_id_map] = GetListOfPatientNames(obj, project_id, group_patients)

            [names1, ids1, short_visible_names1, patient_id_map1] = obj.ImageDatabase.GetListOfPatientNames(project_id, group_patients);
            [names2, ids2, short_visible_names2, patient_id_map2] = obj.MatNatDatabase.GetListOfPatientNames(project_id);

            names = horzcat(names1, names2);
            ids = horzcat(ids1, ids2);
            short_visible_names = horzcat(short_visible_names1, short_visible_names2);

            patient_id_map = patient_id_map1;

            for key = patient_id_map2.keys
                patient_id_map(key{1}) = patient_id_map2(key{1});
            end
        end

        function [names, ids, short_visible_names, num_series, num_patients_combined, patient_id_map] = GetListOfPatientNamesAndSeriesCount(obj, project_id, group_patients)

            [names1, ids1, short_visible_names1, num_series_1, num_patients_combined_1, patient_id_map1] = obj.ImageDatabase.GetListOfPatientNamesAndSeriesCount(project_id, group_patients);
            [names2, ids2, short_visible_names2, num_series_2, num_patients_combined_2, patient_id_map2] = obj.MatNatDatabase.GetListOfPatientNamesAndSeriesCount(project_id);

            names = horzcat(names1, names2);
            ids = horzcat(ids1, ids2);
            short_visible_names = horzcat(short_visible_names1, short_visible_names2);
            num_series = num_series_1 + num_series_2;
            num_patients_combined = num_patients_combined_1 + num_patients_combined_2;

            patient_id_map = patient_id_map1;

            for key = patient_id_map2.keys
                patient_id_map(key) = patient_id_map2(key);
            end
        end
    end

    methods (Access = private)
        function ImageDatabaseChanged(obj, ~, ~)
            notify(obj, 'DatabaseHasChanged');
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function set_font_type(font)
%SET_FONT_TYPE sets the figures font type
%
%   Usage: set_font_type(font)
%
%   Input parameters:
%       font   - font type name
%
%   See also: set_font_size, print_png

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);
isargchar(font)

%% ===== Apply settings ==================================================
% Get handle for active figure
h=gca;
% Set the font type for the axis tics
set(h,'FontName',font);
% set the font type for all text strings
set(findall(h,'type','text'),'FontName',font);
% FIXME: the following produces an segmentation fault under Octave
% Set the font type for colorbar tics
%h=colorbar;
%set(h,'FontName',font);

author: sfstoolbox

repo: sfs-matlab


text:
% Standard test images
classdef standard_images

methods(Static)

    function image_dir = test_image_dir()
        here = which('spx.data.standard_images');
        data = fileparts(here);
        spx = fileparts(data);
        library = fileparts(spx);
        spx = fileparts(library);
        image_dir = fullfile(spx, 'data', 'images');
    end

    function image = barbara_gray_512x512()
        % Returns the Barbara image
        filepath = fullfile(spx.data.standard_images.test_image_dir(), 'barbara_gray_512x512.png');
        image = spx.image.read_gray_image(filepath);
    end

    function image = cameraman_gray_512x512()
        % Returns the Barbara image
        filepath = fullfile(spx.data.standard_images.test_image_dir(), 'cameraman_gray_512x512.tif');
        image = spx.image.read_gray_image(filepath);
    end

    function image = jetplane_gray_512x512()
        % Returns the Barbara image
        filepath = fullfile(spx.data.standard_images.test_image_dir(), 'jetplane_gray_512x512.tif');
        image = spx.image.read_gray_image(filepath);
    end

    function image = lake_gray_512x512()
        % Returns the Barbara image
        filepath = fullfile(spx.data.standard_images.test_image_dir(), 'lake_gray_512x512.tif');
        image = spx.image.read_gray_image(filepath);
    end

    function image = lena_gray_512x512()
        % Returns the Barbara image
        filepath = fullfile(spx.data.standard_images.test_image_dir(), 'lena_gray_512x512.tif');
        image = spx.image.read_gray_image(filepath);
    end

    function image = peppers_gray_512x512()
        % Returns the Barbara image
        filepath = fullfile(spx.data.standard_images.test_image_dir(), 'peppers_gray_512x512.tif');
        image = spx.image.read_gray_image(filepath);
    end

    function image = pirate_gray_512x512()
        % Returns the Barbara image
        filepath = fullfile(spx.data.standard_images.test_image_dir(), 'pirate_gray_512x512.tif');
        image = spx.image.read_gray_image(filepath);
    end

end

end

author: indigits

repo: sparse-plex


text:
% SFR_Bif_Diag_Second_Method - Bifurcation diagram with history.
% Copyright Springer 2013 Stephen Lynch.
clear
format long;
halfN=19999;N=2*halfN+1;N1=1+halfN;
E(1)=0.4;B=0.15;Pmax=16;

% Ramp the power up
for n=1:halfN
    E(n+1)=sqrt(n*Pmax/N1)+B*E(n)*exp(1i*abs(E(n))^2);
    Esqr(n+1)=abs(E(n+1))^2;
end

% Ramp the power down
for n=N1:N
    E(n+1)=sqrt(2*Pmax-n*Pmax/N1)+B*E(n)*exp(1i*abs(E(n))^2);
    Esqr(n+1)=abs(E(n+1))^2;
end

for n=1:halfN
    Esqr1(n)=Esqr(N+1-n);
    ptsup(n)=n*Pmax/N1;
end

% Plot the bifurcation diagrams
fsize=15;
hold on
set(gca,'xtick',0:4:16,'FontSize',fsize)
set(gca,'ytick',0:5:25,'FontSize',fsize)
plot(ptsup(1:halfN),Esqr(1:halfN),'.','MarkerSize',1);
plot(ptsup(1:halfN),Esqr1(1:halfN),'.','MarkerSize',1);
xlabel('Input Power','FontSize',fsize);
ylabel('Output Power','FontSize',fsize);
hold off


text:
function dynamics = dynamicsTest(type);

% DYNAMICSTEST Run some tests on the specified dynamics model.

% FGPLVM

type = [type 'Dynamics'];
X = randn(10, 2);
model = modelCreate(type, 2, 2, X);

% Set the parameters randomly.
params = modelExtractParam(model);
params = randn(size(params))./sqrt(randn(size(params)).^2);
model = modelExpandParam(model, params);

epsilon = 1e-6;
params = modelExtractParam(model);
if ~isempty(params)
  origParams = params;
  for i = 1:length(params);
    params = origParams;
    params(i) = origParams(i) + epsilon;
    model = modelExpandParam(model, params);
    Lplus(i) = modelLogLikelihood(model);
    params(i) = origParams(i) - epsilon;
    model = modelExpandParam(model, params);
    Lminus(i) = modelLogLikelihood(model);
  end
  params = origParams;
  model = modelExpandParam(model, params);

  gLDiff = .5*(Lplus - Lminus)/epsilon;
  g = modelGradientParam(model);
  nameBase = 'Parameter ';
  for i = 1:length(params)
    names{i} = [nameBase num2str(i)];
  end
  paramMaxDiff = max(max(abs(gLDiff-g)));
  if paramMaxDiff > 2*epsilon
    l = 0;
    for i = 1:length(names)
      if l < length(names{i})
        l = length(names{i});
      end
    end

    fprintf([char(repmat(32, 1, l)) '\tanalytic   diffs     delta\n']);
    for i = 1:length(names)
      spaceLen = l - length(names{i});
      space = char(repmat(32, 1, spaceLen));
      fprintf([space names{i} ':\t%4.6f\t%4.6f\t%4.6f\n'], ...
              g(i), gLDiff(i), gLDiff(i) - g(i));
    end
  end
end
Lplus = zeros(size(X));
Lminus = zeros(size(X));
origX = X;
epsilon = 0.001;
for i = 1:size(X, 1)
  for j = 1:size(X, 2)
    X = origX;
    X(i, j) = origX(i, j) + epsilon;
    model = modelSetLatentValues(model, X);
    params = modelExtractParam(model);
    model = modelExpandParam(model, params);
    Lplus(i, j) = modelLogLikelihood(model);
    X(i, j) = origX(i, j) - epsilon;
    model = modelSetLatentValues(model, X);
    params = modelExtractParam(model);
    model = modelExpandParam(model, params);
    Lminus(i, j) = modelLogLikelihood(model);
  end
end
X = origX;
model = modelSetLatentValues(model, X);
gXDiff = .5*(Lplus - Lminus)/epsilon;
g = modelLatentGradients(model);
XMaxDiff = max(max(abs(g-gXDiff)));
fprintf('X max diff: %2.6f.\n', XMaxDiff)
if XMaxDiff > 1e-6
  disp(g)
  disp(gXDiff)
  disp(g - gXDiff);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% RES = reconSFpyr(PYR, INDICES, LEVS, BANDS, TWIDTH)
%
% Reconstruct image from its steerable pyramid representation, in the Fourier
% domain, as created by buildSFpyr.
%
% PYR is a vector containing the N pyramid subbands, ordered from fine
% to coarse.  INDICES is an Nx2 matrix containing the sizes of
% each subband.  This is compatible with the MatLab Wavelet toolbox.
%
% LEVS (optional) should be a list of levels to include, or the string
% 'all' (default).  0 corresonds to the residual highpass subband.
% 1 corresponds to the finest oriented scale.  The lowpass band
% corresponds to number spyrHt(INDICES)+1.
%
% BANDS (optional) should be a list of bands to include, or the string
% 'all' (default).  1 = vertical, rest proceeding anti-clockwise.
%
% TWIDTH is the width of the transition region of the radial lowpass
% function, in octaves (default = 1, which gives a raised cosine for
% the bandpass filters).

% Eero Simoncelli, 5/97.

function res = reconSFpyr(pyr, pind, levs, bands, twidth)

%%------------------------------------------------------------
%% DEFAULTS:

if (exist('levs') ~= 1)
  levs = 'all';
end

if (exist('bands') ~= 1)
  bands = 'all';
end

if (exist('twidth') ~= 1)
  twidth = 1;
elseif (twidth <= 0)
  fprintf(1,'Warning: TWIDTH must be positive.  Setting to 1.\n');
  twidth = 1;
end

%%------------------------------------------------------------

nbands = spyrNumBands(pind);

maxLev =  1+spyrHt(pind);
if strcmp(levs,'all')
  levs = [0:maxLev]';
else
  if (any(levs > maxLev) | any(levs < 0))
    error(sprintf('Level numbers must be in the range [0, %d].', maxLev));
  end
  levs = levs(:);
end

if strcmp(bands,'all')
  bands = [1:nbands]';
else
  if (any(bands < 1) | any(bands > nbands))
    error(sprintf('Band numbers must be in the range [1,3].', nbands));
  end
  bands = bands(:);
end

%----------------------------------------------------------------------

dims = pind(1,:);
ctr = ceil((dims+0.5)/2);

[xramp,yramp] = meshgrid( ([1:dims(2)]-ctr(2))./(dims(2)/2), ...
    ([1:dims(1)]-ctr(1))./(dims(1)/2) );
angle = atan2(yramp,xramp);
log_rad = sqrt(xramp.^2 + yramp.^2);
log_rad(ctr(1),ctr(2)) =  log_rad(ctr(1),ctr(2)-1);
log_rad  = log2(log_rad);

%% Radial transition function (a raised cosine in log-frequency):
[Xrcos,Yrcos] = rcosFn(twidth,(-twidth/2),[0 1]);
Yrcos = sqrt(Yrcos);
YIrcos = sqrt(abs(1.0 - Yrcos.^2));

if (size(pind,1) == 2)
  if (any(levs==1))
    resdft = fftshift(fft2(pyrBand(pyr,pind,2)));
  else
    resdft = zeros(pind(2,:));
  end
else
  resdft = reconSFpyrLevs(pyr(1+prod(pind(1,:)):size(pyr,1)), ...
      pind(2:size(pind,1),:), ...
      log_rad, Xrcos, Yrcos, angle, nbands, levs, bands);
end

lo0mask = pointOp(log_rad, YIrcos, Xrcos(1), Xrcos(2)-Xrcos(1), 0);
resdft = resdft .* lo0mask;

%% residual highpass subband
if any(levs == 0)
  hi0mask = pointOp(log_rad, Yrcos, Xrcos(1), Xrcos(2)-Xrcos(1), 0);
  hidft = fftshift(fft2(subMtx(pyr, pind(1,:))));
  resdft = resdft + hidft .* hi0mask;
end

res = real(ifft2(ifftshift(resdft)));

author: ricedsp


text:
function vbmc_plot(vp_array,stats)

if nargin < 2; stats = []; end

Nsamples = 1e5;

if ~iscell(vp_array)
    temp{1} = vp_array;
    vp_array = temp;
end

if numel(vp_array) == 1 && vbmc_isavp(vp_array{1})
    X = vbmc_rnd(vp_array{1},Nsamples);
    for d = 1:size(X,2); names{d} = ['x_{' num2str(d) '}']; end
    cornerplot(X,names);
else
    Nbins = 40;
    Nvps = numel(vp_array);
    D = vp_array{1}.D;
    mm = zeros(Nvps,D);
    cmap = colormap;
    cmap = cmap(mod((1:27:(1+27*64))-1,64)+1,:);

    plotmat = [1 1; 1 2; 1 3; 2 2; 2 3; 2 3; 2 4; 2 4; 3 3; 3 4; 3 4; 3 4; 3 5; 3 5; 3 5; 4 4; 4 5; 4 5; 4 5];
    nrows = plotmat(D,1);
    ncols = plotmat(D,2);

    for i = 1:Nvps
        if ~isempty(stats) && stats.idx_best == i; best_flag = true; else; best_flag = false; end
        ltext{i} = ['vp #' num2str(i)];
        if best_flag; ltext{i} = [ltext{i} ' (best)']; end

        X = vbmc_rnd(vp_array{i},Nsamples);
        mm(i,:) = median(X);

        for d = 1:D
            subplot(nrows,ncols,d);
            if best_flag; lw = 3; else; lw = 1; end
            hst(i)=histogram(X(:,d),Nbins,'Normalization','probability','Displaystyle','stairs','LineWidth',lw,'EdgeColor',cmap(i,:));
            hold on;
        end
    end

    for i = 1:Nvps
        if ~isempty(stats) && stats.idx_best == i; best_flag = true; else; best_flag = false; end
        for d = 1:D
            subplot(nrows,ncols,d);
            if best_flag; lw = 3; else; lw = 1; end
            hln(i)=plot(mm(i,d)*[1 1],ylim,'-','LineWidth',lw,'Color',cmap(i,:));
            hold on;
        end
    end

    for d = 1:D
        subplot(nrows,ncols,d);

        xlabel(['x_{' num2str(d) '}']);
        set(gca,'TickDir','out');
        box off;

        if d == D
            hleg = legend(hln,ltext{:});
            set(hleg,'box','off','location','best');
        end

    end
    set(gcf,'Color','w');

end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function [A,B]= useconvhulln(Z2)

% This file is part of SOSTOOLS - Sum of Squares Toolbox ver 3.00.
%
% Copyright (C)2002, 2004, 2013  A. Papachristodoulou (1), J. Anderson (1),
%                                G. Valmorbida (1), S. Prajna (2),
%                                P. Seiler (3), P. A. Parrilo (4)
% (1) Department of Engineering Science, University of Oxford, Oxford, U.K.
% (2) Control and Dynamical Systems - California Institute of Technology,
%     Pasadena, CA 91125, USA.
% (3) Aerospace and Engineering Mechanics Department, University of
%     Minnesota, Minneapolis, MN 55455-0153, USA.
% (4) Laboratory for Information and Decision Systems, M.I.T.,
%     Massachusetts, MA 02139-4307
%
% Send bug reports and feedback to: sostools@cds.caltech.edu
%
% 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 2 of the License, or
% (at your option) any later version.
%
% This program 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 this program. If not, see <http://www.gnu.org/licenses/>.

% AP Apr 03 02
% AP Feb 01 03

% First, find the convex hull of Z2

convh = convhulln(Z2);
facets = size(convh,1);
[nZ1,nZ2] = size(Z2);

% Form the hyperplanes
% A hyperplane of the form [x1 ... xn -ones][A';b'] = 0
% [A';b'] is in the nullspace of the matrix [x1,...,xn -ones].
% Ignore it if it has dimension greater than 1
% Normalise if possible.

for i = 1:facets
    vert = [Z2(convh(i,:),:) -ones(nZ2,1)];
    nullvert = null(vert,'r');
    if size(nullvert,2) == 1
        if nullvert(end) ~=0
            nullvert = nullvert./nullvert(end);
        end
        coeff(:,i) = nullvert(1:end-1);
        cold(i,1) = nullvert(end);
    end
end

coeff = coeff';

% Condition the matrix a bit

coeffcold = [coeff cold];
coeffcold = round(coeffcold*100000)/100000;

% Discard same hyperplanes (due to convhulln result)
[coeff2cold,Ix,Iy] = unique(coeffcold,'rows');

%Remove a possible zero row

yind = find(sum(coeff2cold.^2,2) == 0);
if ~isempty(yind)
    Ix = [Ix([1:yind-1],1);Ix(yind+1:end,1)];
end

coeff2 = coeff(Ix,:);
convhnew = convh(Ix,:);
cnew = cold(Ix);
facetsnew = size(convhnew,1);

% Make inequalities out of them by testing a point not on the hyperplane
% Notation: convex hull is now Ax-b<=0
for fac = 1:facetsnew
    for ind = 1:nZ1
        matr = find(convhnew(fac,:) - ind*ones(1,nZ2) == 0);
        tests(fac) = coeff2(fac,:)*Z2(ind,:)'-cnew(fac);
        if isempty(matr) & abs(tests(fac)) > 1e-8
            break
        end
    end
    if tests(fac)>0
        coeff2(fac,:)=-coeff2(fac,:);
        cnew(fac)=-cnew(fac);
    end
end

A = coeff2;
B = cnew;

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function h=m_patch(long,lat,C,varargin)
% M_PATCH Create patches on a map
%    M_PATCH(LONG,LAT,C) is a drop-in replacement for PATCH that uses
%    longitude/latitude coordinates to draw a patch on the current map.
%    See PATCH for more details about the way in which patch colours and
%    properties should be specified.
%
%    Currently you cannot specify C to be other than a string or 1x3 RGB
%    vector.
%
%    See also M_LINE, M_LL2XY

% Rich Pawlowicz (rich@ocgy.ubc.ca) 3/Sep/98
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.

%  10/Mar/99 - changed order of calls ('c' not handled correctly in mu_coast otherwise)
% 6/Nov/00 - eliminate returned stuff if ';' neglected (thx to D Byrne)

[m,n]=size(long);

if m==1 && n>1
  h=mu_coast('vector',[long' lat';long(1) lat(1)],'patch',C,'tag','m_patch',varargin{:});
elseif m>1 && n==1
  h=mu_coast('vector',[long lat;long(1) lat(1)],'patch',C,'tag','m_patch',varargin{:});
else
  h=mu_coast('vector',[reshape([long;long(1,:);NaN+ones(1,n)],(m+2)*n,1),...
                     reshape([lat;lat(1,:);NaN+ones(1,n)],(m+2)*n,1)],'patch',C,'tag','m_patch',varargin{:});
end

if nargout==0
 clear h
end

author: goGPS-Project


text:
function fourthorder
%%
%     laplace^2 u = f;   [0,1]^2
%               u = g_D;
%       Du\cdot n = g_N.
%
%more information , you can check fourorderdata.m,biharmonicP1.m,biharmonicP2.m
%biharmonicP3.m.
% Created by Jie Zhou.
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.
close all
%% Parameters
%maxN = 2e3;     theta = 0.5;
maxIt = 6;
N = zeros(maxIt,1);     errL2 = zeros(maxIt,1);     errH1 = zeros(maxIt,1);

%%  Generate an initial mesh
node = [0,0; 0,1; 1,0;1,1];        % nodes
elem = [1,3,2; 4,2,3];                 % elements
elem = label(node,elem);               % label the mesh
bdEdge = setboundary(node,elem,'Dirichlet');    % Dirichlet boundary condition

% showmesh(node,elem);                            % plot mesh
% findelem(node,elem,'all','index','color','g');  % plot element indices
% findnode(node,'all','index','color','r');       % plot node indices

%%  Get a fine mesh by uniform bisection
for k = 1:2
     [node,elem,bdEdge] = uniformbisect(node,elem,bdEdge);
%    [node,elem,bdEdge] = bisect(node,elem,'all',bdEdge);
end

showmesh(node,elem);
findnode(node,'all','index','color','r');
findelem(node,elem,'all','index','color','r');

%% Set up PDE data
pde = fourorderdata;

%%       Adaptive Finite Element Method
% *SOLVE* -> *ESTIMATE* -> *MARK* -> *REFINE*
 for k = 1:maxIt
%                   [w,u] = biharmonicP1(node,elem,pde,bdEdge);
                    [w,u] = biharmonicP2(node,elem,pde,bdEdge);
%                   [w,u] = biharmonicP3(node,elem,pde,bdEdge);
               errL2(k) = getL2error(node,elem,pde.exactu,u);
               errH1(k) = getH1error(node,elem,pde.Du,u);
                   N(k) = size(w,1)+size(u,1);
     [node,elem,bdEdge] = uniformbisect(node,elem,bdEdge);
 end

%  Plot convergence rates
N= N(1:k);  errH1 = errH1(1:k);     errL2 = errL2(1:k);
figure;
showrate2(N,errH1,3,'-*','||Du-Du_h||',...
          N,errL2,3,'k-+','||u-u_h||');

repo: Aero-matlab


text:
function sample = sampleFromCategorical(probabilities)
% sampleFromCategorical   Sample from a categorical distribution
%
%   sample = sampleFromCategorical(probabilities) returns an index sampled
%   from the categorical distribution represented by the input
%   probabilities.
%
%   Input:
%       probabilities   - A numClasses-by-1 vector of probabilities. The
%                         elements of this vector should sum to 1.
%
%   Output:
%       sample          - A number between 1 and numClasses that is sampled
%                         from the input probabilities.

cdf = cumsum(probabilities);
sample = find( cdf > rand );
sample = sample(1);

end

author: matlab-deep-learning

repo: transformer-models


text:
function g = tanh(f)
%TANH   Hyperbolic tangent of a BALLFUN.
%   TANH(F) computes the hyperbolic tangent of the BALLFUN F.
%
% See also TAN, SINH, COSH.

% Copyright 2019 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

g = compose( f, @tanh );
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function results = vl_test_twister(varargin)
% VL_TEST_TWISTER
vl_test_init ;

function test_illegal_args()
vl_assert_exception(@() vl_twister(-1), 'vl:invalidArgument') ;
vl_assert_exception(@() vl_twister(1, -1), 'vl:invalidArgument') ;
vl_assert_exception(@() vl_twister([1, -1]), 'vl:invalidArgument') ;

function test_seed_by_scalar()
rand('twister',1) ; a = rand ;
vl_twister('state',1) ; b = vl_twister ;
vl_assert_equal(a,b,'seed by scalar + VL_TWISTER()') ;

function test_get_set_state()
rand('twister',1) ; a = rand('twister') ;
vl_twister('state',1) ; b = vl_twister('state') ;
vl_assert_equal(a,b,'read state') ;

a(1) = a(1) + 1 ;
vl_twister('state',a) ; b = vl_twister('state') ;
vl_assert_equal(a,b,'set state') ;

function test_multi_dimensions()
b = rand('twister') ;
rand('twister',b) ;
vl_twister('state',b) ;
a=rand([1 2 3 4 5]) ;
b=vl_twister([1 2 3 4 5]) ;
vl_assert_equal(a,b,'VL_TWISTER([M N P ...])') ;

function test_multi_multi_args()
rand('twister',1) ; a=rand(1, 2, 3, 4, 5) ;
vl_twister('state',1) ; b=vl_twister(1, 2, 3, 4, 5) ;
vl_assert_equal(a,b,'VL_TWISTER(M, N, P, ...)') ;

function test_square()
rand('twister',1) ; a=rand(10) ;
vl_twister('state',1) ; b=vl_twister(10) ;
vl_assert_equal(a,b,'VL_TWISTER(N)') ;

author: jianxiongxiao

repo: ProfXkit


text:
function [k,a]=invkini(x,y,phi)
global   l1 l2 l3
 l1=1;l2=1;l3=0.5;
xx=x-l3*cos(phi);
yy=y-l3*sin(phi);
k=(xx^2+yy^2-l1^2-l2^2)/(2*l1*l2);
if k > 1
    a=[];
    return
end
t2=acos(k);
d=l1^2+l2^2+2*l1*l2*cos(t2);
cc=(xx*(l1+l2*cos(t2))+yy*l2*sin(t2))/d;
ss=(-xx*l2*sin(t2)+yy*(l1+l2*cos(t2)))/d;
t1=atan2(ss,cc);
t3=phi-t1-t2;
a=[ t1 t2 t3];


text:
%   AUTHORSHIP
%   Math Lead & Primary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Secondary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Funded by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%
function [rows, cols, vals] = compute_membership_strengths(knn_indices, knn_dists, sigmas, rhos, same_set)
%COMPUTE_MEMBERSHIP_STRENGTHS Construct the membership strength data for
% the 1-skeleton of each local fuzzy simplicial set -- this is formed as a
% sparse matrix where each row is a local fuzzy simplicial set, with a
% membership strength for the 1-simplex to each other data point.
%
% [rows, cols, vals] = COMPUTE_MEMBERSHIP_STRENGTHS(knn_indices, knn_dists, sigmas, rhos, same_set)
%
% Parameters
% ----------
% knn_indices: array of size (n_samples, n_neighbors)
%     The indices on the "n_neighbors" closest points in the dataset.
%
% knn_dists: array of size (n_samples, n_neighbors)
%     The distances to the "n_neighbors" closest points in the dataset.
%
% sigmas: array of size (n_samples, 1)
%     The normalization factor derived from the metric tensor approximation.
%
% rhos: array of size (n_samples, 1)
%     The local connectivity adjustment.
%
% Returns
% -------
% rows: array of size (n_samples*n_neighbors, 1)
%     Row data for the resulting sparse matrix.
%
% cols: array of size (n_samples*n_neighbors, 1)
%     Column data for the resulting sparse matrix.
%
% vals: array of size (n_samples*n_neighbors, 1)
%     Entries for the resulting sparse matrix.

    if nargin < 5
        same_set = true;
    end

    [n_samples, n_neighbors] = size(knn_indices);

    knn_fail = knn_indices == -1;

    rows = repmat((1:n_samples)', 1, n_neighbors);
    rows(knn_fail) = NaN;
    rows = rows';
    rows = rows(:);

    cols = knn_indices;
    cols(knn_fail) = NaN;
    cols = cols';
    cols = cols(:);

    d = knn_dists - rhos;

    vals = exp(-max(0, d./repmat(sigmas, [1 n_neighbors])));
    if same_set
        itself = knn_indices == repmat((1:n_samples)', 1, n_neighbors);
        vals(itself) = 0;
    end
    vals(knn_fail) = NaN;
    vals = vals';
    vals = vals(:);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% (this is just an interface)
function varargout = mydatesod_inv (varargin)
    [varargout{1:nargout}] = mydatesodi (varargin{:});
end

%!test
%! % just check whether it runs:
%! mydatesod_inv(1, 2);

author: Sable


text:
function result = document_cache(context, identifier, fun, varargin)
% document_cache Cache proxy for report generation
%
% Can be used with report context to cache results of a certain function.
%
% Input:
% - context (structure): Report context structure.
% - identifier (string): Caching identifier.
% - fun (handle): Handle to processing function.
% - varargin (cell): Arguments for processing function.
%
% Output:
% - result: Output argument of the processing function.
%
cache_file = fullfile(context.cachedir, sprintf('%s.mat', identifier));

result = {};
if exist(cache_file, 'file')
    load(cache_file);
end;

if isempty(result)

    clear result;

    result = fun(varargin{:});

    save(cache_file, 'result');

end;

author: votchallenge

repo: toolkit-legacy


text:
function value = r8_tan ( x )

%*****************************************************************************80
%
%% R8_TAN evaluates the tangent of an R8 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the tangent of X.
%
  persistent nterms
  persistent pi2rec
  persistent sqeps
  persistent tancs
  persistent xmax
  persistent xsml

  pi2rec = 0.011619772367581343075535053490057;

  if ( isempty ( nterms ) )
    tancs = [ ...
      +0.22627932763129357846578636531752, ...
      +0.43017913146548961775583410748067E-01, ...
      +0.68544610682565088756929473623461E-03, ...
      +0.11045326947597098383578849369696E-04, ...
      +0.17817477903926312943238512588940E-06, ...
      +0.28744968582365265947529646832471E-08, ...
      +0.46374854195902995494137478234363E-10, ...
      +0.74817609041556138502341633308215E-12, ...
      +0.12070497002957544801644516947824E-13, ...
      +0.19473610812823019305513858584533E-15, ...
      +0.31417224874732446504614586026666E-17, ...
      +0.50686132555800153941904891733333E-19, ...
      +0.81773105159836540043979946666666E-21, ...
      +0.13192643412147384408951466666666E-22, ...
      +0.21283995497042377309866666666666E-24, ...
      +0.34337960192345945292800000000000E-26, ...
      +0.55398222121173811200000000000000E-28, ...
      +0.89375227794352810666666666666666E-30, ...
      +0.14419111371369130666666666666666E-31 ]';
    nterms = r8_inits ( tancs, 19, 0.1 * r8_mach ( 3 ) );
    xmax = 1.0 / r8_mach ( 4 );
    xsml = sqrt ( 3.0 * r8_mach ( 3 ) );
    sqeps = sqrt ( r8_mach ( 4 ) );
  end

  y = abs ( x );

  if ( xmax < y )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_TAN - Warning!\n' );
    fprintf ( 1, '  No precision because |X| is big.\n' );
    value = 0.0;
    return
  end

  ainty = r8_aint ( y );
  yrem = y - ainty;
  prodbg = 0.625 * ainty;
  ainty = r8_aint ( prodbg );
  y = ( prodbg - ainty ) + 0.625 * yrem + pi2rec * y;
  ainty2 = r8_aint ( y );
  ainty = ainty + ainty2;
  y = y - ainty2;

  ifn = r8_aint ( mod ( ainty, 2.0 ) );

  if ( ifn == 1 )
    y = 1.0 - y;
  end

  if ( 1.0 - y < abs ( x ) * sqeps )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_TAN - Warning!\n' );
    fprintf ( 1, '  Answer < half precision.\n' );
    fprintf ( 1, '  |X| big or X near pi/2 or 3*pi/2.\n' );
  end

  if ( y == 1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_TAN - Fatal error!\n' );
    fprintf ( 1, '  X is pi/2 or 3*pi/2.\n' );
    value = 0.0
    error ( 'R8_TAN - Fatal error!' )
  end

  if ( y <= 0.25 )

    value = y;
    if ( xsml < y )
      value = y * ( 1.5 + r8_csevl ( 32.0 * y * y - 1.0, tancs, nterms ) );
    end

  elseif ( y <= 0.5 )

    value = 0.5 * y * ( 1.5 + r8_csevl ( ...
      8.0 * y * y - 1.0, tancs, nterms ) );
    value = 2.0 * value / ( 1.0 - value * value );

  else

    value = 0.25 * y * ( 1.5 + r8_csevl ( ...
      2.0 * y * y - 1.0, tancs, nterms ) );
    value = 2.0 * value / ( 1.0 - value * value );
    value = 2.0 * value / ( 1.0 - value * value );

  end

  if ( x < 0.0 )
    value = - abs ( value );
  elseif ( 0.0 < x )
    value = + abs ( value );
  end

  if ( ifn == 1 )
    value = - value;
  end

  return
end


text:
% contourlet2tree.m
% written by: Duncan Po
% Date: August 24, 2002
% convert the contourlet coefficients from subband structure to tree structure
% Usage: [tree, scaling] = contourlet2tree(coef, dir)
% Inputs:   coef        - contourlet coefficients
%           dir         - the directional coefficients to be processed
% Output:   tree        - the resulting tree structure
%           scaling     - the scaling function of the coefficients

function [tree, scaling] = contourlet2tree(coef, dir)

M = length(coef{2});
nlevel = length(coef)-1;

for l = 2:nlevel+1
    levndir(l-1) = log2(length(coef{l}));
end;

scaling = coef{1};

for r = 1:(length(levndir)-1)
    if levndir(r+1)==levndir(r)
        split(r) = 0;
    elseif levndir(r+1)-levndir(r)==1
        split(r) = 1;
    else
        return;
    end;
end;

for l=2:nlevel+1
    if dir> M/2
      for i=1+2^(levndir(l-1)-levndir(1))*(dir-1):2^(levndir(l-1)-levndir(1))*dir
          coef{l}{i} = coef{l}{i}.';
      end;
    end;
end;

for l=2:nlevel+1
    if (l>2) & (levndir(l-1)~=levndir(1))
        j=1;
        if split(l-2)==1
            for i=1+2^(levndir(l-1)-levndir(1))*(dir-1):2:2^(levndir(l-1)-levndir(1))*dir
                coef{l}{j} = type3transform(coef{l}{i},coef{l}{i+1});
                j = j + 1;
            end;
        else
            for i=1+2^(levndir(l-1)-levndir(1))*(dir-1):2:2^(levndir(l-1)-levndir(1))*dir
                coef{l}{j} = type4transform(coef{l}{i},coef{l}{i+1});
                j = j + 1;
            end;
        end;

        i=length(coef{l})/length(coef{2})/2;
        while i>1
            j=1;
            for k=1:2:i
                coef{l}{j} = type4transform(coef{l}{k}, coef{l}{k+1});
                j = j+1;
            end;
            i = i/2;
        end;
        interimstructure{l-1} = coef{l}{1};
    else;
        interimstructure{l-1} = coef{l}{dir};
    end;
end;

tree{1} = interimstructure{1}(:);
for l = 2:nlevel
    tree{l} = [];
    for col = 1:2:size(interimstructure{l}, 2)
        temp = interimstructure{l}(:,col:col+1);
        temp = temp.';
        tree{l} = [tree{l}; temp(:)];
    end;
end;

repo: mcbench-benchmarks


text:
function startup

if ~ isdeployed
    % add root
    addpath(fileparts(mfilename('fullpath')),0);
end

% startup MTEX
startup_mtex;

author: mtex-toolbox

repo: mtex

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


text:
function unicycle_test07 ( )

%*****************************************************************************80
%
%% UNICYCLE_TEST07 tests UNICYCLE_UNRANK.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 June 2012
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'UNICYCLE_TEST07\n' );
  fprintf ( 1, '  UNICYCLE_UNRANK, given a rank, computes the\n' );
  fprintf ( 1, '  corresponding unicycle.\n' );
  fprintf ( 1, '\n' );
  rank = 6;
  fprintf ( 1, '  The requested rank is %d\n', rank );

  u = unicycle_unrank ( n, rank );

  unicycle_print ( n, u, '  The unicycle:' );

  return
end


text:
% -------------------------------------------------------------------------------------------------
function [bboxes,fps, ground_truth] = tracker2(varargin)
%TRACKER
%   is the main function that performs the tracking loop
%   Default parameters are overwritten by VARARGIN
%
%   Luca Bertinetto, Jack Valmadre, Joao F. Henriques, 2016
% -------------------------------------------------------------------------------------------------
    % These are the default hyper-params for SiamFC-3S
    % The ones for SiamFC (5 scales) are in params-5s.txt
    p.numScale = 3;
    p.scaleStep = 1.0375;
    p.scalePenalty = 0.9745;
    p.scaleLR = 0.59; % damping factor for scale update
    p.responseUp = 16; % upsampling the small 17x17 response helps with the accuracy
    p.windowing = 'cosine'; % to penalize large displacements
    p.wInfluence = 0.176; % windowing influence (in convex sum)
    p.net = '2016-08-17.net.mat';
    %% execution, visualization, benchmark
    p.video = 'vot15_bag';
    p.visualization = false;
    p.gpus = 1;
    p.bbox_output = false;
    p.fout = -1;
    %% Params from the network architecture, have to be consistent with the training
    p.exemplarSize = 127;  % input z size
    p.instanceSize = 255;  % input x size (search region)
    p.scoreSize = 17;
    p.totalStride = 8;
    p.contextAmount = 0.5; % context amount for the exemplar
    p.subMean = false;
    %% SiamFC prefix and ids
    p.prefix_z = 'a_'; % used to identify the layers of the exemplar
    p.prefix_x = 'b_'; % used to identify the layers of the instance
    p.prefix_join = 'xcorr';
    p.prefix_adj = 'adjust';
    p.id_feat_z = 'a_feat';
    p.id_score = 'score';
    % Get environment-specific default paths.
    p.net_base_path = '../models/';
    p.seq_base_path = '../demo-sequences/';
    p.seq_vot_base_path = '/path/to/VOT/evaluation/sequences/'; % (optional)
    p.stats_path = '/home/ethan/dxp/ours/Siam_adapt/models/ILSVRC2015.stats.mat'; % (optional)
    % added parameters
    p.load_video_info = @load_video_info_kcf;
    % Overwrite default parameters with varargin
    p = vl_argparse(p, varargin);
% -------------------------------------------------------------------------------------------------

load_video_info_fun = p.load_video_info;
% p.video = video;
%     % Get environment-specific default paths.
%     p = env_paths_tracking(p);
    % Load ImageNet Video statistics
    if exist(p.stats_path,'file')
        stats = load(p.stats_path);
    else
        warning('No stats found at %s', p.stats_path);
        stats = [];
    end
    % Load two copies of the pre-trained network
    net_z = load_pretrained([p.net_base_path p.net], p.gpus);
    net_x = load_pretrained([p.net_base_path p.net], []);
    [imgFiles, targetPosition, targetSize, ground_truth] = load_video_info_fun(p.seq_base_path, p.video);
    nImgs = numel(imgFiles);
    startFrame = 1;
    % Divide the net in 2
    % exemplar branch (used only once per video) computes features for the target
    remove_layers_from_prefix(net_z, p.prefix_x);
    remove_layers_from_prefix(net_z, p.prefix_join);
    remove_layers_from_prefix(net_z, p.prefix_adj);
    % instance branch computes features for search region x and cross-correlates with z features
    remove_layers_from_prefix(net_x, p.prefix_z);
    zFeatId = net_z.getVarIndex(p.id_feat_z);
    scoreId = net_x.getVarIndex(p.id_score);
    % get the first frame of the video
    im = gpuArray(single(imgFiles{startFrame}));
    % if grayscale repeat one channel to match filters size
	if(size(im, 3)==1)
        im = repmat(im, [1 1 3]);
    end
    % Init visualization
    videoPlayer = [];
    if p.visualization && isToolboxAvailable('Computer Vision System Toolbox')
        videoPlayer = vision.VideoPlayer('Position', [100 100 [size(im,2), size(im,1)]+30]);
    end
    % get avg for padding
    avgChans = gather([mean(mean(im(:,:,1))) mean(mean(im(:,:,2))) mean(mean(im(:,:,3)))]);

    wc_z = targetSize(2) + p.contextAmount*sum(targetSize);
    hc_z = targetSize(1) + p.contextAmount*sum(targetSize);
    s_z = sqrt(wc_z*hc_z);
    scale_z = p.exemplarSize / s_z;
    % initialize the exemplar
    [z_crop, ~] = get_subwindow_tracking(im, targetPosition, [p.exemplarSize p.exemplarSize], [round(s_z) round(s_z)], avgChans);
    if 0
        imwrite(uint8(gather(z_crop)),'test_z.png');
    end
    if p.subMean
        z_crop = bsxfun(@minus, z_crop, reshape(stats.z.rgbMean, [1 1 3]));
    end
    d_search = (p.instanceSize - p.exemplarSize)/2;
    pad = d_search/scale_z;
    s_x = s_z + 2*pad;
    % arbitrary scale saturation
    min_s_x = 0.2*s_x;
    max_s_x = 5*s_x;

    switch p.windowing
        case 'cosine'
            window = single(hann(p.scoreSize*p.responseUp) * hann(p.scoreSize*p.responseUp)');
        case 'uniform'
            window = single(ones(p.scoreSize*p.responseUp, p.scoreSize*p.responseUp));
    end
    % make the window sum 1
    window = window / sum(window(:));
    scales = (p.scaleStep .^ ((ceil(p.numScale/2)-p.numScale) : floor(p.numScale/2)));
    % evaluate the offline-trained network for exemplar z features
    net_z.eval({'exemplar', z_crop});
    z_features = net_z.vars(zFeatId).value;
    z_features = repmat(z_features, [1 1 1 p.numScale]);

    bboxes = zeros(nImgs, 4);
    % start tracking
time = 0;
    tic;
    for i = startFrame:nImgs
        if i>startFrame
            % load new frame on GPU
            im = gpuArray(single(imgFiles{i}));
   			% if grayscale repeat one channel to match filters size
    		if(size(im, 3)==1)
        		im = repmat(im, [1 1 3]);
    		end
            scaledInstance = s_x .* scales;
            scaledTarget = [targetSize(1) .* scales; targetSize(2) .* scales];
            % extract scaled crops for search region x at previous target position
            x_crops = make_scale_pyramid(im, targetPosition, scaledInstance, p.instanceSize, avgChans, stats, p);
            if 0
                imwrite(uint8(gather(x_crops(:,:,:,2))),'test_x.png');
            end
            % evaluate the offline-trained network for exemplar x features
%             tic
            [newTargetPosition, newScale] = tracker_eval(net_x, round(s_x), scoreId, z_features, x_crops, targetPosition, window, p);
%             toc
            targetPosition = gather(newTargetPosition);
            % scale damping and saturation
            s_x = max(min_s_x, min(max_s_x, (1-p.scaleLR)*s_x + p.scaleLR*scaledInstance(newScale)));
            targetSize = (1-p.scaleLR)*targetSize + p.scaleLR*[scaledTarget(1,newScale) scaledTarget(2,newScale)];
        else
            % at the first frame output position and size passed as input (ground truth)
        end

        rectPosition = [targetPosition([2,1]) - targetSize([2,1])/2, targetSize([2,1])];
        % output bbox in the original frame coordinates
        oTargetPosition = targetPosition; % .* frameSize ./ newFrameSize;
        oTargetSize = targetSize; % .* frameSize ./ newFrameSize;
        bboxes(i, :) = [oTargetPosition([2,1]) - oTargetSize([2,1])/2, oTargetSize([2,1])];

        if p.visualization
            if isempty(videoPlayer)
                figure(1), imshow(im/255);
                figure(1), rectangle('Position', rectPosition, 'LineWidth', 4, 'EdgeColor', 'y');
                drawnow
                fprintf('Frame %d\n', startFrame+i);
            else
                im = gather(im)/255;
                im = insertShape(im, 'Rectangle', rectPosition, 'LineWidth', 4, 'Color', 'yellow');
                % Display the annotated video frame using the video player object.
                step(videoPlayer, im);
            end
        end

        if p.bbox_output
            fprintf(p.fout,'%.2f,%.2f,%.2f,%.2f\n', bboxes(i, :));
        end

    end

    time = time + toc;
    bboxes = bboxes(startFrame : i, :);
    fps = nImgs/time;
end

author: shenjianbing

repo: TripletTracking


text:
function seg = segCloseMesh(seg);
% seg = segCloseMesh(seg);
%
% Close a mesh in a segmentation. Also closes any display windows.
%
% ras, 10/2006.
m = seg.settings.mesh;
if meshGet(seg.mesh{m}, 'id') > 0
    % close mesh display window
    mrmSet(seg.mesh{m}, 'Close');
end
ok = setdiff(1:length(seg.mesh), m);
seg.mesh = seg.mesh(ok);
if isempty(ok)
	seg.settings.mesh = 0;
else
	seg.settings.mesh = max(m-1, 1); % either prev. mesh in list, or mesh 1
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [f,fnorm]=setnorm(f,varargin)
%SETNORM  Set a norm of the input signal to a specified value
%   Usage:  h=setnorm(f,...);
%
%
%   `setnorm(f,...)` will set the input signal *f* to a specified norm.
%
%   `[f,fnorm]=setnorm(f,...)` does the same thing, but in addition
%   returns norm *fnorm* of a signal *f*.
%
%   The norm is specified as a string and may be one of:
%
%     '1'       Normalize the $l^1$ norm to be *1*.
%
%     'area'    Normalize the area of the signal to be *1*. This is exactly the same as `'1'`.
%
%     '2'       Normalize the $l^2$ norm to be *1*.
%
%     'energy'  Normalize the energy of the signal to be *1*. This is exactly
%               the same as `'2'`.
%
%     'inf'     Normalize the $l^{\inf}$ norm to be *1*.
%
%     'peak'    Normalize the peak value of the signal to be *1*. This is exactly
%               the same as `'inf'`.
%
%     'rms'     Normalize the Root Mean Square (RMS) norm of the
%               signal to be *1*.
%
%     's0'      Normalize the S0-norm to be *1*.
%
%     'wav'     Normalize to the $l^{\inf}$ norm to be *0.99* to avoid
%               possible clipping introduced by the quantization procedure
%               when saving as a wav file. This only works with floating
%               point data types.
%
%     'null'    Do NOT normalize, output is identical to input.
%
%
%   It is possible to specify the value to which the norm is set and the dimension:
%
%     'val',v    Value to set the norm to
%
%      'dim',d
%                Work along specified dimension. The default value of `[]`
%                means to work along the first non-singleton one.
%
%   See also: rms, s0norm

if ~isnumeric(f)
  error('%s: Input must be numerical.',upper(mfilename));
end;

if nargin<1
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.import={'setnorm'};
definput.keyvals.dim=[];
definput.keyvals.val=[];

[flags,kv]=ltfatarghelper({},definput,varargin,'setnorm');

if flags.do_null || flags.do_norm_notset || isempty(f);
  return
end;

if isa(f,'integer') && ~flags.do_wav
   error('%s: Integer data types are unsupported.',upper(mfilename));
end

if isempty(kv.val)
    kv.val=1;
end
%% ------ Computation --------------------------

[f,L,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,[],kv.dim, ...
                                                  upper(mfilename));
fnorm = zeros(W,1);

for ii=1:W

  if flags.do_1 || flags.do_area
    fnorm(ii) =  norm(f(:,ii),1);
  end;

  if flags.do_2 || flags.do_energy
    fnorm(ii) = norm(f(:,ii),2);
  end;

  if flags.do_inf || flags.do_peak
    fnorm(ii) = norm(f(:,ii),Inf);
  end;

  if flags.do_rms
    fnorm(ii) = rms(f(:,ii));
  end;

  if flags.do_s0
    fnorm(ii) = s0norm(f(:,ii));
  end;

  if flags.do_wav
    if isa(f,'float')
       fnorm(ii) = (1/0.99)*norm(f(:,ii),Inf);
    else
       error(['%s: TO DO: Normalizing integer data types not supported ',...
              'yet.'],upper(mfilename));
    end
  end;

  fnorm = 1/kv.val*fnorm;

  if fnorm(ii) > 0
     f(:,ii)=f(:,ii)/fnorm(ii);
  end
end;

f=assert_sigreshape_post(f,kv.dim,permutedsize,order);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function Data = normalize_data ( Data, options )

n = length(Data);

%% noramlizing the first image
shape          = normalize_first_shape( Data(1) , options );
Data(1).shape  = shape;

%% using the first to noramlizing others.
for i = 2 : n
    [shape] = normalize_rest_shape( Data(1), Data(i), options );
    Data(i).shape  = shape;
end

end

author: tntrung


text:
%	Bloch simulator using mex file
%   [mx,my,mz] = blochCim(b1,gr,tp,T1,T2,freq,pr,mode,sens,mx0,my0,mz0);
%
%	Bloch simulation of rotations due to B1, gradient and
%	off-resonance, including relaxation effects.  At each time
%	point, the rotation matrix and decay matrix are calculated.
%	Simulation can simulate the steady-state if the sequence
%	is applied repeatedly.This code can also simulate the parallel transmit.
%
%   INPUT:
%   b1 = (ntime * ncoil) RF pulse, complex number in (Gauss).
%   gr = (ntime * 3) gradient in (Gauss/cm).
%   tp = (ntime * 1) time duration of each b1 and gr point, in (sec),
%				or (1 * 1) time step if constant for all points
%				or monotonically INCREASING endtime of each
%				interval..
%   T1, T2 = (1 * 1) relaxation time in (sec)
%   freq = (nx * ny * nz) or (npos * 1) off-resonance freq. in (Hz)
%   pr = (npos * 3) array of spatial positions in (cm)
%   mode= Bitmask simulation mode: (default 0)
%		Bit 0:  0-Simulate from start, 1-Steady State
%		Bit 1:  1-Record m at time points.  0-just end time.%
%       Bit 2:  0-display simulation information. 1-only print out warnings.
%
%
%   (OPTIONAL)
%   sens = (ncoils * npos or ncoils * nx*ny*nz) sensitivities matrix.
%   default is 1 for all position. Set to default if you put 1 here.
%   mx0,my0,mz0 = (npos * 1) initial magnetization
%   default is (0,0,1) for all position
%
%	OUTPUT:
%		mx,my,mz = npos x 1 arrays of the resulting magnetization
%				components at each position.
%
%	The code is originally downloaded from Brian Hargreaves's website.
%
%   Hao Sun made the following main changes:
%   (1)match field map with position
%   (2)correct the minus sign error in the 'b1imag' in the blochsimfz sub
%   function in Brian's code
%   (3)add parallel simulation function
%
%
%   Hao Sun, the University of Michigan, Jul 22,2012

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% TRAIN_CNNSKETCH2IMAGERETRIEVAL Code to train the methods presented in the paper:
% F. Radenovic, G. Tolias, O. Chum, Deep Shape Matching, ECCV 2018
%
% Note: The method has been re-coded since our ECCV 2018 paper and minor differences in performance might appear.
%
% Authors: F. Radenovic, G. Tolias, O. Chum. 2018.

clear;

%-------------------------------------------------------------------------------
% Set data folder
%-------------------------------------------------------------------------------

% Set data folder, change if you have downloaded the data somewhere else
data_root = fullfile(get_root_cnnimageretrieval(), 'data');
% Check, and, if necessary, download train data (db with edgemaps), and pre-trained imagenet networks
download_train_sketch(data_root);

%-------------------------------------------------------------------------------
% Reproduce training from ECCV18 paper: Deep Shape Matching ...
%-------------------------------------------------------------------------------

% Set architecture and initialization parameters
opts.init.model = 'VGG'; % (ALEX | VGG | GOOGLENET | RESNET101)
opts.init.modelDir = fullfile(data_root, 'networks', 'imagenet');
opts.init.method = 'edgefilter_mac';
opts.init.objectiveType = {'contrastiveloss', 0.7};
opts.init.errorType = {'batchmap'};
opts.init.averageImageScale = 0;
opts.init.imageChannels = 1;

% Set train parameters
opts.train.dbPath = fullfile(data_root, 'train', 'dbs', 'retrieval-SfM-30k-edgemap.mat');
opts.train.batchSize = 20;
opts.train.numSubBatches = 4;
opts.train.numEpochs = 20;
opts.train.learningRate = 0.001 .* exp(-(0:99)*0.1);
opts.train.numNegative = 5;
opts.train.numRemine = 3;
opts.train.gpus = [1];

opts.train.augment.jitterFlip = true;
opts.train.augment.jitterQueryBinarize = true;

% Trial name (to name a save directory)
trialName = 'test';

% Export directory expDir named after model, method and trialName
opts.init.method = [opts.init.method, '_', trialName];
opts.train.expDir = fullfile(data_root, 'networks', 'exp', [lower(opts.init.model) '_' lower(opts.init.method)]);
if ~exist(opts.train.expDir); mkdir(opts.train.expDir); end % create folder if its not there

% Load opts by respecting added opts
opts = load_opts_train(opts);

% Initialize and train the network
fprintf('>> Experiment folder is set to %s\n', opts.train.expDir);
net = init_network(opts.init);
[net, state, stats] = train_network(net, @(o,i,n,b,s,m,e) get_batch(o,i,n,b,s,m,e), opts.train);

author: filipradenovic

repo: cnnimageretrieval


text:
function sphere_grid_test01 ( )

%*****************************************************************************80
%
%% SPHERE_GRID_TEST01 tests SPHERE_ICOS_POINT_NUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 October 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_GRID_TEST01\n' );
  fprintf ( 1, '  SPHERE_ICOS_POINT_NUM determines the size\n' );
  fprintf ( 1, '  (number of vertices, edges and faces) in a grid\n' );
  fprintf ( 1, '  on a sphere, made by subdividing an initial\n' );
  fprintf ( 1, '  projected icosahedron.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  N determines the number of subdivisions of each\n' );
  fprintf ( 1, '  edge of the icosahedral faces.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         N         V         E         F\n' );
  fprintf ( 1, '  --------  --------  --------  --------\n' );
  fprintf ( 1, '\n' );

  for factor = 1 : 20
    point_num = sphere_icos_point_num ( factor );
    edge_num = sphere_icos_edge_num ( factor );
    face_num = sphere_icos_face_num ( factor );
    fprintf ( 1, '  %8d  %8d  %8d  %8d\n', ...
      factor, point_num, edge_num, face_num );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Repeat, but using N constrained by doubling:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         N         V         E         F\n' );
  fprintf ( 1, '  --------  --------  --------  --------\n' );
  fprintf ( 1, '\n' );

  factor = 1;
  for factor_log = 0 : 10
    point_num = sphere_icos_point_num ( factor );
    edge_num = sphere_icos_edge_num ( factor );
    face_num = sphere_icos_face_num ( factor );
    fprintf ( 1, '  %8d  %8d  %8d  %8d\n', ...
      factor, point_num, edge_num, face_num );
    factor = factor * 2;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Demo for Structured Edge Detector (please see readme.txt first).

error(['This is a demo for Stuctured Edges and not for our RGB-D edge detector.', ...
  ' Please look at run_all.m for how to run our RGB-D edge detector and look at scripts/script_edges.m' ...
  ' for how to train the RGB-D edge detector.']);

%% set opts for training (see edgesTrain.m)
opts=edgesTrain();                % default options (good settings)
opts.modelDir='models/';          % model will be in models/forest
opts.modelFnm='modelBsds';        % model name
opts.nPos=5e5; opts.nNeg=5e5;     % decrease to speedup training
opts.useParfor=0;                 % parallelize if sufficient memory

%% train edge detector (~20m/8Gb per tree, proportional to nPos/nNeg)
tic, model=edgesTrain(opts); toc; % will load model if already trained

%% set detection parameters (can set after training)
model.opts.multiscale=0;          % for top accuracy set multiscale=1
model.opts.sharpen=2;             % for top speed set sharpen=0
model.opts.nTreesEval=4;          % for top speed set nTreesEval=1
model.opts.nThreads=4;            % max number threads for evaluation
model.opts.nms=0;                 % set to true to enable nms

%% evaluate edge detector on BSDS500 (see edgesEval.m)
if(0), edgesEval( model, 'show',1, 'name','' ); end

%% detect edge and visualize results
I = imread('peppers.png');
tic, E=edgesDetect(I,model); toc
figure(1); im(I); figure(2); im(1-E);

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function [ReactionAbundance,ReactionPresence] = fastCalculateReactionAbundance(abundancePath, modelPath, rxnsList, numWorkers)
% Part of the Microbiome Modeling Toolbox. This function calculates and
% plots the total abundance of reactions of interest in a given microbiome
% sample based on the strain-level composition.
% Reaction presence or absence in each strain is derived from the reaction content
% of the respective AGORA model.
%
% USAGE
%
%    ReactionAbundance = fastCalculateReactionAbundance(abundancePath, modelPath, rxnsList, numWorkers)
%
% INPUTS:
%    abundancePath:          Path to the .csv file with the abundance data.
%                            Example: 'cobratoolbox/papers/018_microbiomeModelingToolbox/examples/normCoverage.csv'
%    modelPath:              Folder containing the strain-specific AGORA models
% OPTIONAL INPUTS:
%    rxnsList:               List of reactions for which the abundance
%                            should be calculated (if left empty: all
%                            reactions in all models)
%    numWorkers:             Number of workers used for parallel pool. If
%                            left empty, the parallel pool will not be
%                            started. Parallellization is recommended if
%                            all reactions are computed.
%
% OUTPUT:
%    ReactionAbundance       Table with total abundance for each microbiome
%                            and reaction
%    ReactionPresence:       Table with absolute reaction presence for
%                            all reactions in all microbiome samples
%
% .. Author: - Almut Heinken, 04/2021

% read the csv file with the abundance data
abundance = readInputTableForPipeline(abundancePath);
if isnumeric(abundance{2, 1})
    abundance(:, 1) = [];
end

% load the models
for i = 2:size(abundance, 1)
    model = readCbModel([modelPath filesep abundance{i, 1} '.mat']);
    modelsList{i, 1} = model;
end

if ~exist('rxnsList', 'var') || isempty(rxnsList)  % define reaction list if not entered
    fprintf('No reaction list entered. Abundances will be calculated for all reactions in all models. \n')
    % get model list from abundance input file
    for i = 2:size(abundance, 1)
        model = modelsList{i, 1};
        rxnsList = vertcat(model.rxns, rxnsList);
    end
    rxnsList = unique(rxnsList);
end

% load the models found in the individuals and extract which reactions are
% in which model
for i = 2:size(abundance, 1)
    model = modelsList{i, 1};
    ReactionPresenceInModels{i, 1} = abundance{i, 1};
    for j = 1:length(rxnsList)
        ReactionPresenceInModels{1, j + 1} = rxnsList{j};
        if ~isempty(find(ismember(model.rxns, rxnsList{j})))
            ReactionPresenceInModels{i, j + 1} = '1';
        else
            ReactionPresenceInModels{i, j + 1} = '0';
        end
    end
end
ReactionPresenceInModels{1,1}='Strains';

% prepare table for the total abundance
ReactionAbundance = {};
for i = 2:size(abundance, 2)
    ReactionAbundance{i, 1} = abundance{1, i};
end

% use parallel pool if workers specified as input
if exist('numWorkers', 'var') && numWorkers > 0
    poolobj = gcp('nocreate');
    if isempty(poolobj)
        parpool(numWorkers)
    end
end

clear abundance
clear modelsList

totalAbun={};
parfor i = 2:size(ReactionAbundance, 1)
    % reload the file to avoid running out of memory
    abundance = readInputTableForPipeline(abundancePath);
    if isnumeric(abundance{2, 1})
        abundance(:, 1) = [];
    end

    % temporarily store reaction abundances
    totalAbun{i} = zeros(length(rxnsList), 1);

    for j = 2:size(abundance, 1)
        % find all reactions present in the strain
        presentRxns = find(strcmp(ReactionPresenceInModels(j,2:end),'1'));

        for k = 1:length(presentRxns)
            % summarize total abundance
            if contains(version,'(R202') % for Matlab R2020a and newer
                totalAbun{i}(presentRxns(k),1) = totalAbun{i}(presentRxns(k),1) + abundance{j,i};
            else
                totalAbun{i}(presentRxns(k),1) = totalAbun{i}(presentRxns(k),1) + str2double(abundance{j,i});
            end
        end
    end
end

for i = 1:length(rxnsList)
    ReactionAbundance{1, i + 1} = rxnsList{i};
end

% collect the temporarily stored abundances to put together the table
for i = 2:size(ReactionAbundance, 1)
    ReactionAbundance(i,2:end) = num2cell(totalAbun{i});
end

% export the absolute reaction presence
tol=0.0000001;

ReactionPresence = ReactionAbundance;
for i=2:size(ReactionPresence,1)
    for j=2:size(ReactionPresence,2)
        if contains(version,'(R202') % for Matlab R2020a and newer
            val = ReactionPresence{i,j};
        else
            val = str2double(ReactionPresence{i,j});
        end
        if val > tol
            ReactionPresence{i,j} = 1;
        else
            ReactionPresence{i,j} = 0;
        end
    end
end

end

author: opencobra

repo: cobratoolbox


text:
function outfilenames = SaveFigLotsOfWays(hfig,basename)

if exist('savefig1','var'),
  savefigfun = @savefig1;
else
  savefigfun = @savefig;
end

% formats = {'pdf','tiff','svg','fig','eps'};
formats = {'pdf','tiff','png','fig','eps'};
[path,basename] = myfileparts(basename);

outfilenames = {};
if ismember('svg',formats),
  filename = [basename,'.svg'];
  plot2svg(fullfile(path,filename),hfig,'png');
  outfilenames{end+1} = fullfile(path,filename);
end
formats = setdiff(formats,{'svg'});

if ismember('ai',formats),
  filename = [basename,'.ai'];
  saveas(hfig,fullfile(path,filename),'ai');
  outfilenames{end+1} = fullfile(path,filename);
end
formats = setdiff(formats,{'ai'});

if ismember('fig',formats),
  filename = [basename,'.fig'];
  saveas(hfig,fullfile(path,filename),'fig');
  outfilenames{end+1} = fullfile(path,filename);
end
formats = setdiff(formats,{'fig'});

tmpbasename = regexprep(basename,'[^a-zA-Z_0-9]','_');
for i = 1:numel(formats),
  filename = [basename,'.',formats{i}];
  tmpfilename = [tmpbasename,'.',formats{i}];
  try
    savefigfun(tmpfilename,hfig,formats{i});
    if ~isempty(path),
      movefile(tmpfilename,fullfile(path,filename));
    end
  catch ME,
    warning('Error while creating %s: %s',filename,getReport(ME));
    continue;
  end
  outfilenames{end+1} = fullfile(path,filename);
end

author: kristinbranson

repo: JAABA


text:
%BOX Collision checking class for box-type primitive shapes
%
% Subclass of the Shape class. The x, y and z axes of the transform
% describes the direction of the box's edges, and the translation
% component describes the centre of the box. The scale vector defines
% the edges' half-lengths.
%
% Copyright (C) Bryan Moutrie, 2013-2014
% Licensed under the GNU Lesser General Public License
% see full file for full statement
%
% Syntax:
%  (1) box = Box(T, s)
%  (2) box = Box(T, s, ...)
%  (3) box = Box()
%
%  (2) is as per (1) but with other properties set in name-value pairs
%  (3) returns a box with default values
%
% Outputs:
%  box : Box object
%
% Inputs:
%  T   : Transformation matrix. If empty, defaults to eye(4).
%  s   : Scale vector. If empty, defaults to [1 1 1].
%  ... : Options - other properties in name-value pairs
%
% See documentation for information on properties and methods
% (type doc Shape into the command window)
%
% See also Box.faces Box.n Cone Curvilinear Cylinder.Cylinder
% Ellipsoid.Ellipsoid Shape Sphere.Sphere

% LICENSE STATEMENT:
%
% This file is part of pHRIWARE.
%
% pHRIWARE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as
% published by the Free Software Foundation, either version 3 of
% the License, or (at your option) any later version.
%
% pHRIWARE 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 Lesser General Public
% License along with pHRIWARE.  If not, see <http://www.gnu.org/licenses/>.

classdef Box < Shape

    properties (Dependent)
        volume
    end

    properties
        %FACES  Logical vector of which faces to plot
        % For Boxes, faces is a 1x6 vector, whose elements correspond
        % to the faces in the negative y-z, x-z, x-y and positive y-z,
        % x-z, x-y planes, with respect to that shape's frame, in that
        % order
        faces = true(1,6)

        %N "Resolution" of the shape, for plotting
        % For Boxes, each face will comprise of an nxn grid of points
        n = 2
    end

    methods
        function box = Box(T, s, varargin)
            if ~nargin
                T = [];
                s = [];
            end

            box = box@Shape(T, s, varargin{:});
        end

        function vargout = plot(box, varargin)
            [box, frames] = parseopts(box,varargin{:});

            [X, Y, Z] = box.cloud;

            h = zeros(1,6);
            for f = 1: 6
                if box.faces(f)
                    h(f) = surface(X(:,:,f), Y(:,:,f), Z(:,:,f));
                end
            end
            h = h(box.faces);

            set(h, ...
                'FaceColor', box.FaceColor, ...
                'FaceAlpha', box.FaceAlpha, ...
                'EdgeColor', box.EdgeColor, ...
                'EdgeAlpha', box.EdgeAlpha);

            if ~isempty(frames), h = ellipsoid.animate(h, frames); end
            if nargout, vargout = h; end
        end

        function [X, Y, Z] = cloud(box)
            ptsAlongEdge = linspace(-1,1,box.n);
            [y, z] = meshgrid(ptsAlongEdge,ptsAlongEdge);
            x = -ones(size(z));
            xa = cat(3, x  , y  , z  , x+2, y  , z  );
            ya = cat(3, y  , x  , y  , y  , x+2, y  );
            za = cat(3, z  , z  , x  , z  , z  , x+2);

            [X, Y, Z] = box.sat(xa, ya, za);
        end

        function inside = check(box, varargin)
            switch length(varargin)
                case 0
                    if box.sym == 8
                        syms x y z real
                        points = [x, y, z, 1];
                    else
                        inside = @box.check;
                        return;
                    end
                case 1
                    points = [varargin{1},ones(size(varargin{1},1),1)];
                case 3
                    x = varargin{1}(:);
                    y = varargin{2}(:);
                    z = varargin{3}(:);
                    points = [x, y, z, ones(size(x))];
            end

            tpts = box.transform \ points';
            tpts = abs(tpts');

            inX = tpts(:,1) < box.scale(1);

            inY = tpts(:,2) < box.scale(2);

            inZ = tpts(:,3) < box.scale(3);

            inside = inX & inY & inZ;

            if isa(inside,'sym'), inside = sym2func(inside); end
        end

        function v = get.volume(box)
            v = 8*prod(box.scale);
        end

        function box = set.faces(box, f)
            box.faces = f & [1 1 1 1 1 1];
        end

    end

end

author: Allopart

repo: rbpf-gmapping


text:
%CODEGENERATOR.GENFKINE Generate code for forward kinematics
%
% T = cGen.genfkine() generates a symbolic homogeneous transform matrix (4x4) representing
% the pose of the robot end-effector in terms of the symbolic joint coordinates q1, q2, ...
%
% [T, ALLT] = cGen.genfkine() as above but also generates symbolic homogeneous transform
% matrices (4x4xN) for the poses of the individual robot joints.
%
% Notes::
% - Side effects of execution depends on the cGen flags:
%   - saveresult: the symbolic expressions are saved to
%     disk in the directory specified by cGen.sympath
%   - genmfun: ready to use m-functions are generated and
%     provided via a subclass of SerialLink stored in cGen.robjpath
%   - genslblock: a Simulink block is generated and stored in a
%     robot specific block library cGen.slib in the directory
%     cGen.basepath
%   - genccode: generates C-functions and -headers in the directory
%     specified by the ccodepath property of the CodeGenerator object.
%   - mex: generates robot specific MEX-functions as replacement for the
%     m-functions mentioned above. Access is provided by the SerialLink
%     subclass. The MEX files rely on the C code generated before.
%
% Author::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.CodeGenerator, CodeGenerator.geninvdyn, CodeGenerator.genjacobian.

% Copyright (C) 2012-2014, by Joern Malzahn
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

function [t,allT] = genfkine(CGen)

%% Derivation of symbolic expressions
CGen.logmsg([datestr(now),'\tDeriving forward kinematics']);

q = CGen.rob.gencoords;
[t, allT] = CGen.rob.fkine(q);

CGen.logmsg('\t%s\n',' done!');

%% Save symbolic expressions
if CGen.saveresult
    CGen.logmsg([datestr(now),'\tSaving symbolic forward kinematics up to end-effector frame']);

    CGen.savesym(t,'fkine','fkine.mat')

    CGen.logmsg('\t%s\n',' done!');

    CGen.logmsg([datestr(now),'\tSaving symbolic forward kinematics for joint']);

    for iJoint = 1:CGen.rob.n
        CGen.logmsg(' %s ',num2str(iJoint));
        tName = ['T0_',num2str(iJoint)];
        eval([tName,' = allT(',num2str(iJoint),');']);
        CGen.savesym(eval(tName),tName,[tName,'.mat']);
    end

    CGen.logmsg('\t%s\n',' done!');
end

%% M-Functions
if CGen.genmfun
    CGen.genmfunfkine;
end

%% Embedded Matlab Function Simulink blocks
if CGen.genslblock
    genslblockfkine(CGen);
end

%% C-Code
if CGen.genccode
    CGen.genccodefkine;
end

%% MEX
if CGen.genmex
    CGen.genmexfkine;
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function output = prtPath(varargin)
% prtPath Adds necessary directories for the PRT to your path.

if isdeployed
    return;
end

addpath(fullfile(prtRoot,'util'));
origPath = prtUtilGenCleanPath(prtRoot);
addpath(origPath);

for iArg = 1:length(varargin)
    cArg = varargin{iArg};
    cDir = fullfile(prtRoot,cat(2,']',cArg));
    assert(logical(exist(cDir,'file')),']%s is not a directory in %s',cArg,prtRoot);
    P = prtUtilGenCleanPath(cDir,'removeDirStart',{'.'});
    addpath(P);
    origPath = cat(2,origPath,P);
end
if nargout > 0
    output = origPath;
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT