text:
%SOLVER_INITIALIZE	TFOCS helper script
%	Performs the initializations common to all of the first-order solvers

% Process the options structure and string/value pairs. First, we replace the
% default values with any values specified by the user in the 'opts' structure.
% We ignore any fields in the 'opts' structure that do not match ours, so that
% users can re-use the opts structure for other purposes.

odef = struct( ...
    'maxIts',     Inf   ,...
    'maxCounts',  Inf   ,...
    'countOps',   false ,... % adjusted if maxCounts < Inf
    'saveHist',   true  ,...
    'adjoint',    false ,...
    'saddle',     false ,...
    'tol',        1e-8  ,...
    'errFcn',     {{}}  ,...
    'stopFcn',    {{}}  ,...
    'printEvery', 100   ,...
    'maxmin',     1     ,...
    'beta',       0.5   ,...
    'alpha',      0.9   ,... % See below for special CG stuff
    'L0',         1     ,... % See below
    'Lexact',     Inf   ,...
    'mu',         0     ,...
    'fid',        1     ,...
    'stopCrit',   1     ,... % the TFOCS paper used stopCrit = 2
    'alg',        'AT'  ,...
    'restart',    Inf   ,...
    'printStopCrit', false,...
    'cntr_reset',  -50   ,... % how often to explicitly recompute A*x and A*y (set to Inf if you want )
    'cg_restart', Inf   ,... % for CG only
    'cg_type',    'pr'  ,...  % for CG only
    'stopping_criteria_always_use_x',   false, ...
    'data_collection_always_use_x',     false, ...
    'output_always_use_x',              false,  ...
    'autoRestart',  'gra', ... % function or gradient
    'printRestart', true, ...
    'debug',      false ....
    );

% Calling the solver with a no arguments returns the default options structure
if narginn < 1 || ( narginn ==1 && isstruct( smoothF ) )
    opts = odef;
    % remove the "CG" options for now, since they are undocumented
    opts    = rmfield(opts,'cg_restart');
    opts    = rmfield(opts,'cg_type');
    % if any default options are passed in (e.g. by tfocs_SCD), add those now:
    if ( narginn ==1 && isstruct( smoothF ) )
        for f = fieldnames( smoothF )'
            opts.( f{1} )   = smoothF.( f{1} );
        end
    end

    out = [];
    x = opts;
    if nargoutt == 0,
        disp('====== TFOCS default options =======');
%         disp('Format is    fieldname: { [default]  ''Usage type.
%         Description''}');

        % add some helpful description
        desc = [];
        desc.alg    = 'Algorithm, e.g. GRA for gradient descent. Values: AT, GRA, LLM, N07, N83, TS';
        desc.maxIts = 'Basic. Maximum number of allowed iterations';
        desc.saveHist   = 'Basic. Record history of all iterations';
        desc.restart    = 'Basic. Restart parameter. Make negative for "no-regress" feature';
        desc.tol    = 'Basic. Tolerance for stopping criteria';
        desc.L0     = 'Basic. Initial estimate of Lipschitz constant';
        desc.Lexact = 'Basic. Known bound of Lipschitz constant';
        desc.printEvery = 'Basic. How often to print info to the terminal; set to 0 or Inf to suppress output';
        desc.maxmin = 'Basic. +1 for convex minimization, -1 for concave maximization';

        desc.errFcn = 'Medium. User-specified error function. See user guide';
        desc.beta   = 'Medium. Backtracking parameter, in (0,1). No line search if >= 1';
        desc.alpha  = 'Medium. Line search increase parameter, in (0,1)';
        desc.autoRestart= 'Medium. Set to ''gra'' or ''fun'' to choose behavior when restart<0';

        desc.maxCounts  = 'Advanced. Vector that fine-tunes various types of iteration limits; same form as countOps';
        desc.countOps   = 'Advanced. Record number of linear multiplies, etc.; form [fcn, grad, linear, nonsmth, proj]';
        desc.mu         = 'Advanced. Strong convexity parameter.';
        desc.fid        = 'Advanced. File id, e.g. via fopen. All output is sent to this file';
        desc.stopFcn    = 'Advanced. User-supplied stopping criteria. See user guide';
        desc.stopCrit   = 'Advanced. Controls which stopping criteria to use; 1,2, 3 or 4.';
        desc.printStopCrit = 'Advanced. Controls whether to display the value used in the stopping criteria';
        desc.printRestart  = 'Advanced. Whether to signal when a restart happens';
        desc.cntr_reset = 'Advanced. Controls how often to reset some numerical computations to avoid roundoff';
        desc.debug = 'Advanced.  Turns on more useful error messages';

        desc.stopping_criteria_always_use_x = 'Advanced. Forces usage of x, never y, in stopping crit.';
        desc.data_collection_always_use_x   = 'Advanced. Forces usage of x, nevery y, in recording errors.';
        desc.output_always_use_x            = 'Advanced. Forces output of x, never y. Default: uses whichever is better';

        desc.adjoint = 'Internal.';
        desc.saddle = 'Internal.  Used by TFOCS_SCD';

        disp( opts );
        disp('====== Description of TFOCS options =======');
        disp( desc );
    end
    return % from this script only; the calling function does not return
end

% [smoothF, affineF, projectorF, x0, opts ] = deal(varargin{:});
% Check for incorrect types
F_types = {'function_handle','cell','double','single'};
assert( ismember(class(smoothF),F_types),'smoothF is of wrong type' );
assert( ismember(class(affineF),F_types),'affineF is of wrong type' );
assert( ismember(class(projectorF),F_types),'projectorF 3 is of wrong type' );
if narginn >= 4
    x0_types = {'cell','double','single','tfocs_tuple'};
    assert( ismember(class(x0),x0_types),'x0 is of wrong type' );
end
if narginn >= 5 && ~isempty(opts)
    assert( ismember(class(opts),{'struct'}),'opts is of wrong type' );
end

% Some parameters defaults depend on whether the user supplies other
% options. These will be updated later
L0_default  = odef.L0;
odef.L0     = Inf;

alpha_default    = odef.alpha; % typically 0.9
alpha_default_CG = 1;
odef.alpha  = Inf;

% Process the options structure, merging it with the default options
%   (merge "opts" into "odef")
def_fields = fieldnames(odef)';              % default options
if narginn > 4 && ~isempty(opts),
    use_fields = zeros(size(def_fields));
    opt_fields = fieldnames(opts)';          % user-supplied options
    for k = opt_fields,
        k = k{1};
        ndx = find(strcmpi(def_fields,k));
        if ~isempty(ndx)
            if ~isempty(opts.(k)) && ~use_fields(ndx),
                odef.(def_fields{ndx}) = opts.(k);
            end
            use_fields(ndx) = use_fields(ndx) + 1;
        else
            % Warn if the field is not found in our options structure
            warnState = warning('query','backtrace');
            warning off backtrace
            warning('TFOCS:OptionsStructure',' Found extraneous field "%s" in options structure', k );
            warning(warnState);
        end
    end
    % Warn if fields appear twice due to capitalization; e.g., maxIts/maxits
    if any(use_fields>1),
        warnState = warning('query','backtrace');
        warning off backtrace
        warning('TFOCS:OptionsStructure',' Some fieldnames of the options structure are identical up to capitalization: unpredictable behavior');
        warning(warnState);
    end
end

% Remove unnecessary options
if ~strcmpi( odef.alg, 'cg' )
    odef = rmfield( odef, 'cg_restart' );
    odef = rmfield( odef, 'cg_type' );
    % update our list of fieldnames
    def_fields = fieldnames(odef)';
end

% If opts.alpha wasn't supplied, use a default value:
if isinf(odef.alpha)
    if strcmpi( odef.alg, 'cg' )
        odef.alpha = alpha_default_CG;
    else
        odef.alpha = alpha_default;
    end
end

% If opts.printEvery is Inf, set it to 0
if isinf(odef.printEvery)
    odef.printEvery = 0;
end

% The default value of L0 is set to Lexact, if it is supplied
if isinf(odef.L0),
    if isinf(odef.Lexact),
        if odef.beta >= 1,
            error( 'For a fixed step size, L0 or Lexact must be supplied.' );
        end
        odef.L0 = L0_default;
    else
        odef.L0 = odef.Lexact;
    end
end
% If maxCounts is set, set countOps to true
if any(odef.maxCounts<Inf),
    odef.countOps = true;
end
% If cntr_reset is not set (signfied by being negative), set to defaults
if odef.cntr_reset < 0
    odef.cntr_reset = round(abs(odef.cntr_reset));
    % and if we requested high precision, change the default
    if odef.tol < 1e-12
        odef.cntr_reset = 10;
    end
end

% Now move the options into the current workspace
for k = def_fields,
    assignin('caller',k{1},odef.(k{1}));
end
opts = odef;

%
% Smooth function, pass 1
%

if isempty( smoothF ),
	error( 'Must supply a smooth function specification.' );
elseif isa( smoothF, 'function_handle' ),
    smoothF = { smoothF };
elseif ~isa( smoothF, 'cell' ),
    error( 'smoothF must be a function handle, or a cell array of them.' );
end
n_smooth = numel(smoothF);
saddle_ndxs = 1 : n_smooth;
% Adding Feb 2011: a get_dual function (see below)

%
% Projector, pass 1
%

if isa( projectorF, 'function_handle' ),
    projectorF = { projectorF };
elseif ~isa( projectorF, 'cell' ) && ~isempty( projectorF ),
    error( 'projectorF must be a function handle, or a cell array of them.' );
end
n_proj = numel(projectorF);

%
% Linear functions and affine offsets
%

% If the affine operator is anything *but* a cell array, convert it to one.
maxmin = sign(maxmin);
if isempty( affineF )
    affineF = { @linop_identity };
elseif isnumeric( affineF ),
    if ndims(affineF) > 1,
        error( 'Multidimensional arrays are not permitted.' );
    end
    if numel(affineF) == 1,
        identity_linop = affineF == 1;
        affineF = { linop_scale( affineF ) };
    else
        identity_linop = false;
        affineF = { linop_matrix( affineF ) };
    end
elseif isa( affineF, 'function_handle' ),
    affineF = { affineF };
elseif ~isa( affineF, 'cell' ),
    error( 'Invalid affine operator specification.' );
end

% If adjoint mode is specified, temporarily transpose the affineF cell
% array so that the rows and columns match the number of smooth functions
% and projectors, respectively. Then verify size compatibility.
if adjoint, affineF = affineF'; end
[ m_aff, n_aff ] = size( affineF );
% if all( m_aff ~= n_smooth + [0,1] ) || n_proj && all( n_aff ~= n_proj + [0,1] ),
%     error( 'The affine operation matrix has incompatible dimensions.' );
% May 16, 2011: making error messages more informative
if all( m_aff ~= n_smooth + [0,1] )
    if fid
        fprintf(fid,'Detected error: inputs are of wrong size\n');
        fprintf(fid,'  Found %d smooth functions, so expect that many primal variables\n',...
            n_smooth );
        fprintf(fid,'  So the affine operator should have %d (or %d if there is an offset) entries\n',...
            n_smooth, n_smooth+1 );
        fprintf(fid,'  but instead found %d affine operators\n', m_aff );
        fprintf(fid,' Perhaps you need the transpose of the affine operator?\n');
    end
    error( 'The affine operation matrix has incompatible dimensions.' );
elseif n_proj && all( n_aff ~= n_proj + [0,1] ),
    if fid
        fprintf(fid,'Detected error: inputs are of wrong size\n');
        fprintf(fid,'  Found %d nonsmooth functions, so expect %d or %d sets of affine operators\n',...
            n_proj, n_proj, n_proj + 1 );
        fprintf(fid,'  but instead found %d affine operators\n', n_aff );
    end
    error( 'The affine operation matrix has incompatible dimensions.' );
elseif n_proj == 0,
    n_proj = max( 1, m_aff - 1 );
end
inp_dims = cell( 1, n_proj );
otp_dims = cell( 1, n_smooth );

% If an additional affine portion of the objective has been specified,
% create an additional smooth function to contain it.
if m_aff == n_smooth + 1,
    offset          = true; % note: saddle_ndxs is 1:n_smooth
    otp_dims{end+1} = [1,1];
    smoothF{end+1}  = smooth_linear( maxmin );
    n_smooth = n_smooth + 1;
    for k = 1 : n_proj,
        offX = affineF{end,k};
        if isempty(offX),
            affineF{end,k} = 0;
        elseif isa(offX,'function_handle'),
            if adjoint, pos = 'row'; else pos = 'column'; end
            error( 'The elements in the last %s must be constants.', pos );
        elseif isnumeric(offX) && numel(offX) == 1 && offX == 0,
            affineF{end,k} = 0;
        elseif nnz(offX),
            affineF{end,k} = linop_dot( affineF{end,k}, adjoint );
            % add a case if ~nnz(offX)... Jan 2012
        elseif ~nnz(offX)
            affineF{end,k} = 0;
        end
    end
else
    offset      = false;
end

% If an affine offset has been specified, integrate those offsets into
% each smooth function, then remove that portion of the array.
if n_aff == n_proj + 1,
    for k = 1 : n_smooth,
        offX = affineF{k,end};
        if isempty(offX),
            continue;
        elseif isa(offX,'function_handle'),
            if adjoint, pos = 'column'; else pos = 'row'; end
            error( 'The elements in the last %s must be constant matrices.', pos );
        elseif isnumeric(offX) && numel(offX) == 1 && offX == 0,
            continue;
        else
            otp_dims{k} = size(offX);
            if nnz(offX),
                smoothF{k} = @(x)smoothF{k}( x + offX );
            end
        end
    end
    n_aff = n_aff - 1;
    affineF(:,end) = [];
end

% -- Todo:
%   If opts.debug = true, then before calling linop_stack,
%   we should print a message showing the sizes of everything
%   (this is useful when the linear portion is entered as "1" or "[]"
%    and we have automatically determined the size ).

% Transpose back, if necessary; then check dimensions
if adjoint,
    affineF = affineF';
    [ linearF, otp_dims, inp_dims ] = linop_stack( affineF, otp_dims, inp_dims, debug );
    linearF = linop_adjoint( linearF );
else
    [ linearF, inp_dims, otp_dims ] = linop_stack( affineF, [], [], debug );
end
identity_linop = isequal( linearF, @linop_identity ); % doesn't always catch identities
square_linop = identity_linop || isequal( inp_dims, otp_dims );
adj_arr = [0,2,1];
if countOps,
    apply_linear = @(x,mode)solver_apply( 3, linearF, x, mode );
else
    apply_linear = linearF;
end

%
% Smooth function, pass 2: integrate scale, counting
%

smoothF = smooth_stack( smoothF );
if maxmin < 0,
    smoothF = tfunc_scale( smoothF, -1 );
end
if countOps,
    apply_smooth = @(x)solver_apply( 1:(1+(nargoutt>1)), smoothF, x );
else
    apply_smooth = smoothF;
end

%
% Projector, pass 2: supply default, stack it up, etc.
%

if isempty( projectorF ),
    n_proj = 0;
    projectorF = proj_Rn;
else
    projectorF = prox_stack( projectorF );
end
if countOps,
    apply_projector = @(varargin)solver_apply( 4:(4+(nargoutt>1)), projectorF, varargin{:} );
else
    apply_projector = projectorF;
end

%
% Initialize the op counts
%

if countOps,
    global tfocs_count___
    tfocs_count___ = [0,0,0,0,0];
    maxCounts = maxCounts(:)';
end

%
% Construct the common initial values
%

L = L0;
theta = Inf;
f_v_old = Inf;
x = []; A_x = []; f_x = Inf; C_x = Inf; g_x = []; g_Ax = [];
restart_iter = 0;
warning_lipschitz = 0;
backtrack_simple = true;
backtrack_tol = 1e-10;
backtrack_steps = 0;

%
% Try to determine the size of the input, and construct the initial point.
%

% Attempt 1: From x0 itself
zero_x0 = true;
if ~isempty( x0 ),
    if isa( x0, 'tfocs_tuple' )
        x0  = cell(x0);
    end
    if isa( x0, 'cell' ),
        n_x0 = numel( x0 );
        if n_x0 == 1,
            x0 = x0{1};
        else
            x0 = tfocs_tuple( x0 );
        end
    else
        n_x0 = 1;
    end
    if n_proj && n_proj ~= n_x0,
        error( 'Size conflict detected between the projector and x0.' );
    end
    zero_x0 = ~nnz(x0);
% Attempt 2: From the linear operator dimensions
elseif ~size_ambig( inp_dims ),
    x0 = tfocs_zeros( inp_dims );
elseif ~size_ambig( otp_dims ),
    A_x = tfocs_zeros( otp_dims );
    x0 = apply_linear( A_x, 2 );
end
if isempty( x0 ),
    error( 'Could not determine the dimensions of the problem. Please supply an explicit value for x0.' );
end
x = x0;
if isinf( C_x ),
    C_x = apply_projector( x );
    if isinf( C_x ),
        zero_x0 = false;
        [ C_x, x ] = apply_projector( x, 1 );
    end
end
if isempty( A_x ),
    if identity_linop || zero_x0 && square_linop,
        A_x = x;
    elseif ~zero_x0 || size_ambig( otp_dims ),  % Jan 2012: todo: give size_ambig the 'offset' information
        A_x = apply_linear( x, 1 ); % celldisp( size(A_x) )
    else
        A_x = tfocs_zeros( otp_dims );
    end
end

% New, Jan 2012
if debug
    if ~adjoint, str1 = 'primal'; str2 = 'dual';
    else, str1 = 'dual'; str2 = 'primal';
    end
    fprintf(fid,'------- DEBUG INFO -----------\n');
    offsetPossible  = offset && ~adjoint;
    fprintf(fid,'Size of %s variable, via method 1\n', str1);
    print_cell_size( size(x0), fid, offsetPossible);
    fprintf(fid,'Size of %s variable, via method 2\n', str1);
    print_cell_size( inp_dims, fid, offsetPossible);

    offsetPossible  = offset && adjoint;
    fprintf(fid,'Size of %s variable, via method 1\n', str2);
    print_cell_size( size(A_x), fid, offsetPossible);
    fprintf(fid,'Size of %s variable, via method 2\n', str2);
    print_cell_size( otp_dims, fid, offsetPossible);
    fprintf(fid,'------------------------------\n');
end
% Added Dec 2012, check for bad sizes that are not integers
if ~isequal( otp_dims, tfocs_round(otp_dims) )
    error('Output dimensions must be integer valued');
end
if ~isequal( inp_dims, tfocs_round(inp_dims) )
    error('Input dimensions must be integer valued');
end

% Final size check
[ isOK1, inp_dims ] = size_compat( size(x0), inp_dims );
[ isOK2, otp_dims ] = size_compat( size(A_x), otp_dims );
if ~isOK1 || ~isOK2,
    if debug
        if ~isOK1
            fprintf(fid,'Debug message: size of %s variables did not line up\n',str1);
        end
        if ~isOK2
            fprintf(fid,'Debug message: size of %s variables did not line up\n',str2);
        end
    end
    error( 'Could not determine the dimensions of the problem. Please supply an explicit value for x0.' );
end
[ f_x, g_Ax ] = apply_smooth( A_x );
if isinf( f_x ),
    error( 'The initial point lies outside of the domain of the smooth function.' );
end
% Adding Feb 2011:
if isa(g_Ax,'tfocs_tuple')
    get_dual = @(g_Ax) get( g_Ax, saddle_ndxs );
else
    get_dual = @(g_Ax) g_Ax;
end

% Theta advancement function
% if mu > 0 && Lexact > mu
if mu > 0 && ~isinf(Lexact) && Lexact > mu, % fixed Dec 2011
    if ~strcmp(alg,'N83')
        warnState = warning('query','backtrace');
        warning off backtrace
        warning('TFOCS:OptionsStructure',' Lexact and mu>0 specifications only give guaranteed convergence with N83 algorithm');
        if strcmp(alg,'AT')
            warning('TFOCS:OptionsStructure',' (With AT algorithm, known to give wrong solutions sometimes when mu>0)');
        end
        warning(warnState);
    end
    % Note that we have not yet derived theory to adapt this to
    %   a local Lipschitz constant but that it should be possible
    theta_scale = sqrt(mu / Lexact);
    theta_scale = ( 1 - theta_scale ) / ( 1 + theta_scale );
    advance_theta = @(theta_old,L,L_old) min(1,theta_old*theta_scale);
else
    advance_theta = @(theta_old,L,L_old) 2/(1+sqrt(1+4*(L/L_old)/theta_old.^2));
end

% Preallocate the arrays for the output structure
out.alg = alg;
out.algorithm = algorithm;
if ~isempty(errFcn) && ~iscell(errFcn),
   errFcn = { errFcn };
end
if ~isempty(stopFcn) && ~iscell(stopFcn),
   stopFcn = { stopFcn };
end
errs = zeros(1,length(errFcn));
if nargoutt == 1,
    saveHist = false;
end
if saveHist,
    [ out.f, out.normGrad, out.stepsize, out.theta ] = deal( zeros(0,1) );
    if countOps,
        out.counts = zeros(0,length(tfocs_count___));
    end
    if ~isempty(errFcn),
        out.err = zeros(0,length(errs));
    end
end
if saddle,
    out.dual = [];
end
n_iter = 0;
status = '';

% Initialize the iterate values
y    = x;    z    = x;
A_y  = A_x;  A_z  = A_x;
C_y  = Inf;  C_z  = C_x;
f_y  = f_x;  f_z  = f_x;
g_y  = g_x;  g_z  = g_x;
g_Ay = g_Ax; g_Az = g_Ax;
norm_x = sqrt( tfocs_normsq( x ) );

% for recomputing linear operators
cntr_Ay     = 0;
cntr_Ax     = 0;

% Special setup for constant step sizes
if beta >= 1,
    beta = 1;
	alpha = 1;
end

% Print the opening text
if fid && printEvery,
	fprintf(fid,'%s\n',algorithm);
	fprintf(fid,'Iter    Objective   |dx|/|x|    step');
    if countOps, fprintf( fid, '       F     G     A     N     P' ); end
    if ~isempty(errFcn)
        nBlanks = max( [0, length(errFcn)*9 - 9] );
        fprintf( fid, '      errors%s',repmat(' ',nBlanks,1) );
    end
    if printStopCrit, fprintf( fid, '    stopping criteria' ); end
    fprintf( fid, '\n' );
	fprintf(fid,'----+----------------------------------' );
    if countOps, fprintf( fid, '+-------------------------------' ); end
%     if ~isempty(errFcn), fprintf( fid, '+-------------------' ); end
    if ~isempty(errFcn)
        fprintf( fid, '+%s', repmat('-',1+length(errFcn)*9,1) );
    end

    if printStopCrit, fprintf( fid, '+%s', repmat('-',19,1) ); end

    fprintf(fid,'\n');
end

% Initialize some variables
just_restarted = false;

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function x = p17_start ( n )

%*****************************************************************************80
%
%% P17_START returns a starting point for optimization for problem 17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables X.
%
%    Output, real X(N), a starting point for the optimization.
%
  x = [ -3.0, -1.0, -3.0, -1.0 ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% -------------------------------------------------------------------------------------------------------------------------
function [init_sz, final_sz] = ideal_size(net, max_sz)
%   Luca Bertinetto, Jack Valmadre, Joao F. Henriques, 2016
% -------------------------------------------------------------------------------------------------------------------------
final_sz = forward(net, max_sz);
init_sz = backward(net, final_sz);
while ~all(init_sz <= max_sz)
    final_sz = final_sz - 1;
    init_sz = backward(net, final_sz);
end

end

function n = forward(net, n)
    for i = 1:numel(net.layers)
        l = net.layers{i};
        switch l.type
        case 'conv'
            m = [size(l.weights{1}, 1), size(l.weights{1}, 2)];
            n = filter(n, l.pad, m, l.stride);
        case 'pool'
            n = filter(n, l.pad, l.pool, l.stride);
        end
    end
end

function n = backward(net, n)
    for i = numel(net.layers):-1:1
        l = net.layers{i};
        switch l.type
        case 'conv'
            m = [size(l.weights{1}, 1), size(l.weights{1}, 2)];
            n = unfilter(n, l.pad, m, l.stride);
        case 'pool'
            n = unfilter(n, l.pad, l.pool, l.stride);
        end
    end
end

function n = filter(n, pad, m, k)
assert(numel(pad) == 1);
n = floor((n + 2*pad - m) / k) + 1;
end

function n = unfilter(n, pad, m, k)
assert(numel(pad) == 1);
n = k*(n - 1) + m - 2*pad;
end

author: bertinetto

repo: siamese-fc


text:
function layer_vec = get_bit_layer(N)
layer_vec = zeros(N, 1);
for phi = 0 : N - 1
    psi = floor(phi/2);
    layer = 0;
    while(mod(psi, 2) == 1)
        psi = floor(psi/2);
        layer = layer + 1;
    end
    layer_vec(phi + 1) = layer;
end
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
function [matching, objective] = greedy_matching_rowwise(iou_matrix)
  assert(size(iou_matrix, 1) <= size(iou_matrix, 2));
  n = size(iou_matrix, 1);
  matching = zeros(n, 1);
  objective = 0;
  for i = 1:n
    % find max element int matrix
 % [max_per_row, max_col_per_row] = max(iou_matrix, [], 2);
%       size(max_per_row)
%       size( max_col_per_row)
   [max_per_row, max_col_per_row] = max(iou_matrix')
    [max_iou,row] = max(max_per_row);
    if max_iou == -inf
      break
    end

    objective = objective + max_iou;
    col = max_col_per_row(row);
    matching(row) = col;
    iou_matrix(row,:) = -inf;
    iou_matrix(:,col) = -inf;
  end
end

author: Cloud-CV

repo: object-proposals


text:
function p_homo = homogenize(p,y);
%HOMOGENIZE Homogenize polynomial
%
% f = homogenize(p,t)

deg   = degree(p);
deg_y = degree(y);
if rem(deg,deg_y)~=0
    error('The degree of the homogenizer is not an even fraction of deg(p).');
end

if 0
    error('The homogenizer must be homogenious.');
end

p_variables = getvariables(p);
p_homo = getbasematrix(p,0)*y^(deg/deg_y);
for i = 1:length(p_variables);
    monom = recover(p_variables(i));
    if degree(monom)<deg
        power = (deg-(degree(monom)))/deg_y;
        p_homo = p_homo + getbasematrix(p,p_variables(i))*monom*y^power;
    else
        p_homo = p_homo + getbasematrix(p,p_variables(i))*monom;
    end;
end
% Reset info about conic terms
p_homo.conicinfo = [0 0];

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function J = computeCost(X, y, theta)
%COMPUTECOST Compute cost for linear regression
%   J = COMPUTECOST(X, y, theta) computes the cost of using theta as the
%   parameter for linear regression to fit the data points in X and y

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta
%               You should set J to the cost.

## for i = 1:m,
##   J = J + ((X(i, :) * theta) - y(i)) ^ 2;
## end;

J = sum(((X * theta) - y) .^ 2);

J = 1 / (2 * m) * J;

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

end

author: fewtime

repo: ML

save_path: github-repos/MATLAB/fewtime-ML


text:
function [um] = pm2um(pm)
% Convert length from picometers to micrometers (or microns).
% Chad A. Greene 2012
um = pm*1e-6;

author: Sable

repo: mcbench-benchmarks


text:
function jed = yjf_to_jed_gregorian ( y, j, f )

%*****************************************************************************80
%
%% YJF_TO_JED_GREGORIAN converts a Gregorian YJF date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, J, real F, the YJF date.
%
%    Output, real JED, the Julian Ephemeris Date.
%

%
%  Copy the input.
%
  y1 = y;
  j1 = j;
  f1 = f;
%
%  Check the input.
%
  [ y1, j1, ierror ] = yj_check_gregorian ( y1, j1 );

  if ( ierror ~= 0 )
    jed = -1.0;
    return
  end
%
%  Convert the input.
%
  [ y2, m2, d2, f2 ] = yjf_to_ymdf_gregorian ( y1, j1, f1 );

  jed = ymdf_to_jed_gregorian ( y2, m2, d2, f2 );

  return
end


text:
   function [mean_pos,ma_pos] = mean_posit_track(t,s,tau,train_len)
%| function [mean_pos,ma_pos] = mean_posit_track(t,s,tau,train_len)
%| INPUT
%|	t: time indices
%|	s: position sequence (should be same length as t)
%|	tau: delay parameter,
%|		set to be the equivalence index to about 0.5sec for breathing.
%|	train_Len: training length (starting up stage),
%|		approximately one period's equivalence
%|
%| OUTPUT
%|	mean_pos: a N*2 matrix,
%|		for practical purposes, take the first column mean_pos(:,1)
%|	ma_pos: naive moving average.
%|
%| Dan Ruan, University of Michigan, 2007
%|
%| Based on the paper by Dan Ruan, J A Fessler, James M Balter
%| Mean position tracking of respiratory motion
%| Med. Phys. 35(2):782-92, Feb. 2008
%| doi 10.1118/1.2825616

if nargin < 1, help(mfilename), error(mfilename), end
if streq(t, 'test'), self_test, return, end

s_complete = s;
s_d = s(1:end-tau);
s = s(tau+1:end);
for tt = 1:length(t)-train_len-tau;
        s_temp = s(tt:tt+train_len);
        sd_temp = s_d(tt:tt+train_len);
        params = fitellipse(s_temp,sd_temp);
        mean_pos(tt,:) =  [params(1),params(2)];
        ma_pos(tt) = mean(s(tt:tt+train_len));
end

plot(t,s_complete,'linewidth',2);
hold on
plot(t(1+tau+train_len:end),ma_pos,'g:','linewidth',2);
plot(t(1+tau+train_len:end),mean_pos(:,1),'r--','linewidth',2);
hold off
set(gca,'DataAspectRatio',[5 5 1]);
set(gca,'fontsize',14,'FontWeight','demi');
h = axis;
axis(h + [0 0 -1 1]);
xlabel('time (second)');
ylabel('displacement');
legend('Observation','Estimation with MA','Estimate with Ellipse Center','Location','SouthWest');

%
function a = fitellipse(X,Y)
% fit ellipse with sample locations in X,Y
% normalize data
mx = mean(X);
my = mean(Y);
sx = (max(X)-min(X))/2;
sy = (max(Y)-min(Y))/2;

x = (X-mx)/sx;
y = (Y-my)/sy;

% Force to column vectors
x = x(:);
y = y(:);

% Build design matrix
D = [ x.*x  x.*y  y.*y  x  y  ones(size(x)) ];

% Build scatter matrix
S = D'*D;

% Build 6x6 constraint matrix
C(6,6) = 0; C(1,3) = -2; C(2,2) = 1; C(3,1) = -2;

  % Break into blocks
  tmpA = S(1:3,1:3);
  tmpB = S(1:3,4:6);
  tmpC = S(4:6,4:6);
  tmpD = C(1:3,1:3);
  tmpE = inv(tmpC)*tmpB';
  [evec_x, eval_x] = eig(inv(tmpD) * (tmpA - tmpB*tmpE));

  % Find the positive (as det(tmpD) < 0) eigenvalue
  I = find(real(diag(eval_x)) < 1e-8 & ~isinf(diag(eval_x)));

  % Extract eigenvector corresponding to negative eigenvalue
  A = real(evec_x(:,I));

  % Recover the bottom half...
  evec_y = -tmpE * A;
  A = [A; evec_y];

% unnormalize
par = [
  A(1)*sy*sy,   ...
      A(2)*sx*sy,   ...
      A(3)*sx*sx,   ...
      -2*A(1)*sy*sy*mx - A(2)*sx*sy*my + A(4)*sx*sy*sy,   ...
      -A(2)*sx*sy*mx - 2*A(3)*sx*sx*my + A(5)*sx*sx*sy,   ...
      A(1)*sy*sy*mx*mx + A(2)*sx*sy*mx*my + A(3)*sx*sx*my*my   ...
      - A(4)*sx*sy*sy*mx - A(5)*sx*sx*sy*my   ...
      + A(6)*sx*sx*sy*sy   ...
      ]';

% Convert to geometric radii, and centers

thetarad = 0.5*atan2(par(2),par(1) - par(3));
cost = cos(thetarad);
sint = sin(thetarad);
sin_squared = sint.*sint;
cos_squared = cost.*cost;
cos_sin = sint .* cost;

Ao = par(6);
Au =   par(4) .* cost + par(5) .* sint;
Av = - par(4) .* sint + par(5) .* cost;
Auu = par(1) .* cos_squared + par(3) .* sin_squared + par(2) .* cos_sin;
Avv = par(1) .* sin_squared + par(3) .* cos_squared - par(2) .* cos_sin;

% ROTATED = [Ao Au Av Auu Avv]

tuCentre = - Au./(2.*Auu);
tvCentre = - Av./(2.*Avv);
wCentre = Ao - Auu.*tuCentre.*tuCentre - Avv.*tvCentre.*tvCentre;

uCentre = tuCentre .* cost - tvCentre .* sint;
vCentre = tuCentre .* sint + tvCentre .* cost;

Ru = -wCentre./Auu;
Rv = -wCentre./Avv;

Ru = sqrt(abs(Ru)).*sign(Ru);
Rv = sqrt(abs(Rv)).*sign(Rv);

a = [uCentre, vCentre, Ru, Rv, thetarad];

% self_test
function self_test
t = linspace(0, 60, 601);
dt = t(2) - t(1);
s = (8 + cos(2*pi/6*t)).^2 - 8^2;
tau = round(0.5 / dt);
train_len = round(5.1 / dt);

[mean_pos ma_pos] = mean_posit_track(t, s, tau, train_len);

prompt
tt = t(tau+train_len+1:end);
clf, plot(t, s, '-', tt, mean_pos(:,1), '--', tt, ma_pos, ':', 'linewidth', 2)
legend('signal', 'proposed', 'moving avg')
xlabel t

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function spm_normalise_disp(matname,VF)
% Display results of spatial normalisation
% FORMAT spm_normalise_disp(matname)
% matname - name of sn3d.mat file
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id$

fg = spm_figure('FindWin','Graphics');
if isempty(fg), return; end;

if nargin<1, matname = spm_select(1,'.*_sn.mat$','Select parameter file'); end;

if ischar(matname),
    t = load(deblank(matname));
else %assume it is a structure
    t = matname;
end;

if nargin<2, VF = t.VF(1); end;

Q = t.VG(1).mat*inv(t.Affine)/VF.mat;

spm_figure('Clear','Graphics');
ax=axes('Position',[0.1 0.51 0.8 0.45],'Visible','off','Parent',fg);
text(0,0.90, 'Spatial Normalisation','FontSize',16,'FontWeight','Bold',...
    'Interpreter','none','Parent',ax);
text(0,0.75, [ 'Image     : ' VF.fname],'FontSize',12,'FontWeight','Bold',...
    'Interpreter','none','Parent',ax);
%text(0,0.7, [ 'Parameters : ' spm_str_manip(matname,'sd')],'FontSize',12,...
%   'Interpreter','none','Parent',ax);

%str = 'no flipping';
%if det(t.Affine(1:3,1:3))<0, str = 'image flipped'; end;
text(0,0.6, 'Linear {affine} component','FontWeight','Bold',...
    'Interpreter','none','Parent',ax);
text(0,0.55, sprintf('X1 = %0.3f*X %+0.3f*Y %+0.3f*Z %+0.3f',Q(1,:)),...
    'Interpreter','none','Parent',ax);
text(0,0.50, sprintf('Y1 = %0.3f*X %+0.3f*Y %+0.3f*Z %+0.3f',Q(2,:)),...
    'Interpreter','none','Parent',ax);
text(0,0.45, sprintf('Z1 = %0.3f*X %+0.3f*Y %+0.3f*Z %+0.3f',Q(3,:)),...
    'Interpreter','none','Parent',ax);

d = [size(t.Tr) 1 1 1];
d = d(1:3);

if prod(d)>1 && isfinite(t.flags.reg),
    text(0,0.35, sprintf('%d nonlinear iterations',t.flags.nits),...
        'Interpreter','none','Parent',ax);
    text(0,0.30, sprintf('%d x %d x %d basis functions',d),...
        'Interpreter','none','Parent',ax);
else
    text(0,0.35, 'No nonlinear components',...
        'Interpreter','none','Parent',ax);
end;

spm_orthviews('Reset');
spm_orthviews('Image',t.VG(1).fname,[0.01 0.1 .48 .6]);
VN = spm_write_sn(VF,matname);
h2 = spm_orthviews('Image',VN,[.51 0.1 .48 .6]);
spm_orthviews('Space',h2);
spm_print;
drawnow;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function value = index2_row ( i_min, i, i_max, j_min, j, j_max, index_min )

%*****************************************************************************80
%
%% INDEX2_ROW indexes a 2D array by row.
%
%  Discussion:
%
%    Entries of the array are indexed starting at entry (I_MIN,J_MIN),
%    and increasing the column index first.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I_MIN, I, I_MAX, for row indices,
%    the minimum, the index, and the maximum.
%
%    Input, integer J_MIN, J, J_MAX, for column indices,
%    the minimum, the index, and the maximum.
%
%    Input, integer INDEX_MIN, the index of element (I_MIN,J_MIN).
%    Typically, this is 0 or 1.
%
%    Output, integer VALUE, the index of element (I,J).
%
  value = index_min + ( j - j_min ) + ( i - i_min ) * ( j_max + 1 - j_min );

  return
end


text:
classdef PoolingNode < GraphNode
    properties
        poolingDimension = 3;      % poolinig happens along this dimension
        poolingType = 'max';       % [max|min|mean|median|sum]
        selectedIdx = [];          % remember the index of data selected by pooling in [max|min|median]
    end

    methods
        function obj = PoolingNode(dimOut, poolingType, poolingDimension)
            obj = obj@GraphNode('Pooling',dimOut);
            if nargin>=2
                obj.poolingType = poolingType;
            end
            if nargin>=3
                obj.poolingDimension = poolingDimension;
            end
        end

        function obj = forward(obj,prev_layers)
            obj = obj.preprocessingForward(prev_layers);
            input = prev_layers{1}.a;

            switch lower(obj.poolingType)
                case 'max'
                    [obj.a, obj.selectedIdx] = max(input, [], obj.poolingDimension);
                case 'min'
                    [obj.a, obj.selectedIdx] = min(input, [], obj.poolingDimension);
                case 'mean'
                    obj.a = mean(input, obj.poolingDimension);
                case 'median'
                    obj.a = median(input, obj.poolingDimension);
                    % need to implement selectedIdx by myself in the future
                case 'sum'
                    obj.a = sum(input, obj.poolingDimension);
            end

            obj = forward@GraphNode(obj, prev_layers);
        end

        function obj = backward(obj,prev_layers, future_layers)
            if obj.skipGrad || obj.skipBP
                return;
            end

            future_grad = obj.GetFutureGrad(future_layers);
            input = prev_layers{1}.a;
            [D(1), D(2), D(3), D(4)] = size(input);

            switch lower(obj.poolingType)
                case {'max', 'min', 'median'}
                    obj.grad{1} = obj.AllocateMemoryLike(D, future_grad);
                    obj.grad{1}(obj.selectedIdx) = future_grad;
                case {'mean','sum'}
                    if strcmpi(obj.poolingType, 'mean')
                        future_grad = future_grad / D(obj.poolingDimension);
                    end
                    shape = ones(1,4);
                    shape(obj.poolingDimension) = D(obj.poolingDimension);
                    obj.grad{1} = repmat(future_grad, shape);
            end

            obj = backward@GraphNode(obj, prev_layers, future_layers);
        end

    end

end

author: singaxiong

repo: SignalGraph


text:
clear all;

f = @(p)[p(1) p(1)^2]';

% sampling intervals for each parameter
domain = [-2 2];
% grid size: number of grid points for each parameter
gridsize = 19;

% sampling
D = sampling_vec(f, domain, gridsize);

% hosvd
[S U] = hosvd(D, [1 0]);
U = U{1};

[W V] = genhull(U,'box');

Uhat = -ones(size(W));
Uhat(1,1) = 0;
Uhat(8,1) = 0.4;
Uhat(10,1) = 0.6;
Uhat(19,1) = 0;
Uhat(5,2) = 0.6;
Uhat(1,3) = 0;
Uhat(19,3) = 1;
Uhat(1,4) = 1;
Uhat(19,4) = 0;
M = hull_manip(W,Uhat);
Wmod = W*M;
Vmod = M\V;
[i,j] = find(Uhat > -1);
uhat = Uhat(Uhat > -1);
n = length(i);
%for k=1:n
%    [uhat(k) Wmod(i(k),j(k))] %#ok<NOPTS>
%end

close all
set(0, 'DefaultFigureWindowStyle', 'docked')
plothull(U);
plothull(W);
plothull(Wmod);
hold on
plot(i,uhat,'kx')
hold off

figure
hold on
plot(D(:,1), D(:,2))
plot(V*S(:,1), V*S(:,2), 'go')
plot(Vmod*S(:,1), Vmod*S(:,2), 'mo')
hold off


text:
function stroud_test03 ( )

%*****************************************************************************80
%
%% TEST03 tests BALL_UNIT_07_3D, BALL_UNIT_14_3D, BALL_UNIT_15_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  global FUNC_3D_INDEX;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03\n' );
  fprintf ( 1, '  For integrals in the unit ball in 3D:\n' );
  fprintf ( 1, '  BALL_UNIT_07_3D uses a formula of degree 7;\n' );
  fprintf ( 1, '  BALL_UNIT_14_3D uses a formula of degree 14;\n' );
  fprintf ( 1, '  BALL_UNIT_15_3D uses a formula of degree 15.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Unit ball volume = %f\n', ball_unit_volume_nd ( 3 ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Rule:      #7             #14           #15\n' );
  fprintf ( 1, '    F(X)\n' );
  fprintf ( 1, '\n' );

  num = function_3d_num ( );

  for i = 1 : num

    FUNC_3D_INDEX = i;

    result1 = ball_unit_07_3d ( 'function_3d' );
    result2 = ball_unit_14_3d ( 'function_3d' );
    result3 = ball_unit_15_3d ( 'function_3d' );

    fname = function_3d_name ( i );

    fprintf ( 1, '  %7s  %12f  %12f  %12f\n', ...
      fname, result1, result2, result3 );

  end

  return
end


text:
function jd=mjd2jd(mjd)
% MJD2JD  Converts Modified Julian Date to Julian Date.
%   Non-vectorized version. See also CAL2JD, DOY2JD, GPS2JD,
%   JD2CAL, JD2DOW, JD2DOY, JD2GPS, JD2MJD, JD2YR, YR2JD.
% Version: 2010-03-25
% Usage:   jd=mjd2jd(mjd)
% Input:   mjd - Modified Julian date
% Output:  jd  - Julian date

% Copyright (c) 2011, Michael R. Craymer
% All rights reserved.
% Email: mike@craymer.com

if nargin ~= 1
  warning('Incorrect number of arguments');
  return;
end

jd=mjd+2400000.5;


text:
function pred=lengauer_tarjan_dominator_tree(A,u,varargin)
% LENGAUER_TARJAN_DOMINATOR_TREE Compute a dominator tree for a graph.
%
% A dominator tree encodes dominates relations.  A vertex u dominates a
% vertex v if all paths to v must go through u.  In an undirected graph,
% this means that if u dominates v, then v will be disconnected from the
% graph if we remove u and start the search at the root.  A dominator tree
% is rooted at a particular vertex and the entire graph must be reachable
% from that vertex.  (Historically, these were called flowgraphs.)
%
% p = lengauer_tarjan_dominator_tree(A,u) returns the predecessor array for
% the dominator tree rooted at u.
%
% The runtime is O((V+E) log (V+E)) and the algorithm works on unweighted,
% directed graphs.
%
% ... = lengauer_tarjan_dominator_tree(A,...) takes a set of
% key-value pairs or an options structure.  See set_matlab_bgl_options
% for the standard options.
%   There are no additional options for this function.
%
% Note: this function does not depend upon the non-zero values of A, but
% only uses the non-zero structure of A.
%
% Example:
%   load('graphs/dominator_tree_example.mat');
%   p=lengauer_tarjan_dominator_tree(A,1);

% David Gleich
% Copyright, Stanford University, 2007-2008

%% History
%  2007-07-13: Initial version
%%

[trans check full2sparse] = get_matlab_bgl_options(varargin{:});
if full2sparse && ~issparse(A), A = sparse(A); end
if check, check_matlab_bgl(A,struct()); end
if trans, A = A'; end

pred = dominator_tree_mex(A,u);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function z = uminus( x )

%   Disciplined convex programming information for UMINUS (-):
%      Unary minus may be used in DCPs without restrictions---with the
%      understanding, of course, that it produces a result with the
%      opposite cuvature to its argument; i.e., the negative of a convex
%      expression is concave, and vice versa.
%
%   Disciplined geometric programming information for UMINUS(-):
%      Negation of non-constant values may not be used in disciplined
%      geometric programs.

z = cvx( x.size_, -x.basis_ );
tt = cvx_isvalid( z, true );
if ~all( tt ),
    cvx_dcp_error( '-', 'unary', cvx_subsref( x, ~tt ) );
end

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

author: yu-jiang

repo: radpbook

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


text:
function [Ri, RiFull] = syntheticCamera(F, type)

if nargin < 2
    type = 'continuous';
end

Ri = cell(F, 1);
RiFull = cell(F, 1);

switch type
    case 'random'
        for i = 1:F
            R = orth(randn(3));
            Ri{i} = R(2:3, :);
            RiFull{i} = R;
        end

    case 'continuous'
        Ux = @(u) [   0, -u(3),  u(2); ...
                   u(3),     0, -u(1); ...
                  -u(2),  u(1),    0];
        theta = randn(1); % create an angle randomly
        u = randn(3, 1);
        u = u/norm(u); % create a unit axis randomly
        for i = 1:F
            theta = theta + 1;
            % create a rotation matrix
            orthR = cos(theta)*eye(3) + sin(theta)*Ux(u) + (1-cos(theta))*kron(u,u');
            Ri{i} = orthR(2:3, :);
            RiFull{i} = orthR(1:3, :);
        end

    case 'vertical'
        Ux = @(u) [   0, -u(3),  u(2); ...
                   u(3),     0, -u(1); ...
                  -u(2),  u(1),    0];
        %theta = randn(1);
        theta = 360/F;
        u = [0.2;0;1]; % z-axis rotation
        for i = 1:F
            % in degrees
            orthR = cosd(theta*i)*eye(3) + sind(theta*i)*Ux(u) + (1-cosd(theta*i))*kron(u,u');

            Ri{i} = orthR(2:3, :);
            RiFull{i} = orthR(1:3, :);
        end

    case 'real'
        Ux = @(u) [   0, -u(3),  u(2); ...
                   u(3),     0, -u(1); ...
                  -u(2),  u(1),    0];
        %theta = randn(1);
        theta = 2;
        totalRotDeg = 120;
        u = [0;0;1]; % z-axis rotation
        for i = 1:F

            orthR = cosd(theta)*eye(3) + sind(theta)*Ux(u) + (1-cosd(theta))*kron(u,u');

            if i >= round(F/2) && i <= round((3*F)/4)
                camActNum = round(F/2) - round((3*F)/4);
                degrees = totalRotDeg / camActNum;

                % in degrees
                theta = theta + degrees;
                orthR = cosd(theta)*eye(3) + sind(theta)*Ux(u) + (1-cosd(theta))*kron(u,u');

                Ri{i} = orthR(2:3, :);
                RiFull{i} = orthR(1:3, :);
            else
                Ri{i} = orthR(2:3, :);
                RiFull{i} = orthR(1:3, :);
            end
        end

        case '360'
        Ux = @(u) [   0, -u(3),  u(2); ...
                   u(3),     0, -u(1); ...
                  -u(2),  u(1),    0];
        %theta = randn(1);
        theta = 360/F;
        u = [0;0;1]; % z-axis rotation
        for i = 1:F
            % in degrees
            orthR = cosd(theta*i)*eye(3) + sind(theta*i)*Ux(u) + (1-cosd(theta*i))*kron(u,u');
            Ri{i} = orthR(2:3, :);
            RiFull{i} = orthR(1:3, :);
        end

    otherwise
        error('%s is a wrong camera type!', type);
end

author: jhonykaesemodel


text:
function cs_install (do_pause)
%CS_INSTALL: compile and install CXSparse for use in MATLAB.
%   Your current working directory must be CXSparse/MATLAB in order to use this
%   function.
%
%   The directories
%
%       CXSparse/MATLAB/CSparse
%       CXSparse/MATLAB/Demo
%       CXSparse/MATLAB/UFget
%
%   are added to your MATLAB path (see the "pathtool" command to add these to
%   your path permanently, for future MATLAB sessions).
%
%   Next, the MATLAB CXSparse demo program, CXSparse/MATLAB/cs_demo is executed.
%   To run the demo with pauses so you can see the results, use cs_install(1).
%   To run the full MATLAB test programs for CXSparse, run testall in the
%   Test directory.
%
%   Example:
%       cs_install          % install and run demo with no pauses
%       cs_install(1)       % install and run demo with pauses
%
%   See also: cs_demo
%
%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

fprintf ('Compiling and installing CXSparse\n') ;
if (nargin < 1)
    do_pause = 0 ;
end

if (do_pause)
    input ('Hit enter to continue: ') ;
end
addpath ([pwd filesep 'CSparse']) ;
addpath ([pwd filesep 'Demo']) ;

v = getversion ;
if (v >= 7.0)
    addpath ([pwd filesep 'UFget']) ;
else
    fprintf ('UFget not installed (MATLAB 7.0 or later required)\n') ;
end

cd ('CSparse') ;
cs_make (1) ;
cd ('../Demo') ;
cs_demo (do_pause)

%-------------------------------------------------------------------------------
function v = getversion
% determine the MATLAB version, and return it as a double.
v = sscanf (version, '%d.%d.%d') ;
v = 10.^(0:-1:-(length(v)-1)) * v ;


text:
function SaveVARtoMat(hfig)
save_dir = getappdata(hfig,'save_dir');

% copy of VAR files will be saved into this subfolder:
arcmatfolder = fullfile(save_dir, 'arc mat');
if ~exist(arcmatfolder, 'dir')
    mkdir(arcmatfolder);
end

global VAR; %#ok<NUSED>
timestamp  = datestr(now,'mmddyy_HHMM');
matname = [timestamp '.mat'];
save(fullfile(arcmatfolder,matname),'VAR','-v6');

% also save the current VAR file
save(fullfile(save_dir,'VAR_new.mat'),'VAR','-v6');
disp('saved both to workspace and .mat');

end

author: xiuyechen

repo: FishExplorer


text:
function f2 = p00_f2 ( problem, x )

%*****************************************************************************80
%
%% P00_F2 evaluates the second derivative for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the values of the variables.
%
%    Output, real F2, the second derivative.
%
  if ( problem == 1 )
    f2 =  p01_f2 ( x );
  elseif ( problem == 2 )
    f2 =  p02_f2 ( x );
  elseif ( problem == 3 )
    f2 =  p03_f2 ( x );
  elseif ( problem == 4 )
    f2 =  p04_f2 ( x );
  elseif ( problem == 5 )
    f2 =  p05_f2 ( x );
  elseif ( problem == 6 )
    f2 =  p06_f2 ( x );
  elseif ( problem == 7 )
    f2 =  p07_f2 ( x );
  elseif ( problem == 8 )
    f2 =  p08_f2 ( x );
  elseif ( problem == 9 )
    f2 =  p09_f2 ( x );
  elseif ( problem == 10 )
    f2 =  p10_f2 ( x );
  elseif ( problem == 11 )
    f2 =  p11_f2 ( x );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_F2 - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem number PROBLEM = %d\n', problem );
    error ( 'P02_F2 - Fatal error!' );
  end

  return
end


text:
R = input ('Clear data? [y/n] \n','s');
if strcmp(R,'y')
    clear
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function show_reg_warp(I,U,V,grid_size)
%SHOW_REG_WARP Summary of this function goes here
%   Detailed explanation goes here
% Input:
%   I - image matrix or handle to the axes
%   (U,V) - translation field v(x) = (U(x),V(x)), s.t. T(x) = x + v(x)
%   grid_size - number of points per dimension

if nargin < 1
    I = imread('pout.tif');
end

if nargin < 3
    [rows,cols,b] = size(I);
    U = zeros(rows,cols);
    V = zeros(rows,cols);
end

if nargin < 4
    grid_size = 20;
end

[rows,cols] = size(U);

x1 = linspace(1, cols, grid_size);
y1 = linspace(1, rows, grid_size);
[X,Y] = meshgrid(x1,y1);
U1 = interp2((1:cols),(1:rows),U,X,Y);
V1 = interp2((1:cols),(1:rows),V,X,Y);
X_target = X - U1; % X_target + U1 = X;
Y_target = Y - V1; % Y_target + V1 = Y;
x = X_target(:);
y = Y_target(:);

[W,M] = image2graph(ones(length(y1),length(x1)),0.05,1e-9);

% [A,J2,A2] = create_adjacency_graph(X,'epsball',1.01,1);

if ishandle(I)
    ax_h = I;
elseif ~isempty(I)
    figure,imshow(I);
    ax_h = gca;
end
hold(ax_h, 'on');

for i = 1:length(x)
    neighbor_inds = find(M(i,:));
    for j = 1:length(neighbor_inds)
        plot(ax_h,[x(i),x(neighbor_inds(j))],[y(i),y(neighbor_inds(j))],'-b',...
            'linewidth', 1);
    end
end

end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function trajectory_planning()
%%
close all
clear
%% requirements
% given obstacle field O, start configuration z_s with centroid s and destination g for the formation's centroid.
import iris.inflate_region.*
import iris.drawing.*
import iris.thirdParty.polytopes.*
% constants
lcx = 0.6; lcy = 0.6; lcz = 0.2; % size of object
lmx = 0.6; lmy = 0.4; lmz = 0.3; lm_b = 0.15; % size of platform
% obstacles
xo1 = [3.75,4.75,0.5,0,0,0];
xo2 = [8.25,3,0.5,0,0,0];
obst(1) = pkgMechanics.RigidCuboid(1,xo1,[1.5,3.5,1]);
obst(2) = pkgMechanics.RigidCuboid(1,xo2,[1.5,1,1]);
% obstacles cell
obstcell = createobstcell(obst);
% range
range.lb = [0;0];
range.ub = [10;10];
%% initialize an empty graph
nodez = []; % z vector in row
nodepoly = []; % polytope struct
nodezinpoly = {}; % z index
% initialize edges
edgez = []; % z index
%% configuration
a0 = 0.7; w0 = 0.5;
config = [a0,w0];
%% generate start and goal
start = [1.5,1.5];     % start position
goal = [8.5,8.5];     % goal position
[nodez, nodepoly, nodezinpoly]= startxgoal(nodez,nodepoly,nodezinpoly,[start;goal],config,obstcell,range,[lcx,lcy,lcz],[lmx,lmy,lmz]);
%% detemine z
z0 = [start,config];
[xr0,xc0,xm0]= z2x(z0,[lcx,lcy,lcz],[lmx,lmy,lmz]);
%% build object
cub = pkgMechanics.RigidCuboid(1,xc0,[lcx,lcy,lcz]);
%% build platforms
for i=1:4
    plat(i) = pkgMechanics.RigidCuboid(1,xm0(i,:),[lmx,lmy,lmz]);
end
%% build robot arms
Tc2e = lc2Tce([lcx,lcy,lcz]); % => constants
Tm2b = transl([lm_b,0,lmz/2]); % => constants
mdl_puma560
for i=1:4
    rob(i) = SerialLink('name','robot');
    copy(rob(i),p560); rob(i).name = ['robot',num2str(i)];
    rob(i).base = x2T(xm0(i,:))*Tm2b;
    q0(i,:) = rob(i).ikine6s(x2T(xc0)*Tc2e{i},'ru');% right hand elbow up
end
%% visualize map
figure;
draw_region_2d(nodepoly,obstcell.vert,range);
drawplat3d(nodez,[lcx,lcy,lcz],[lmx,lmy,lmz])
%% add goal node to graph
[nodez,nodezinpoly,nodegrid,edgez]=creategraph(nodez,nodepoly,nodezinpoly,edgez,1,config,range,[lcx,lcy,lcz],[lmx,lmy,lmz]);
[noderoute,polyroute] = shortestpath(nodez,nodezinpoly,nodegrid,1,2);
%% generate random nodes
counter = 0;
while isempty(noderoute)
    nodenum = 1;
    [nodez,nodepoly,nodezinpoly] = randnode(nodez,nodepoly,nodezinpoly,nodenum,config,obstcell,range,[lcx,lcy,lcz],[lmx,lmy,lmz]);
    % create undirected graph and its edges
    [nodez,nodezinpoly,nodegrid,edgez]=creategraph(nodez,nodepoly,nodezinpoly,edgez,nodenum,config,range,[lcx,lcy,lcz],[lmx,lmy,lmz]);
    [noderoute,polyroute] = shortestpath(nodez,nodezinpoly,nodegrid,1,2);
    counter = counter + 1;
end
draw_region_2d(nodepoly(polyroute),[],range);
drawplat3d(nodez(noderoute,:),[lcx,lcy,lcz],[lmx,lmy,lmz]);
hold on
x1=nodez(noderoute,1);
y1=nodez(noderoute,2);
plot(x1,y1,'r-');
hold off
if ~exist('shortestpath.mat', 'file')
    save('shortestpath.mat','nodez','nodepoly','edgez','noderoute','polyroute');
end
%%%%%%%%% the above is the whole global planning part %%%%%%%%%%%

%%%%%%%%% the following is the animation %%%%%%%%%%
load('shortestpath2.mat','nodez','nodepoly','edgez','noderoute','polyroute');
%% save data
% trajectory
dt = 0.5; tacc = 1;
[zarray, dzarray, tarray] = calctraj(nodez(noderoute,:),dt,tacc);
%% visualize map
figure;
xlim([range.lb(1) range.ub(1)]);
ylim([range.lb(2) range.ub(2)]);
hold on
% % obstacles
% draw_region_2d([],obstcell.vert,range);
% fig.obst = patch('Vertices',obst.verticesStates.position','Faces',obst.faces,'FaceColor','g','FaceAlpha',0.5);
% pause(0.2)
% % object
% fig.cub = patch('Vertices',cub.verticesStates.position','Faces',cub.faces,'FaceColor','m','FaceAlpha',0.5);
% % robot arm and base
zlim([0 8]); plotaxis = gca;
plotopt = {'noname', 'noraise', 'nobase', 'noshadow', 'notiles', 'nowrist', 'nojaxes', 'delay', 0, 'workspace',[plotaxis.XLim,plotaxis.YLim,plotaxis.ZLim]};
% for i=1:4
%     fig.plat(i) = patch('Vertices',plat(i).vertw,'Faces',plat(i).face,'FaceColor','y','FaceAlpha',0.5);
%     rob(i).plot(q0(i,:),plotopt{:});
% end
% pause(0.2)
% % polytope
% for i=1:length(polyroute)
%     draw_region_2d(nodepoly(polyroute(i)),[],range);
%     pause(0.2)
% end
% % graph
% for kk=1:2:size(edgez,1)
%     x1=[nodez(edgez(kk,1),1);nodez(edgez(kk,2),1)];
%     y1=[nodez(edgez(kk,1),2);nodez(edgez(kk,2),2)];
%     line(x1,y1);
% end
% pause(0.2)
% % path
% plot(nodez(noderoute(end),1),nodez(noderoute(end),2),'kd', 'LineWidth', 2)
% plot(nodez(noderoute,1),nodez(noderoute,2),'k-', 'LineWidth', 2)
% pause(0.2)
% hold off
% % obstacles
% cla
% hold on
draw_region_2d([],obstcell.vert,range);
for i=1:length(obst)
    fig.obst = patch('Vertices',obst(i).verticesStates.position','Faces',obst(i).faces,'FaceColor','g','FaceAlpha',0.5);
end
% object
fig.cub = patch('Vertices',cub.verticesStates.position','Faces',cub.faces,'FaceColor','m','FaceAlpha',0.5);
% robot arm and base
for i=1:4
    fig.plat(i) = patch('Vertices',plat(i).verticesStates.position','Faces',plat(i).faces,'FaceColor','b','FaceAlpha',0.5);
    rob(i).plot(q0(i,:),plotopt{:});
end
% path
plot(nodez(noderoute(end),1),nodez(noderoute(end),2),'rd', 'LineWidth', 2)
plot(nodez(noderoute,1),nodez(noderoute,2),'r-', 'LineWidth', 2)
pause(0.2)
view(3)
hold off
%% animation loop
tic;
playspeed = 1;
while toc<tarray(end)/playspeed
    tnow = toc*playspeed;
    % choose via point in time 'tnow'
    [z,dz] = getviatnow(zarray, dzarray, tarray, tnow);
    [xr,xc,xm] = z2x(z,[lcx,lcy,lcz],[lmx,lmy,lmz]);
    % update object
    cub.updateStates(xc',zeros(6,1));
    % set robot arm state
    for i=1:4
        rob(i).base = x2T(xm(i,:))*Tm2b;
        q(i,:) = rob(i).ikine6s(x2T(xc)*Tc2e{i});
        plat(i).updateStates(xm(i,:)',zeros(6,1));
    end
    % animation
    set(fig.cub,'Vertices',cub.verticesStates.position','Faces',cub.faces);
    for i=1:4
        set(fig.plat(i),'Vertices',plat(i).verticesStates.position','Faces',plat(i).faces);
        rob(i).animate(q(i,:));
    end
    drawnow
end
end

function obstacle = createobstcell(obstcub)
n_obs=length(obstcub);
for i=1:n_obs
    obstacle_pts(:,:,i) = obstcub(i).verticesStates.position(1:2,1:4);
end
% reset obstacles vertices in its convexhull (iris requires obstacles to be convex)
obsidx = cell(1,size(obstacle_pts,3));
obstacle.vert = cell(1,size(obstacle_pts,3));
for j = 1:n_obs
    obs = obstacle_pts(:,:,j);
    if size(obs, 2) > 1
        if size(obs, 2) > 2
            obsidx{j} = convhull(obs(1,:), obs(2,:));
        else
            obsidx{j} = [1,2,1];
        end
    end
    obstacle.vert{j} = obs(:,obsidx{j})';
    obstacle.calc{j} = obs(:,obsidx{j});
end
% obstacles cell
obstacle.poly = obstacle.vert{1}';
for i=2:size(obstacle.vert,1)
    obstacle.poly = [obstacle.poly,NaN(2,1),obstacle.vert{i}'];
end
end

function [nodez, nodepoly, nodezinpoly]= startxgoal(nodez,nodepoly,nodezinpoly,randnode,nodeconfig,obstacle,range,lc,lm)
% set nodez to all zero
znum = size(nodez,1);
polynum = length(nodepoly);
% generate nodes outside obstacles
randnum = size(randnode,1);
for i=1:randnum
    % generate ploytope
    randpoly = polytope(obstacle.calc,randnode(i,:)',range);
    % check whether formation exists in a polytope
    [randz,fg] = formationP2z(randpoly,[randnode(i,:),nodeconfig],range,lc,lm);
    if fg==1
        % add this location to nodelocation list
        nodez = [nodez;randz];
        nodepoly = [nodepoly,randpoly];
        nodezinpoly{polynum+i} = znum+i;
        %%%%%
%         hold on
%         plot(randlocation(1),randlocation(2),'go');
%         plot(randz(1),randz(2),'r*');
%         draw_region_2d(randpoly,[],range);
%         drawplat3d(randz);
%         hold off
        %%%%%
    end
end
end

function [nodez, nodepoly, nodezinpoly]= randnode(nodez,nodepoly,nodezinpoly,nodenum,nodeconfig,obstacle,range,lc,lm)
% set nodez to all zero
znum = size(nodez,1);
polynum = length(nodepoly);
% generate nodes outside obstacles
ctr=1; % counter
while (ctr<=nodenum)
    % generate random two number in range of map's border
    randlocation = [rand* (range.ub(1)-range.lb(1)) + range.lb(1);
        rand* (range.ub(2)-range.lb(2)) + range.lb(2)];
    % if this node is not inside any obstacle
    if ~inpolygon(randlocation(1),randlocation(2),obstacle.poly(1,:),obstacle.poly(2,:))
        % chech if inside a existing polytope
        if ~innodepoly(randlocation,nodepoly)
            % generate ploytope
            randpoly = polytope(obstacle.calc,randlocation,range);
            % check whether formation exists in a polytope
            [randz,fg] = formationP2z(randpoly,[randlocation',nodeconfig],range,lc,lm);
            if fg==1
                % add this location to nodelocation list
                nodez = [nodez;randz];
                nodepoly = [nodepoly,randpoly];
                nodezinpoly{polynum+ctr} = znum+ctr;
                %%%%%
%                 hold on
%                 plot(randlocation(1),randlocation(2),'go');
%                 plot(randz(1),randz(2),'r*');
%                 draw_region_2d(randpoly,[],range);
%                 drawplat3d(randz);
%                 hold off
                %%%%%
                ctr=ctr+1;
            end
        end
    end
end
end

function flag = innodepoly(randlocation,nodepoly)
polynum = length(nodepoly);
count = 0;
for i=1:polynum
    if nodepoly(i).A*randlocation-nodepoly(i).b<=0
        count = count+1;
    end
end
if count==0
    flag = false;
else
    flag = true;
end
end

repo: paper-simulation


text:
function [L1,L1_l,L1_x1] = reanchorPlucker(L0,x1)

% REANCHORPLUCKER  Plucker to anchored Plucker line conversion
%   REANCHORPLUCKER(L,X) reanchors the anchored Plucker line L to the point X.
%
%   [L1,L1_l,L1_x] = REANCHORPLUCKER(L,X) returns the Jacobians.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

x0 = L0(1:3);
n0 = L0(4:6);
v0 = L0(7:9);

L1 = [x1; n0 + cross((x0-x1),v0); v0];

if nargout > 1

    Z33 = zeros(3);
    I33 = eye(3);

    L1_l = [...
        Z33      Z33    Z33
        -hat(v0) I33 hat(x0-x1)
        Z33      Z33    I33];

    L1_x1 = [...
        I33
        hat(v0)
        Z33];

end

return

%%
syms t1 t2 t3 x1 x2 x3 a b c d n1 n2 n3 v1 v2 v3 y1 y2 y3 real

L0 = [x1;x2;x3;n1;n2;n3;v1;v2;v3];
y = [y1;y2;y3];

[L1,L1_l,L1_y] = reanchorPlucker(L0,y);

simplify(L1_l - jacobian(L1,L0))
simplify(L1_y - jacobian(L1,y))

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

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

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [y, w] = pts_tps_map( s, t, x, w, FAST, PROGRESS )
% PTS_TPS_MAP  Interpolate/warp/map N-dimensional points using a thin-plate
% spline transformation
%
% [Y, W] = PTS_TPS_MAP(S, T, X, W)
%
%    S is a (P,Ds,N)-volume where each (:,:,i)-matrix has the coordinates
%    of the source points that define the warp.
%
%    T is a (P,Dt,N)-volume where each (:,:,i)-matrix has the coordinates
%    of the target points that define the warp.
%
%      (P is the number of points, Ds and Dt are the dimension and N is the
%      number of configurations).
%
%    X is a (P,Ds,N)-volume where each (:,:,i)-matrix has the coordinates
%    of the points to be interpolated.
%
%    Y is a (P,Dt,N)-volume with the interpolated points.
%
%    X can also be a struct with 3 volumes: X.in (P1,Ds,N), X.endo (P2,Ds,N),
%    X.epi (P3,Ds,N). In this case, 3 warps are computed: With all points
%    for X.in; with the first half of S and T points for X.endo; and with
%    the second half of S and T points for X.epi.
%
%    W is a D2-column matrix where each column has the weight and affine
%    weights vector computed with PTS_TSP_WEIGHTS (if W is empty or
%    missing, then W is computed internally).
%
%    If X is a struct, then W has to be a similar struct with the weights
%    for each of the warps.
%
% ... = PTS_TPS_MAP(S, T, X, W, FAST)
%
%    To speed up computations, all points in the same configuration are
%    warped together using vector and matrix operations.
%
%    This is more memory expensive. If for some reason memory becomes a
%    problem, using FAST=false it is possible to run the algorithm using
%    a "for" loop to warp 1 point at a time. This is roughly 14 times
%    slower. By default the faster method is used, FAST=true.
%
% ... = PTS_TPS_MAP(S, T, X, W, FAST, PROGRESS)
%
%    Making PROGRESS=true, and if you have function STATUSBAR() by
%    Leutenegger Marcel, you will see a progress bar with an estimate of
%    the time remaining to completion if you are using the slow method.
%
%
% See also: pts_tps_weights.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2006-2011 University of Oxford
% Version: 0.6.4
% $Rev: 434 $
% $Date: 2011-06-03 17:16:33 +0100 (Fri, 03 Jun 2011) $
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
error(nargchk(3, 6, nargin));
error(nargoutchk(0, 2, nargout));

% defaults
if (nargin < 5 || isempty(FAST))
    FAST = true;
end
if (nargin < 6 || isempty(PROGRESS))
    PROGRESS = false;
end

% struct and sizes check
P = size(s, 1); % number of control points
if (P ~= size(t, 1))
    error('S and T must have the same number of points')
end

% compute weights if needed
N = size(s, 3); % number of frames
if (nargin < 4 || isempty(w))
    if isstruct(x) % triple warp

        for I = 1:N
            w(I).in = pts_tps_weights(s(:, :, I), t(:, :, I));
            w(I).endo = pts_tps_weights(s(1:(P/2), :, I), ...
                t(1:(P/2), :, I));
            w(I).epi = pts_tps_weights(s((P/2+1):end, :, I), ...
                t((P/2+1):end, :, I));
        end

    else
        w = pts_tps_weights(s, t);
    end
end

% deal with empty point sets
if isempty(x)
    y = [];
    return;
end

if isstruct(x) % triple warp
    if any(~isfield(x, {'endo', 'epi', 'in'}))
        error('Invalid fields for struct in X')
    end

    % compute each warp separately
    y(N) = struct( ...
        'in', zeros(size(x(1).in)), ...
        'endo', zeros(size(x(1).endo)), ...
        'epi', zeros(size(x(1).epi)));
    for I = 1:N
        y(I).in = pts_tps_map(s(:, :, I), t(:, :, I), ...
            x(I).in, w(I).in, FAST);
        y(I).endo = pts_tps_map(s(1:(P/2), :, I), ...
            t(1:(P/2), :, I), x(I).endo, w(I).endo, FAST);
        y(I).epi = pts_tps_map(s((P/2+1):end, :, I), ...
            t((P/2+1):end, :, I), x(I).epi, w(I).epi, FAST);
    end

    return
end

% dimensionality of points
D = size(s, 2);

if (D ~= size(x, 2))
    error(['Dimensionality of S is ' num2str(D) ', of X is ' num2str(size(x, 2))])
end
if (N ~= size(t, 3) || N ~= size(x, 3))
    error('S, T and X must have the same number of point configurations')
end

% init output
y = zeros(size(x, 1), size(w, 2));

PP = size(x, 1); % PP: number of points to be interpolated

if FAST

    % warp every input point
    for I = 1:N % loop frames

        % init aux matrix
        u = zeros(P * PP, D);

        % interleave points to be warped
        for J = 1:D
            u(:, J) = reshape(repmat(x(:, J, I), 1, P)', P * PP, 1);
        end

        % compute norm(Pi - (x,y)).^2
        u = sum((u - repmat(s(:, :, I), PP, 1)) .^ 2, 2);

        % reshape to have 1 row vector per point to be interpolated
        u = reshape(u, P, PP)';

        % thin-plate spline distance function
        % U(r) = r^2 log10(r)
        % note: it's faster to compute the log10 this way than directly
        warning('off', 'MATLAB:log:logOfZero');
        u = 0.5 * u .* log(u) * (1/log(10));
        u(isnan(u)) = 0;
        warning('on', 'MATLAB:log:logOfZero');

        % factor by weights f(x,y) = a1 + ax*x + ay*y + sum(wi*U(|Pi - (x,y)|))
        y(:, :, I) = [u , ones(PP, 1) , x(:, :, I)] * w(:, :, I);

    end

else % loop throug each point; this is slower but less memory consuming

    % show status bar
    if PROGRESS
        delete(statusbar)
        bar = statusbar('Progress bar...');
    else
        bar = [];
    end

    % warp every input point
    for I = 1:N % loop frames
        for J = 1:PP
            % compute norm(Pi - (x,y)).^2
            % U(r) = r^2 log10(r)
            % note: it's faster to compute the log10 this way than directly
            u = sum(...
                (repmat(x(J, :, I), P, 1) - s(:, :, I)) .^ 2, 2)';
            warning('off', 'MATLAB:log:logOfZero');
            u = 0.5 * u .* log(u) * (1/log(10));
            u(isnan(u)) = 0;
            warning('on', 'MATLAB:log:logOfZero');

            % factor by weights f(x,y) = a1 + ax*x + ay*y + sum(wi*U(|Pi - (x,y)|))
            y(J, :, I) = ...
                [u , ...
                1 , ...
                x(J, :, I)] ... % x, y, z, etc
                * w(:, :, I);

            % update progress bar every 10th of the total number of points
            if mod(J, round(PP/10))
                if (PROGRESS && isempty(statusbar(((I-1)*PP+J)/(N*PP), bar)))
                    break
                end
            end

        end

%         % update progress bar every frame
%         if isempty(statusbar(I/N, bar))
%             break
%         end

    end

    % delete progress bar
    delete(bar);

end

author: tomdoel

repo: pulmonarytoolkit


text:
function c=ref_dgt2(f,g1,g2,a1,a2,M1,M2);
%REF_DGT2  Reference DGT2
%
%  Compute a DGT2 using a DGT along each dimension.

L1=size(f,1);
L2=size(f,2);

N1=L1/a1;
N2=L2/a2;

c=dgt(f,g1,a1,M1);

c=reshape(c,M1*N1,L2);

c=c.';

c=dgt(c,g2,a2,M2);

c=reshape(c,M2*N2,M1*N1);

c=c.';

c=reshape(c,M1,N1,M2,N2);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function varargout = tanh(varargin)
%TANH   Hyperbolic tangent of a CHEBFUN2.
%
%   TANH(F) returns the hyperbolic tangent of a CHEBFUN2.

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

[varargout{1:nargout}] = tanh@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%Analysis of process runtimes
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 1/12/2011
classdef RunTime
    methods (Static)
        function run(fileName)
            profSummary = edu.stanford.covert.cell.sim.analysis.RunTime.calcProcessRunTimes();
            subtotals = ...
                [profSummary.evolveState]' + ...
                [profSummary.calcResourceRequirements_Current]' + ...
                [profSummary.copyToState]' + ...
                [profSummary.copyFromState]';
            total = sum(subtotals);

            fid = fopen(fileName, 'w');
            fprintf(fid, '<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">\n');
            fprintf(fid, '\n');
            fprintf(fid, '\t\n');
            fprintf(fid, '\t\t<title>%s</title>\n', 'Profile Summary');
            fprintf(fid, '\t\t\n');
            fprintf(fid, '\t\n');
            fprintf(fid, '\t\n');
            fprintf(fid, '\t\t
\n');
            fprintf(fid, '\t\t\t
\n');
            fprintf(fid, '\t\t\t\t
\n');
            fprintf(fid, '\t\t\t\t\t%s	\n', 'Process');
            fprintf(fid, '\t\t\t\t\t%s	\n', 'Evolve State');
            fprintf(fid, '\t\t\t\t\t%s	\n', 'Substrate Requirements');
            fprintf(fid, '\t\t\t\t\t%s	\n', 'Copy from State');
            fprintf(fid, '\t\t\t\t\t%s	\n', 'Copy to State');
            fprintf(fid, '\t\t\t\t\t%s	\n', 'Percent');
            fprintf(fid, '\t\t\t\t
\n');
            fprintf(fid, '\t\t\t
\n');
            fprintf(fid, '\t\t\t
\n');
            for i = 1:numel(profSummary)
                fprintf(fid, '\t\t\t\t
\n');
                fprintf(fid, '\t\t\t\t\t%s	\n', profSummary(i).name);
                fprintf(fid, '\t\t\t\t\t%.3f	\n', profSummary(i).evolveState);
                fprintf(fid, '\t\t\t\t\t%.3f	\n', profSummary(i).calcResourceRequirements_Current);
                fprintf(fid, '\t\t\t\t\t%.3f	\n', profSummary(i).copyFromState);
                fprintf(fid, '\t\t\t\t\t%.3f	\n', profSummary(i).copyToState);
                fprintf(fid, '\t\t\t\t\t%.1f	\n', subtotals(i)/total*100);
                fprintf(fid, '\t\t\t\t
\n');
            end
            fprintf(fid, '\t\t\t
\n');
            fprintf(fid, '\t\t\t
\n');
            fprintf(fid, '\t\t\t\t\t%s	\n', 'Percent');
            fprintf(fid, '\t\t\t\t\t%.1f	\n', sum([profSummary.evolveState])/total*100);
            fprintf(fid, '\t\t\t\t\t%.1f	\n', sum([profSummary.calcResourceRequirements_Current])/total*100);
            fprintf(fid, '\t\t\t\t\t%.1f	\n', sum([profSummary.copyFromState])/total*100);
            fprintf(fid, '\t\t\t\t\t%.1f	\n', sum([profSummary.copyToState])/total*100);
            fprintf(fid, '\t\t\t\t\t%.1f	\n', 100);
            fprintf(fid, '\t\t\t
\n');
            fprintf(fid, '\t\t
\n');
            fprintf(fid, '\t\n');
            fprintf(fid, '\n');
            fclose(fid);
        end

        function times = calcProcessRunTimes(profData)
            if ~exist('profData', 'var')
                profData = profile('info');
            end

            functionNames = {profData.FunctionTable.FunctionName};

            names = {};
            for i = 1:numel(profData.FunctionTable)
                processClass = regexp(profData.FunctionTable(i).FunctionName, '(.*?)>.*?\.evolveState$', 'tokens');
                if isempty(processClass); continue;
                end

                names = [names; processClass{1}(1)]; %#ok<AGROW>
            end
            names = [{'Process'}; sort(names)];

            times = struct;
            for i = 1:numel(names)
                times(i).name = names{i};

                %evolveState
                idx = find(strcmp(functionNames, [times(i).name '>' times(i).name '.evolveState']), 1);
                if ~isempty(idx)
                    times(i).evolveState = profData.FunctionTable(idx).TotalTime;
                else
                    times(i).evolveState = 0;
                end

                %copyFromState
                idx = find(strcmp(functionNames, [times(i).name '>' times(i).name '.copyFromState']), 1);
                if ~isempty(idx)
                    times(i).copyFromState = profData.FunctionTable(idx).TotalTime;
                else
                    times(i).copyFromState = 0;
                end

                %copyToState
                idx = find(strcmp(functionNames, [times(i).name '>' times(i).name '.copyToState']), 1);
                if ~isempty(idx)
                    times(i).copyToState = profData.FunctionTable(idx).TotalTime;
                else
                    times(i).copyToState = 0;
                end

                %calcResourceRequirements_Current
                idx = find(strcmp(functionNames, [times(i).name '>' times(i).name '.calcResourceRequirements_Current']), 1);
                if ~isempty(idx)
                    times(i).calcResourceRequirements_Current = profData.FunctionTable(idx).TotalTime;
                else
                    times(i).calcResourceRequirements_Current = 0;
                end
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [Eft, Covft, ljpyt, Eyt, Covyt] = gpia_jpred(gp_array, x, y, varargin)
%GPIA_JPRED  Prediction with Gaussian Process GP_IA solution.
%
%  Description
%    [EFT, COVFT] = GPIA_JPRED(GP_ARRAY, X, Y, XT, OPTIONS)
%    takes a cell array of GP structures together with matrix X of
%    training inputs and vector Y of training targets, and
%    evaluates the predictive distribution at test inputs XT with
%    parameters marginalized out with IA. Returns a posterior mean
%    EFT and covariance COVFT of latent variables.
%
%    [EFT, COVFT, JPYT, EYT, COVYT] = GPIA_JPRED(GP, X, Y, XT, 'yt', YT, ...)
%    returns also logarithm of the predictive joint density PYT of
%    the observations YT at test input locations XT with parameters
%    marginalized out with IA. This can be used for example in the
%    cross-validation. Here Y has to be vector. Returns also
%    posterior predictive mean EYT and covariance COVYT.
%
%    [EFT, COVFT, LJPYT, EYT, COVYT, FT, PFT] = ...
%      GPIA_JPRED(GP_ARRAY, X, Y, XT, OPTIONS)
%    returns also the numerical representation of the marginal
%    posterior of latent variables at each XT. FT is a vector of
%    latent values and PFT_i = p(FT_i) is the posterior density for
%    FT_i.
%
%    [EF, COVF, LJPY, EY, COVY, F, PF] = GPIA_JPRED(GP, X, Y, OPTIONS)
%    evaluates the predictive distribution at training inputs X
%    and logarithm of the predictive density PY of the training
%    observations Y.
%
%    OPTIONS is optional parameter-value pair
%      predcf - index vector telling which covariance functions are
%               used for prediction. Default is all (1:gpcfn). See
%               additional information below.
%      tstind - a vector/cell array defining, which rows of X belong
%               to which training block in *IC type sparse models.
%               Deafult is []. In case of PIC, a cell array
%               containing index vectors specifying the blocking
%               structure for test data. IN FIC and CS+FIC a vector
%               of length n that points out the test inputs that
%               are also in the training set (if none, set TSTIND=[])
%      yt     - optional observed yt in test points (see below)
%      z      - optional observed quantity in triplet (x_i,y_i,z_i)
%               Some likelihoods may use this. For example, in case
%               of Poisson likelihood we have z_i=E_i, that is,
%               expected value for ith case.
%      zt     - optional observed quantity in triplet (xt_i,yt_i,zt_i)
%               Some likelihoods may use this. For example, in case
%               of Poisson likelihood we have z_i=E_i, that is, the
%               expected value for the ith case.
%
%    NOTE! In case of FIC and PIC sparse approximation the
%    prediction for only some PREDCF covariance functions is just
%    an approximation since the covariance functions are coupled in
%    the approximation and are not strictly speaking additive
%    anymore.
%
%    For example, if you use covariance such as K = K1 + K2 your
%    predictions Eft1 = gpia_pred(gp_array, X, Y, X, 'predcf', 1) and
%    Eft2 = gpia_pred(gp_array, x, y, x, 'predcf', 2) should sum up to
%    Eft = gpia_pred(gp_array, x, y, x). That is Eft = Eft1 + Eft2. With
%    FULL model this is true but with FIC and PIC this is true only
%    approximately. That is Eft \approx Eft1 + Eft2.
%
%    With CS+FIC the predictions are exact if the PREDCF covariance
%    functions are all in the FIC part or if they are CS
%    covariances.
%
%    NOTE! When making predictions with a subset of covariance
%    functions with FIC approximation the predictive variance can
%    in some cases be ill-behaved i.e. negative or unrealistically
%    small. This may happen because of the approximative nature of
%    the prediction.
%
%  See also
%    GP_PRED, GP_SET, GP_IA
%

% Copyright (c) 2009 Ville Pietil�inen
% Copyright (c) 2009-2010 Jarno Vanhatalo
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2012 Aki Vehtari

% This software is distributed under the GNU General Public
% Licence (version 3 or later); please refer to the file
% Licence.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'GPIA_JPRED';
  ip.addRequired('gp_array', @iscell);
  ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addOptional('xt', [], @(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))))
  ip.addParamValue('yt', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('zt', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('predcf', [], @(x) isempty(x) || ...
                   isvector(x) && isreal(x) && all(isfinite(x)&x>0))
  ip.addParamValue('tstind', [], @(x) isempty(x) || iscell(x) ||...
                   (isvector(x) && isreal(x) && all(isfinite(x)&x>0)))
  if numel(varargin)==0 || isnumeric(varargin{1})
    % inputParser should handle this, but it doesn't
    ip.parse(gp_array, x, y, varargin{:});
  else
    ip.parse(gp_array, x, y, [], varargin{:});
  end
  xt=ip.Results.xt;
  yt=ip.Results.yt;
  z=ip.Results.z;
  zt=ip.Results.zt;
  predcf=ip.Results.predcf;
  tstind=ip.Results.tstind;
  if isempty(xt)
    xt=x;
    if isempty(tstind)
      if iscell(gp)
        gptype=gp{1}.type;
      else
        gptype=gp.type;
      end
      switch gptype
        case {'FULL' 'VAR' 'DTC' 'SOR'}
          tstind = [];
        case {'FIC' 'CS+FIC'}
          tstind = 1:size(x,1);
        case 'PIC'
          if iscell(gp)
            tstind = gp{1}.tr_index;
          else
            tstind = gp.tr_index;
          end
      end
    end
    if isempty(yt)
      yt=y;
    end
    if isempty(zt)
      zt=z;
    end
  end

  % pass these forward
  options=struct();
  if ~isempty(ip.Results.yt);options.yt=ip.Results.yt;end
  if ~isempty(ip.Results.z);options.z=ip.Results.z;end
  if ~isempty(ip.Results.zt);options.zt=ip.Results.zt;end
  if ~isempty(ip.Results.predcf);options.predcf=ip.Results.predcf;end
  if ~isempty(ip.Results.tstind);options.tstind=ip.Results.tstind;end

  if nargout > 2 && isempty(yt)
    pyt = NaN;
  end

  nGP = numel(gp_array);

  for i=1:nGP
    P_TH(i,:) = gp_array{i}.ia_weight;
  end

  % Make predictions with different models in gp_array
  for j = 1:nGP
    if isempty(yt)
      [Eft_grid(j,:), Covft_grid(:,:,j)]=gp_jpred(gp_array{j},x,y,xt,options);
    else
      [Eft_grid(j,:), Covft_grid(:,:,j), ljpyt_grid(j), Eyt_grid(j,:), Covyt_grid(:,:,j)]=gp_jpred(gp_array{j},x,y,xt, options);
    end
  end

  % Sample from the combined predictive distribution
  N = 100; % Sample size per gp in gp_array
  pft = zeros(N*size(Eft_grid,1), 1);
  ft = [];
  for j = 1 : size(Eft_grid,1)
    ftt = repmat(Eft_grid(j,:),N,1)+(chol(Covft_grid(:,:,j))'*randn(size(Covft_grid(:,:,j),1),N))';
    ft = [ft; ftt];
  end

  % Calculate mean and variance of the distributions
  Eft = mean(ft)';
  Covft = cov(ft);

  % Calculate jpyt with weight given in P_TH.
  if nargout > 2
    if ~isempty(yt)
      ljpyt = log(sum(exp(ljpyt_grid)'.*P_TH));
    else
      error('yt must be provided to get ljpyt');
    end
  end

  if nargout > 3
    Eyt = sum(Eyt_grid.*repmat(P_TH,1,size(Eyt_grid,2)),1);
    Covyt = zeros(size(Covyt_grid));
    for i = 1:size(Covyt_grid,3)
      Covyt(:,:,i) = Covyt_grid(:,:,i).*P_TH(i) + diag((Eyt_grid(i,:) - Eyt)).^2;
    end
    Covyt = sum(Covyt,3);
  end

  if nargout > 6
    for j=1:size(Eft_grid,1)
      for i=1:N
        index = N*(j-1)+i;
        pft(index,:) = mnorm_pdf(ft(index,:), Eft_grid(j,:), Covft_grid(:,:,j));
      end
    end
    % Normalize distributions
    pft = bsxfun(@rdivide,pft,sum(pft,1));
  end


text:
function publish_html(filename, outputDir, stylesheet)

% publish_html - publish a file to HTML format
%
%    publish_html(filename, outputDir, stylesheet);
%
%   Copyright (c) 2008 Gabriel Peyre

if nargin<1
    filename = 'content';
end
if nargin<2
    outputDir = 'html';
end
if nargin<3
    stylesheet = [outputDir '/gpeyre.xsl'];
end

opts.outputDir = outputDir;
opts.stylesheet = stylesheet;
file = publish(filename,opts);

web(file);

author: gpeyre

repo: matlab-toolboxes


text:
function p = square_uniform ( n )

%*****************************************************************************80
%
%% SQUARE_UNIFORM returns sample points from the unit square.
%
%  Discussion:
%
%    This routine returns N points sampled uniformly at random
%    from within the unit square.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 May 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points to generate.
%
%    Output, real P(N,2), the sample points.
%
  p = rand ( n, 2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [P,X] = vgg_PX_from_6pts_3img(x1,x2,x3)
% vgg_PX_from_6pts_3img  Computes camera matrices and world points
% from 6 points across 3 images.
%
%   [P,X] = vgg_PX_from_6pts_3img(x), where
%      x ... double(3,6,3) or cell{3} of double(3,6), 6 homogeneous points in 3 images
%      P ... double(3,4,3), P(:,:,k) is k-th camera matrix
%      X ... double(4,6), homogeneous world points
%   There are 0 to 3 solutions for (P,X). Solutions are pruned by requirement that
%   scalars s in all projective equations s*x==P*X are positive.
%   In case of multiple solutions, P and X have one dimension
%   more such that P(:,:,:,n) and X(:,:,n) is the n-th solution.
%
%   Also the form [P,X] = vgg_PX_from_6pts_3img(x1,x2,x3) is accepted.

% Algorithm in Hartley-Zisserman, Alg 19.1 page 493 in 1st edition,
%                                 Alg 20.1 page 511 in 2nd edition
% Coded by werner@robots.ox.ac.uk, Nov 2002.

if nargin==3
  x = cat(3,x1,x2,x3);
else
  if iscell(x1)
    x = cat(3,x1{:});
  else
    x = x1;
  end
end
if any(size(x)~=[3 6 3])
  error('Wrong size of input points.');
end

for k = 1:3
  % Find homographies H_k mapping first 4 pts in each image to standard projective basis.
  % Now, x(:,1:4,k) = H(:,:,k)*[eye(3) [1;1;1]].
  H(:,:,k) = H_from_4x(x(:,3:6,k));

  % Form transformed points xs(:,1:2,k)
  xs(:,:,k) = inv(H(:,:,k))*x(:,:,k);
end

% Compute dual fundamental matrix
Fd = Fdual_from_x(xs(:,1:2,:));

% Retrieve (non-dual) cameras P and world points X from (each solution for) dual fund. matrix.
P = [];
X = [];
for i = 1:size(Fd,3)

  % Compute canonical Xi and Pi
  Xi = X_from_Fdual(Fd(:,:,i));
  Pi = P_from_Xx_canonical(Xi,xs);

  % return from canonical to original image bases
  for k = 1:3
    Pi(:,:,k) = H(:,:,k)*Pi(:,:,k);
  end

  % Compute signs of P and X, if possible. If impossible, Pi==Xi==[].
  [Pi,Xi] = vgg_signsPX_from_x(Pi,Xi,x);
  if isempty(Pi), continue, end
  %for k = 1:3, Pi(:,:,k)*Xi ./ x(:,:,k), end  % test code

  P = cat(4,P,Pi);
  X = cat(3,X,Xi);

end

return

%%%%%%%%%%%%%%%% auxiliary functions

% Solve for dual fundamental matrix.
function F = Fdual_from_x(x)

% Linear step. After that,
% for k=1:3 and i=1:2, it is x(:,1,k)'*F_i*x(:,2,k)==0.
A = [];
for k = 1:3
  x1 = x(1,1,k);  y1 = x(2,1,k);  z1 = x(3,1,k);
  x2 = x(1,2,k);  y2 = x(2,2,k);  z2 = x(3,2,k);
  A = [A; [x1*y2-z1*y2,...
           x1*z2-z1*y2,...
           y1*x2-z1*y2,...
           y1*z2-z1*y2,...
           z1*x2-z1*y2] ];
end
[u,s,v] = svd(A,0);
v1 = v(:,end-1);
v2 = v(:,end);
FF{1} = [0 v1(1:2)'; v1(3) 0 v1(4); v1(5) -sum(v1) 0];
FF{2} = [0 v2(1:2)'; v2(3) 0 v2(4); v2(5) -sum(v2) 0];

% Non-linear step. Find linear combination of F_i having zero determinant.
% Dual fund. matrix is now  F = a*FF{1} + (1-a)*FF{2}, for each element of a.
a = vgg_singF_from_FF(FF);
for i = 1:length(a)
  F(:,:,i) = a(i)*FF{1} + (1-a(i))*FF{2};
end
return

% Given dual fundamental matrix Fd, computes cameras P and world points X such that
% for each k, P(:,:,k)*X ~ x(:,:,k).
function X = X_from_Fdual(Fd)

% Retrieve second dual camera matrix Pd2.
% It is done by considering Fdual in form
%[0 b*(d-c) -c*(d-b)
% -a*(d-c) 0 c*(d-a)
% a*(d-b) -b*(d-a) 0].

Fd_aux = reshape( Fd([4 7 1 2 1 8 1 3 6]), [3 3] );
[u,s,v] = svd(Fd_aux,0);
ABC = v(:,3); % It is a:b:c = A:B:C
[u,s,v] = svd(Fd',0);
KLM = v(:,3); % It is K:L:M = ((d-a):(d-b):(d-c)

G = [0 -ABC(3) ABC(2) 0
     ABC(3) 0 -ABC(1) 0
     -ABC(2) ABC(1) 0 0
     KLM(2) -KLM(1) 0 KLM(1)-KLM(2)
     0 KLM(3) -KLM(2) KLM(2)-KLM(3)
     -KLM(3) 0 KLM(1) KLM(3)-KLM(1) ];
[u,s,v] = svd(G,0);
abcd = v(:,end);

% The test code:
%Pd1 = [eye(3) [1;1;1]];
%Pd2 = [diag(abcd(1:3)) [1;1;1]*abcd(4)];
%vgg_F_from_P(Pd1,Pd2) ./ Fd

% Retrieve 6 world points X(:,1:6)
X = [ abcd [1;1;1;1] eye(4) ];

return

% Computes camera P from world points X and image points x, everything in canonical form.
function P = P_from_Xx_canonical(X,x)
X = X(:,1);
for k = 1:3
  A = [contreps(x(:,1,k))*[X(1) 0 0 X(4)
                           0 X(2) 0 X(4)
                           0 0 X(3) X(4)]
       contreps(x(:,2,k))*[eye(3) [1;1;1]]];
  [u,s,v] = svd(A,0);
  a = v(:,end);
  P(:,:,k) = [a(1) 0 0 a(4)
              0 a(2) 0 a(4)
              0 0 a(3) a(4)];
end
return

% A bit faster vgg_contreps.
function X = contreps(x)
X = [0 x(3) -x(2)
     -x(3) 0 x(1)
     x(2) -x(1) 0];
return

% H_from_4x  Having four point matches, the homography relating them is given by
% H = H_from_4x(x2)*inv(H_from_4x(x1)).
function H = H_from_4x(x)
H = x(:,1:3) * diag(inv(x(:,1:3))*x(:,4));
return

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

% Test code:
P = randn(3,4,3); X = randn(4,6);
for k = 1:3, x(:,:,k) = P(:,:,k)*X; end
[Pn,Xn] = vgg_PX_from_6pts_3img(x);

% check image points x predicted up to scale
for k=1:3, Pn(:,:,k)*Xn(:,:) ./ x(:,:,k), end

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
%% Systems of linear equations
% At its heart, the subject of linear algebra is concerned with how to
% solve the seemingly simple equation:
%
% $${\bf{A}}x=b$,
% where $${\bf{A}}$ is a matrix and $$x$ and $$b$ are vectors.
%%
%
% This might seem surprising, but this problem is right at the heart of
% almost all problems in scientific computing, from differential equations
% to eigenvalue problems to interpolation and curve fitting.

%%
% Of course, for a single scalar variable, x, we know what that means:
% $$ a\times x = y$
% from which we can deduce that
% $$ x = y/a$, for non-zero $$a$.
%%
% The question is how do we solve such a problem when we have more than one
% unknown variable, or, equivalently, when x is a vector? To investigate this
% question, let's consider the following electrical circuit:
%
[c,m]=imread('circuit2.gif','gif');

figure('color',[1 1 1]),
image(c), colormap(m), axis equal, axis off

%%
% We want to calculate the currents $$I_1,I_2,I_3$
%%
% Using Kirchoff's voltage laws we get the following three equations for
% the three unknown currents:
%%
% <latex>
% $\left\{\begin{array}{ccc}
% I_1 + 25(I_1-I_2)+50(I_1-I_3) & = & 10 \\
% 25(I_2-I_1) + 30I_2 + I_2 -I_3 & = & 0 \\
% 50(I_3-I_1)+I_3-I_2+55I_3 & = & 0 \end{array}\right . $
% </latex>
%%
% which, with a little re-arranging becomes:
%%
% <latex>
% $\left\{\begin{array}{ccc}
% 76I_1 - 25I_2-50I_3 & = & 10 \\
% -25I_1 +56I_2 -I_3 & = & 0 \\
% -50I_1-I_2+106I_3 & = & 0 \end{array}\right .$
% </latex>
%%
% Finally, we can gather all the coefficients into a matrix, and the
% unknowns into a vector, to write the single matrix linear equation:
%%
% <latex>
% $\left( {\begin{array}{ccc}
% 76 & -25 & -50 \\
% -25 & 56 & -1  \\
% -50& -1 & 106 \end{array}} \right)
% \times \left(\begin{array}{c}
% I_1\\
% I_2\\
% I_3 \end{array}\right) = \left( \begin{array}{c}
% 10 \\
% 0\\
% 0 \end{array}\right)$
% </latex>
%%
% In this way we have got the system of equations into the form that we
% described above, namely:
% $${\bf{A}}x=b$
%%
% Of course, for a small system (3 by 3 in this case) we can solve this by
% hand. For a larger system, however, we need a computer. The basic
% approach to solving such a system is an extension of that used by hand,
% namely we add and subtract multiples of one row to one another until we
% finally end up with an upper triangular matrix. This action
% of adding and subtrating multiples of one row to another is known as a
% basic row operation. Let us illustrate the process with MATLAB. Before we
% begin, we augment the right hand side vector to be an extra column in our
% matrix $$\bf A$, so it is now three rows by four columns.
A=[76 -25 -50 10;-25 56 -1 0; -50 -1 106 0]
rrefexample(A)
%%
% What is the advantage of reducing $${\bf A}$ to be upper triangular?
% Starting from the last row, we see that we now have an equation that
% depends only on $$I_3$, that is $$I_3 = 0.117$. The row above depends on
% $$I_3$ and $$I_2$ only, and we now know what $$I_3$ is! In this way, we
% work our way up through the rows substituting as we go, and in each row
% there is only a single unknown variable. This process is called _back
% substitution_. This entire process - that of reducing the matrix to a
% triangular form and solving via back substitution - is the way most
% computer packages solve these linear systems. In MATLAB, the shorthand
% way of solving such systems is with the backslash character ``\''
A=[76 -25 -50;-25 56 -1; -50 -1 106]
b =[10;0;0]
I= A\b
%%
% In the example above, we have performed a special case of what is a more
% general idea - that of transforming a matrix into a special form, so that
% the equations can be solved more easily. In general, the idea of
% transforming a matrix through a series of elementary row operations is a
% powerful one. May special types of _matrix factorisations_ exist, and two
% of the most important examples are:
%%
% *QR decomposition *
% In this case, a general matrix is written as the product of two matrices,
% one orthogonal and one upper triangular:
% $${\bf A = QR}$
% As the inverse of an orthogonal matrix is equal to its transpose, we find
% the solution to the general system $${\bf A}x=b$ can be calculated as:
%%
% <latex>
% $\begin{array}{l}
% {\bf A}x = b \\
% {\bf QR}x = b \\
% {\bf Q}({\bf R}x)=b \\
% {\bf R}x = {\bf Q^T} \times b \\
% x = {\bf R}\setminus{\bf Q^T} \times b \end{array}$
% </latex>
%%
% *LU Decomposition*
% In this case the matrix is decomposed into two, a lower and an upper
% triangular matrix. As both the matrices are triangular, the back (and
% forward) substitution process is very quick:
%%
% <latex>
% $\begin{array}{l}
% {\bf A}x = b \\
% {\bf LU}x = b \\
% {\bf U}({\bf L}x)=b \\
% {\bf U}x = {\bf L}\setminus b \\
% x = {\bf U}\setminus({\bf L}\setminus b) \end{array}$
% </latex>

%   Copyright 2008-2009 The MathWorks, Inc.
%   $Revision: 35 $  $Date: 2009-05-29 15:27:34 +0100 (Fri, 29 May 2009) $

author: Sable

repo: mcbench-benchmarks