text:
classdef MergeRobertson < handle
    %MERGEROBERTSON  Merge exposure sequence to a single image
    %
    % The resulting HDR image is calculated as weighted average of the
    % exposures considering exposure values and camera response.
    %
    % For more information see [RB99].
    %
    % ## References
    % [RB99]:
    % > Mark A Robertson, Sean Borman, and Robert L Stevenson.
    % > "Dynamic range improvement through multiple exposures".
    % > In Image Processing, 1999. ICIP 99. Proceedings. 1999 International
    % > Conference on, volume 3, pages 159-163. IEEE, 1999.
    %
    % See also: cv.MergeDebevec, cv.MergeMertens, makehdr
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    %% MergeRobertson
    methods
        function this = MergeRobertson()
            %MERGEROBERTSON  Creates MergeRobertson object
            %
            %     obj = cv.MergeRobertson()
            %
            % See also: cv.MergeRobertson.process
            %
            this.id = MergeRobertson_(0, 'new');
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.MergeRobertson
            %
            if isempty(this.id), return; end
            MergeRobertson_(this.id, 'delete');
        end
    end

    %% MergeExposures
    methods
        function dst = process(this, src, etimes, varargin)
            %PROCESS  Merges images
            %
            %     dst = obj.process(src, etimes)
            %     dst = obj.process(src, etimes, response)
            %
            % ## Input
            % * __src__ vector of input images, all of the same size and
            %   `uint8` type.
            % * __etimes__ vector of exposure time values for each image.
            % * __response__ 256x1xCN `single` matrix with inverse camera
            %   response function (CRF) for each pixel value, it should have
            %   the same number of channels as images `src{i}`.
            %
            % ## Output
            % * __dst__ result image, same size as `src{i}` and `single` type.
            %
            % The function has a short version, that doesn't take the extra
            % `response` argument.
            %
            % See also: cv.MergeRobertson.MergeRobertson
            %
            dst = MergeRobertson_(this.id, 'process', src, etimes, varargin{:});
        end
    end

    %% Algorithm
    methods (Hidden)
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.MergeRobertson.empty, cv.MergeRobertson.load
            %
            MergeRobertson_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Returns true if the algorithm is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.MergeRobertson.clear, cv.MergeRobertson.load
            %
            b = MergeRobertson_(this.id, 'empty');
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.MergeRobertson.save, cv.MergeRobertson.load
            %
            name = MergeRobertson_(this.id, 'getDefaultName');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.MergeRobertson.load
            %
            MergeRobertson_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.MergeRobertson.save
            %
            MergeRobertson_(this.id, 'load', fname_or_str, varargin{:});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function res = parallelEdge(edge, dist)
%PARALLELEDGE Edge parallel to another edge.
%
%   EDGE2 = parallelEdge(EDGE, DIST)
%   Computes the edge parallel to the input edge EDGE and located at the
%   given signed distance DIST.
%
%   Example
%     obox = [30 40 80 40 30];
%     figure; hold on; axis equal;
%     drawOrientedBox(obox, 'LineWidth', 2);
%     edge1 = centeredEdgeToEdge(obox([1 2 3 5]));
%     edge2 = centeredEdgeToEdge(obox([1 2 4 5])+[0 0 0 90]);
%     drawEdge(edge1, 'LineWidth', 2, 'color', 'g');
%     drawEdge(edge2, 'LineWidth', 2, 'color', 'g');
%     drawEdge(parallelEdge(edge1, -30), 'LineWidth', 2, 'color', 'k');
%     drawEdge(parallelEdge(edge2, -50), 'LineWidth', 2, 'color', 'k');
%
%   See also
%     edges2d, parallelLine, drawEdge, centeredEdgeToEdge, edgeToLine

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2012-07-31, using Matlab 7.9.0.529 (R2009b)
% Copyright 2012-2022 INRA - Cepia Software Platform

% compute the line parallel to the supporting line of edge
line = parallelLine(edgeToLine(edge), dist);

% result edge is given by line positions 0 and 1.
res = [line(:, 1:2) line(:, 1:2)+line(:, 3:4)];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function geometry_test0325 ( )

%*****************************************************************************80
%
%% TEST0325 tests ICOS_SIZE, ICOS_SHAPE, SHAPE_PRINT_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0325\n' );
  fprintf ( 1, '  For the icosahedron,\n' );
  fprintf ( 1, '  ICOS_SIZE returns dimension information;\n' );
  fprintf ( 1, '  ICOS_SHAPE returns face and order information.\n' );
  fprintf ( 1, '  SHAPE_PRINT_3D prints this information.\n' );

  [ point_num, edge_num, face_num, face_order_max ] = icos_size ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of points =   %d\n', point_num );
  fprintf ( 1, '  Number of edges =    %d\n', edge_num );
  fprintf ( 1, '  Number of faces =    %d\n', face_num );
  fprintf ( 1, '  Maximum face order = %d\n', face_order_max );

  [ point_coord, edge_point, face_order, face_point ] = icos_shape ( ...
    point_num, edge_num, face_num, face_order_max );

  shape_print_3d ( point_num, face_num, face_order_max, ...
    point_coord, face_order, face_point );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [poses] = InterpolatePoses(ins_file, pose_timestamps, origin_timestamp, use_rtk)

% InterpolatePoses - interpolate INS poses to find poses at given timestamps
%
% [poses] = InterpolatePoses(ins_file, pose_timestamps, origin_timestamp)
%
% INPUTS:
%   ins_file: csv file containing INS data
%   pose_timestamps: array of UNIX timestamps at which interpolated poses are
%     required
%   origin_timestamp: timestamp for origin frame, relative to which poses are
%     reported
%
% OUTPUTS:
%   poses: cell array of 4x4 matrices, representing SE3 poses at the times
%     specified in pose_timestamps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (c) 2016 University of Oxford
% Authors:
%  Geoff Pascoe (gmp@robots.ox.ac.uk)
%  Will Maddern (wm@robots.ox.ac.uk)
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License.
% To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to
% Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  ins_file_id = fopen(ins_file);
  if ~use_rtk
      header_count = 15;
  else
      header_count = 23;
  end
  headers = textscan(ins_file_id, '%s', header_count, 'Delimiter',',');
  if ~use_rtk
      format_str = '%u64 %s %f %f %f %f %f %f %s %f %f %f %f %f %f';
  else
      format_str = '%u64 %f %f %f %f %f %f %s %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f';
  end
  ins_data = textscan(ins_file_id, ...
      format_str,'Delimiter',',');
  fclose(ins_file_id);

  ins_timestamps = ins_data{1};

  lower_index = max(find(ins_timestamps(:,1)<=min(pose_timestamps), 1, ...
      'last')-1,1);
  upper_index = find(ins_timestamps(:,1)>max(pose_timestamps), 1, 'first')+1;

  lower_index = max(lower_index, 1);
  upper_index = min(upper_index, numel(ins_timestamps));

  ins_poses = cell(1, upper_index - lower_index + 1);
  ins_quaternions = cell(1, upper_index - lower_index + 1);

  if ~use_rtk
    northing_col = 6;
    easting_col = 7;
    down_col = 8;
    roll_col = 13;
    pitch_col = 14;
    yaw_col = 15;
  else
    northing_col = 5;
    easting_col = 6;
    down_col = 7;
    roll_col = 12;
    pitch_col = 13;
    yaw_col = 14;
  end
  northings = ins_data{northing_col};
  eastings = ins_data{easting_col};
  downs = ins_data{down_col};
  rolls = ins_data{roll_col};
  pitches = ins_data{pitch_col};
  yaws = ins_data{yaw_col};

  pose_timestamps = [origin_timestamp pose_timestamps];

  for i=1:upper_index
    ins_poses{i} = SE3MatrixFromComponents(...
      northings(i), eastings(i), downs(i), rolls(i), pitches(i), yaws(i));
    ins_quaternions{i} = SO3ToQuaternion(ins_poses{i}(1:3,1:3))';
  end

  [lower_index_rows, lower_index_cols] = ...
    find(and(bsxfun(@le,ins_timestamps,pose_timestamps)',...
    circshift(bsxfun(@gt,ins_timestamps,pose_timestamps)',[0 -1])));

  lower_indices = zeros(size(lower_index_rows));
  lower_indices(lower_index_rows) = lower_index_cols;
  lower_indices = max(lower_indices, 1);

  ins_timestamps = cast(ins_timestamps, 'double');
  pose_timestamps = cast(pose_timestamps, 'double');
  fractions = (pose_timestamps - ins_timestamps(lower_indices)')./...
    (ins_timestamps(lower_indices+1)'-ins_timestamps(lower_indices)');

  quaternions_lower = [ins_quaternions{lower_indices}];
  quaternions_upper = [ins_quaternions{lower_indices+1}];

  % interpolate quaternions
  d_array = sum(quaternions_lower.*quaternions_upper,1);

  linear_interp_indices = find(~(abs(d_array)<1.0));
  sin_interp_indices = find(abs(d_array)<1.0);

  scale0_array = zeros(size(d_array));
  scale1_array = zeros(size(d_array));

  scale0_array(linear_interp_indices) = 1-fractions(linear_interp_indices);
  scale1_array(linear_interp_indices) = fractions(linear_interp_indices);

  theta_array = acos(abs(d_array(sin_interp_indices)));

  scale0_array(sin_interp_indices) = sin((1-fractions(sin_interp_indices)).*...
      theta_array)./sin(theta_array);
  scale1_array(sin_interp_indices) = sin(fractions(sin_interp_indices).*...
      theta_array)./sin(theta_array);

  negative_d_indices = find(d_array < 0);
  scale1_array(negative_d_indices) = -scale1_array(negative_d_indices);

  quaternions_interp = repmat(scale0_array,4,1).*quaternions_lower + ...
      repmat(scale1_array,4,1).*quaternions_upper;

  % interpolate positions
  ins_poses_array = [ins_poses{:}];
  positions_lower = ins_poses_array(1:3,4*lower_indices);
  positions_upper = ins_poses_array(1:3,4*(lower_indices+1));
  positions_interp = repmat((1-fractions),3,1).*positions_lower + ...
      repmat(fractions,3,1).*positions_upper;

  poses_array = zeros(4, 4*numel(pose_timestamps));
  poses_array(1,1:4:end) = 1-2*quaternions_interp(3,:).^2 ...
      -2*quaternions_interp(4,:).^2;
  poses_array(1,2:4:end) = 2*quaternions_interp(2,:).*...
    quaternions_interp(3,:) - 2*quaternions_interp(4,:).*...
    quaternions_interp(1,:);
  poses_array(1,3:4:end) = 2*quaternions_interp(2,:).*...
      quaternions_interp(4,:) + 2*quaternions_interp(3,:).*...
      quaternions_interp(1,:);

  poses_array(2,1:4:end) = 2*quaternions_interp(2,:).*...
      quaternions_interp(3,:) + 2*quaternions_interp(4,:).*...
      quaternions_interp(1,:);
  poses_array(2,2:4:end) = 1-2*quaternions_interp(2,:).^2 ...
      -2*quaternions_interp(4,:).^2;
  poses_array(2,3:4:end) = 2*quaternions_interp(3,:).*...
      quaternions_interp(4,:) - 2*quaternions_interp(2,:).*...
      quaternions_interp(1,:);

  poses_array(3,1:4:end) = 2*quaternions_interp(2,:).*...
      quaternions_interp(4,:) - 2*quaternions_interp(3,:).*...
      quaternions_interp(1,:);
  poses_array(3,2:4:end) = 2*quaternions_interp(3,:).*...
      quaternions_interp(4,:) + 2*quaternions_interp(2,:).*...
      quaternions_interp(1,:);
  poses_array(3,3:4:end) = 1-2*quaternions_interp(2,:).^2 ...
      -2*quaternions_interp(3,:).^2;

  poses_array(1:3,4:4:end) = positions_interp;
  poses_array(4,4:4:end) = 1;
  poses_array = poses_array(1:4,1:4) \ poses_array;

  poses = cell(numel(pose_timestamps)-1,1);
  for i=2:numel(pose_timestamps)
    poses{i-1} = poses_array(:,(i-1)*4+1:i*4);
  end

end

author: ori-mrg

repo: robotcar-dataset-sdk


text:
function prediction = VGG8(varargin)
%VGG8 Returns a VGG-S/M/F (VGG-8) model for ImageNet
%   M = models.VGG8() returns the model VGG-F proposed in:
%
%     Chatfield et al., "Return of the Devil in the Details: Delving Deep
%     into Convolutional Networks", BMVC 2014.
%
%   models.VGG8(..., 'option', value, ...) accepts the following options:
%
%   `variant`:: 'f'
%     Model variant: s, m or f (slow/medium/fast).
%
%   `pretrained`:: false
%     If true, returns a model pre-trained on ImageNet (using the
%     MatConvNet example code).
%
%   `input`:: default input
%     Specifies an input (images) layer for the network. If unspecified, a
%     new one is created.
%
%   `numClasses`:: 1000
%     Number of output classes.
%
%   `batchNorm`:: true
%     Whether to use batch normalization.
%
%   `normalization`:: [5 1 0.0001/5 0.75]
%     Parameters for vl_nnnormalize layer (only used without batch-norm).
%
%   Any other options will be passed to models.ConvBlock(), and can be used
%   to change the activation function, weight initialization, etc.
%
%   Suggested SGD training options are also returned in the struct M.meta.

% Copyright (C) 2018 Joao F. Henriques, Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

  % parse options. unknown arguments will be passed to ConvBlock (e.g.
  % activation).
  opts.pretrained = false ;  % whether to fetch a pre-trained model
  opts.input = Input('name', 'images', 'gpu', true) ;  % default input layer
  opts.numClasses = 1000 ;  % number of predicted classes
  opts.variant = 'f' ;  % choose between variants s/m/f (slow/medium/fast)
  opts.batchNorm = true ;  % whether to use batch normalization
  opts.normalization = [5 1 0.0001/5 0.75] ;  % for LRN layer (vl_nnnormalize)
  [opts, convBlockArgs] = vl_argparse(opts, varargin, 'nonrecursive') ;

  % default training options for this network (returned as output.meta)
  switch lower(opts.variant)
  case 's'
    meta.batchSize = 128 ;
  case 'm'
    meta.batchSize = 196 ;
  case 'f'
    meta.batchSize = 256 ;
  otherwise
    error('Unknown variant.') ;
  end
  meta.imageSize = [224, 224, 3] ;
  meta.augmentation.crop = 224 / 256;
  meta.augmentation.location = true ;
  meta.augmentation.flip = true ;
  meta.augmentation.brightness = 0.1 ;
  meta.augmentation.aspect = [2/3, 3/2] ;
  meta.weightDecay = 0.0005 ;

  % the default learning rate schedule
  if ~opts.pretrained
    if ~opts.batchNorm
      meta.learningRate = logspace(-2, -4, 60) ;
    else
      meta.learningRate = logspace(-1, -4, 20) ;
    end
    meta.numEpochs = numel(meta.learningRate) ;
  else  % fine-tuning has lower LR
    meta.learningRate = 1e-5 ;
    meta.numEpochs = 20 ;
  end

  % return a pre-trained model
  if opts.pretrained
    if opts.batchNorm
      warning('The pre-trained model does not include batch-norm (set batchNorm to false).') ;
    end
    if opts.numClasses ~= 1000
      warning('Model options are ignored when loading a pre-trained model.') ;
    end
    prediction = models.pretrained(['imagenet-vgg-' opts.variant]) ;

    % return prediction layer (not softmax)
    assert(isequal(prediction{1}.func, @vl_nnsoftmax)) ;
    prediction = prediction{1}.inputs{1} ;

    % replace input layer with the given one
    input = prediction.find('Input', 1) ;
    prediction.replace(input, opts.input) ;

    prediction.meta = meta ;
    return
  end

  % get conv block generator with the given options. default activation is
  % ReLU, with pre-activation batch normalization (can be overriden).
  conv = models.ConvBlock('batchNorm', opts.batchNorm, ...
    'preActivationBatchNorm', true, convBlockArgs{:}) ;

  % build network
  images = opts.input ;

  % implement the 3 variants: S, M and F (from paper)
  switch lower(opts.variant)
  case 's'
    % first conv block
    x = conv(images, 'size', [7, 7, 3, 96], 'stride', 2) ;
    if ~opts.batchNorm
      x = vl_nnnormalize(x, opts.normalization) ;
    end
    x = vl_nnpool(x, 3, 'stride', 3, 'pad', [0 2 0 2]) ;

    % second conv block
    x = conv(x, 'size', [5, 5, 96, 256]) ;
    if ~opts.batchNorm
      x = vl_nnnormalize(x, opts.normalization) ;
    end
    x = vl_nnpool(x, 2, 'stride', 2, 'pad', [0 1 0 1]) ;

    % conv blocks 3-5
    x = conv(x, 'size', [3, 3, 256, 512], 'pad', 1) ;
    x = conv(x, 'size', [3, 3, 512, 512], 'pad', 1) ;
    x = conv(x, 'size', [3, 3, 512, 512], 'pad', 1) ;
    x = vl_nnpool(x, 3, 'stride', 3, 'pad', [0 1 0 1]) ;

    % first fully-connected block
    x = conv(x, 'size', [6, 6, 512, 4096]) ;

  case 'm'
    % first conv block
    x = conv(images, 'size', [7, 7, 3, 96], 'stride', 2) ;
    if ~opts.batchNorm
      x = vl_nnnormalize(x, opts.normalization) ;
    end
    x = vl_nnpool(x, 3, 'stride', 2) ;

    % second conv block
    x = conv(x, 'size', [5, 5, 96, 256], 'stride', 2, 'pad', 1) ;
    if ~opts.batchNorm
      x = vl_nnnormalize(x, opts.normalization) ;
    end
    x = vl_nnpool(x, 3, 'stride', 2, 'pad', [0 1 0 1]) ;

    % conv blocks 3-5
    x = conv(x, 'size', [3, 3, 256, 512], 'pad', 1) ;
    x = conv(x, 'size', [3, 3, 512, 512], 'pad', 1) ;
    x = conv(x, 'size', [3, 3, 512, 512], 'pad', 1) ;
    x = vl_nnpool(x, 3, 'stride', 2) ;

    % first fully-connected block
    x = conv(x, 'size', [6, 6, 512, 4096]) ;

  case 'f'
    % first conv block
    x = conv(images, 'size', [11, 11, 3, 64], 'stride', 4) ;
    if ~opts.batchNorm
      x = vl_nnnormalize(x, opts.normalization) ;
    end
    x = vl_nnpool(x, 3, 'stride', 2, 'pad', [0 1 0 1]) ;

    % second conv block
    x = conv(x, 'size', [5, 5, 64, 256], 'pad', 2) ;
    if ~opts.batchNorm
      x = vl_nnnormalize(x, opts.normalization) ;
    end
    x = vl_nnpool(x, 3, 'stride', 2) ;

    % conv blocks 3-5
    x = conv(x, 'size', [3, 3, 256, 256], 'pad', 1) ;
    x = conv(x, 'size', [3, 3, 256, 256], 'pad', 1) ;
    x = conv(x, 'size', [3, 3, 256, 256], 'pad', 1) ;
    x = vl_nnpool(x, 3, 'stride', 2) ;

    % first fully-connected block
    x = conv(x, 'size', [6, 6, 256, 4096]) ;
  end

  % finish first fully-connected block
  if ~opts.batchNorm
    x = vl_nndropout(x) ;
  end

  % second fully-connected block
  x = conv(x, 'size', [1, 1, 4096, 4096]) ;
  if ~opts.batchNorm
    x = vl_nndropout(x) ;
  end

  % prediction layer
  prediction = conv(x, 'size', [1, 1, 4096, opts.numClasses], ...
    'batchNorm', false, 'activation', 'none') ;

  prediction.meta = meta ;

end

repo: MFT


text:
function plotgm(X, model)
% Plot 2d Gaussian mixture model.
% Written by Mo Chen (sth4nth@gmail.com).
level = 64;
n = 256;

spread(X);
x_range = xlim;
y_range = ylim;

x = linspace(x_range(1),x_range(2), n);
y = linspace(y_range(2),y_range(1), n);

[a,b] = meshgrid(x,y);
z = exp(loggmpdf([a(:)';b(:)'],model));

z = z-min(z);
z = floor(z/max(z)*(level-1));

figure;
image(reshape(z,n,n));
colormap(jet(level));
set(gca, 'XTick', [1 256]);
set(gca, 'XTickLabel', [min(x) max(x)]);
set(gca, 'YTick', [1 256]);
set(gca, 'YTickLabel', [min(y) max(y)]);
axis off

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function st_to_mm_test ( )

%*****************************************************************************80
%
%% ST_TO_MM_TEST tests the ST_TO_MM library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2008
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ST_TO_MM_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the ST_TO_MM library.\n' );

  st_to_mm_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ST_TO_MM_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function mappedX = diffusion_maps(X, no_dims, t, sigma)
%DIFFUSION_MAPS Runs the diffusion map algorithm
%
%   mappedX = diffusion_maps(X, no_dims, t, sigma)
%
% The functions runs the diffusion map algorithm on dataset X to reduce it
% to dimensionality no_dims. The variable sigma is the variance of the Gaussian
% used in the affinity computation (default = 1). The variable alpha
% determines the operator that is applied on the graph (default = 1).
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    % Give memory warning
    if size(X, 1) > 3000
        warning(['Due to the large number of instances (' num2str(size(X, 1)) '), diffusion maps may run out of memory.']);
    end

    % Normalize data
    X = double(X);
    X = X - min(X(:));
    X = X / max(X(:));

    % Compute Gaussian kernel matrix
    disp(['Compute Markov forward transition probability matrix with ' num2str(t) ' timesteps...']);
    sumX = sum(X .^ 2, 2);
    K = exp(-bsxfun(@plus, sumX, bsxfun(@plus, sumX', -2 * (X * X'))) ./ (2 .* sigma ^ 2));

    % Compute Markov probability matrix with t timesteps
    p = sum(K, 1)';
    K = K ./ ((p * p') .^ t);
    p = sqrt(sum(K, 1))';
    K = K ./ (p * p');

    % Perform economy-size SVD
    disp('Perform eigendecomposition...');
    [U, S, V] = svd(K, 0);
    U = bsxfun(@rdivide, U, U(:,1));
    mappedX = U(:,2:no_dims + 1);

author: epfl-lasa


text:
% edgeNL: Edge detection function based on nonlinear derivatives
%   (elementary demo):
%   This function differientiates the image to estimate the gradient.
%   The principle is based on polarized derivatives to automatically select
%   the best edge localization.
%   The main benefits are
%     - univocal edge localization for synthetic and real images
%     - noise reduction with no regularization: the noise level is weaker than
%       the noise level in the original image
%     - better direction estimation of the gradient
%     - can product a confident edge reference map for synthetic images
%     - extremely efficient on salt noise OR pepper noise (this last case needs
%       a change in the nonlinear derivatives)
%     - still noise reduction with regularized schemes (Canny, Demigny, ...,
%       can also be adapted to the asymetrical filters (Prewitt, Sobel, ...)
%   Drawback
%     - no detection of vertical and horizontal "white" thin (1 pixel) lines
%   Rk: this demo only performs edge detection and does not include edge
%       extraction (local maxima) and other steps to obtain a binary
%       edge map.
%   Written by O. Laligant - 2009, University of Burgundy
%   Ref: A Nonlinear Derivative Scheme Applied to Edge Detection,
%       Olivier Laligant, Frederic Truchetet, IEEE Transactions on Pattern
%       Analysis and Machine Intelligence - PAMI , vol. 32, no. 2,
%       pp. 242-257, 2010
function edgeNL()

%  ------------- Edge detection on a simple synthetic image  -------------
Im = example();
figure(1), imagesc(Im);
title ('Original image');

% edge detection (nonlinear gradient estimate)
% better localization and direction estimation than the linear scheme
gI = algoNL(Im);

% For a synthetic object, the edges are localized inside the object shape
figure(2), imagesc(gI);
title('Edge detection');

% ------------- Additive gaussian white noise ------------------------

Imn = imnoise(Im, 'gaussian', 0, 0.01);
figure(3), imagesc(Imn);
title ('Noisy image');

gIn = algoNL(Imn);

figure(4), imagesc(gIn);
title('Edge detection on the noisy image');

end
% ----------------------------------------------------------------------------

%
%  --------------------------- functions  -------------------------------------
%

% ------------- computation of the nonlinear derivatives  -------------
% polarized derivatives
% lead to an univocal localization of edges
function [gm, gh, gv] = algoNL(I);
% for classical regularization schemes, I can be replaced by a regularized
%   version
% for asymetrical schemes (Prewitt, Sobel, etc), I must be replaced by two
%   regularized versions (the regularization is different on columns and rows)
dph = thresh0(conv2(I, [0 1 -1], 'same'));
dnh = -thresh0(-conv2(I, [1 -1 0], 'same'));
gh = dph+dnh;
dpv = thresh0(conv2(I, [0; 1; -1], 'same'));
dnv = -thresh0(-conv2(I, [1;  -1; 0], 'same'));
gv = dpv + dnv;
gm = sqrt(gh.*gh + gv.*gv);

end

% ------------- threshold -------------
function st = thresh0(s)
	st = s.*(sign(s)+1)/2;
end%function

% ------------- test image -------------
% an example of image to illustrate the localization of the method
function Im = example()
Im =[
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0   255   255   255     0     0     0     0     0     0     0   255     0     0
     0     0   255   255   255     0     0     0     0     0     0   255   255     0     0
     0     0   255   255   255     0     0     0     0     0   255   255   255     0     0
     0     0     0     0     0     0     0     0     0   255   255   255   255     0     0
     0     0     0     0     0     0     0     0   255   255   255   255   255     0     0
     0     0     0     0     0     0     0   255   255   255   255     0     0     0     0
     0     0     0     0     0     0   255   255   255   255     0     0     0     0     0
     0     0     0     0     0   255   255   255   255   255     0     0     0     0     0
     0     0     0     0   255   255   255     0   255   255     0     0     0     0     0
     0     0     0   255   255   255     0     0     0   255   255     0     0     0     0
     0     0   255   255   255   255   255     0   255   255   255   255   255     0     0
     0     0     0   255   255   255   255   255   255   255   255   255     0     0     0
     0     0     0     0   255   255   255   255   255   255   255   255   255     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
];
end

author: Sable

repo: mcbench-benchmarks


text:
function M = mult(A, f, lambda)
%MULT   Multiplication operator for the ultraspherical spectral method.
%   M = MULT(A, F, lambda) returns the multiplication operator that represents
%   u(x) -> F(x)u(x), in the C^{(lambda)} ultraspherical polynomial basis.
%
%   If lambda = 0, then the operator is Toeplitz-plus-Hankel-plus-rank-1 and
%   represents multiplication in Chebyshev T coefficients.
%
%   If lambda = 1, then the operator is Toeplitz-plus-Hankel and represents
%   multiplication in Chebyshev U or C^{(1)} coefficients.
%
%   If lambda > 1, then the operator does not have any Toeplitz/Hankel structure
%   and is constructed using a three-term recurrence.

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

% Obtaining some useful information:
n = A.dimension;
d = A.domain;
f = restrict(f, d);
numIntervals = length(d) - 1;

% Find the diagonal blocks;
blocks = cell(numIntervals);
for k = 1:numIntervals
    blocks{k} = ultraS.multmat(n(k), f.funs{k}, lambda);
end

% Assemble:
M = blkdiag(blocks{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function mcrop = matcrop_withpad(m,varargin)

subs = varargin;

ndim = length(subs);
sz = size(m);
if sz(1) == 1 && ndim == 1,
  subs = [1,subs];
  ndim = 2;
elseif length(sz) ~= ndim,
  error('size(m) does not match number of input indices');
end

mcrop = m;
for dim = 1:ndim,
  if subs{dim}(1) > subs{dim}(2),
    mcrop = eval(['mcrop(',repmat(':,',[1,dim-1]),'[]',repmat(',:',[1,ndim-dim]),');']);
    continue;
  end
  npadprev = max(0,1 - subs{dim}(1));
  npadafter = max(0,subs{dim}(2)-sz(dim));
  subs{dim}(1) = max(1,subs{dim}(1));
  subs{dim}(2) = min(sz(dim),subs{dim}(2));
  mcrop = eval(['mcrop(',repmat(':,',[1,dim-1]),sprintf('%d:%d',subs{dim}(1),subs{dim}(2)),...
    repmat(',:',[1,ndim-dim]),');']);
  if npadprev > 0,
    tmppad = size(mcrop);
    tmppad(dim) = npadprev;
    mcrop = cat(dim,zeros(tmppad),mcrop);
  end
  if npadafter > 0,
    tmppad = size(mcrop);
    tmppad(dim) = npadafter;
    mcrop = cat(dim,mcrop,zeros(tmppad));
  end
end

author: kristinbranson

repo: JAABA


text:
function [S, V, D, Sigma2] = MySVD(A)
[m, n] = size(A);
if 2*m < n
    AAT = A*A';
    [S, Sigma2, D] = svd(AAT);
    Sigma2 = diag(Sigma2);
    V = sqrt(Sigma2);
    tol = max(size(A)) * eps(max(V));
    R = sum(V > tol);
    V = V(1:R);
    S = S(:,1:R);
    D = A'*S*diag(1./V);
    V = diag(V);
    return;
end
if m > 2*n
    [S, V, D, Sigma2] = MySVD(A');
    mid = D;
    D = S;
    S = mid;
    return;
end
[S,V,D] = svd(A);
Sigma2 = diag(V).^2;

author: andrewssobral


text:
function [nam] = get_name(dat)

%   get_name(DATA)      returns name describing the object
  nam=dat.algorithm.name;

%%%% SGE support
if isdeferred(dat), nam = get_name(dat.algorithm.deferred, nam, 1); end
%%%%

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
    function [yM] = model(xm, a)
    %
    % Model (linear and nonlinear) definitions
    % for the function numerFminS() as a 'fun'
    %
    % Linear (dataLRM.txt):
    %yM = zeros(2,1);
    %yM=a(1)*xm+a(2);
    % Second degree polynomial (dataNRM.txt):
    yM = zeros(3,1);
    yM=a(1)*xm.^2+a(2)*xm+a(3);
    % Power exponet (dataERM.txt):
    %yM = zeros(2,1);
    %yM=a(1)*xm.^a(2);


text:
% Test file for chebtech/isempty.m

function pass = test_isempty(varargin)

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    f = testclass.make();
    pass(n, 1) = isempty(f);

    f = testclass.make(@sin);
    pass(n, 2) = ~isempty(f);

    f = testclass.make(@(x) [sin(x), cos(x)]);
    pass(n, 3) = ~isempty(f);

    f = [ testclass.make(@sin), testclass.make(@sin) ];
    pass(n, 4) = ~isempty(f);

    f = [ testclass.make() testclass.make() ];
    pass(n, 5) = isempty(f);
end

end


text:
function [ y, m, d, f ] = jed_to_ymdf_ethiopian ( jed )

%*****************************************************************************80
%
%% JED_TO_YMDF_ETHIOPIAN converts a JED to an Ethiopian YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Algorithm F,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, pages 324-325.
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, integer Y, M, D, real F, the YMDF date.
%

%
%  Determine the computational date (Y'/M'/D').
%
  j = floor ( jed + 0.5 );
  f = ( jed + 0.5 ) - j;

  j_prime = j + 124;

  y_prime = floor ( ( 4 * j_prime + 3 ) / 1461 );
  t_prime = floor ( mod ( 4 * j_prime + 3, 1461 ) / 4 );
  m_prime = floor ( t_prime / 30 );
  d_prime = mod ( t_prime, 30 );
%
%  Convert the computational date to a calendar date.
%
  d = d_prime + 1;
  m = mod ( m_prime, 13 ) + 1;
  y = y_prime - 4720 + floor ( ( 13 - m ) / 13 );

  return
end


text:
 function back = Gsub_back(Gsub, yi)
%
%	backprojector for Gsub object
%	in:
%		Gsub		Gsubset object
%		yi	[?,na]		data
%	out:
%		back	[np,2]	backprojection G'y
%
%	Copyright 2002-1-30	Jeff Fessler	The University of Michigan

if nargin < 2, ir_usage, end

na = Gsub.na;
np = Gsub.np;
nsubset = Gsub.nsubset;

yi = reshape(yi, [numel(yi)/na na]);	% [?,na]

%
%	loop over subsets
%
back = zeros(np, 1);
for iset=1:nsubset
	ia = Gsub.starts(iset):nsubset:na;
%	ia = Gsub.ia{iset};
	back = back + Gsub.cell{iset} * col(yi(:,ia));
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function varargout = visualization(mode,image_dir,varargin)

switch mode

  % init figure
  case 'init'

    % create figure using size of first image in repository
    fig = figure(1);
    img = imread(sprintf('%s/%010d.png',image_dir,0));
    set(fig,'position',[100,100,0.8*size(img,2),0.8*2*size(img,1)]);
    h(1).axes = axes('position',[0,0.5,1,0.5]);
    h(2).axes = axes('position',[0,0,1,0.5]);
    varargout{1} = h;

  % update figure
  case 'update'

    % unpack input arguments
    h        = varargin{1};
    img_idx  = varargin{2};
    nimages  = varargin{3};

    % read image
    img = imread(sprintf('%s/%010d.png',image_dir,img_idx));

    % clear axes, draw image
    cla(h(1).axes); cla(h(2).axes);
    imshow(img,'parent',h(1).axes); axis(h(1).axes,'image','off'); hold(h(1).axes, 'on');
    imshow(img,'parent',h(2).axes); axis(h(2).axes,'image','off'); hold(h(2).axes, 'on');

    % title
    text(size(img,2)/2,3,sprintf('2D Bounding Boxes'),'parent',h(1).axes,'color','g','HorizontalAlignment','center','VerticalAlignment','top','FontSize',14,'FontWeight','bold','BackgroundColor','black');
    text(size(img,2)/2,3,sprintf('3D Bounding Boxes'),'parent',h(2).axes,'color','g','HorizontalAlignment','center','VerticalAlignment','top','FontSize',14,'FontWeight','bold','BackgroundColor','black');

    % legend
    text(0,00,'Not occluded','parent',h(1).axes,'color','g','HorizontalAlignment','left','VerticalAlignment','top','FontSize',14,'FontWeight','bold','BackgroundColor','black');
    text(0,30,'Partly occluded','parent',h(1).axes,'color','y','HorizontalAlignment','left','VerticalAlignment','top','FontSize',14,'FontWeight','bold','BackgroundColor','black');
    text(0,60,'Fully occluded','parent',h(1).axes,'color','r','HorizontalAlignment','left','VerticalAlignment','top','FontSize',14,'FontWeight','bold','BackgroundColor','black');
    text(0,90,'Unknown','parent',h(1).axes,'color','w','HorizontalAlignment','left','VerticalAlignment','top','FontSize',14,'FontWeight','bold','BackgroundColor','black');

    % frame number
    text(size(img,2),0,sprintf('frame %d/%d',img_idx,nimages-1), 'parent', h(1).axes,'color','g','HorizontalAlignment','right','VerticalAlignment','top','FontSize',14,'FontWeight','bold','BackgroundColor','black', 'Interpreter','none');

    % usage instructions
    text(size(img,2)/2,size(img,1),sprintf('''SPACE'': Next Image  |  ''-'': Previous Image  |  ''x'': +100  |  ''y'': -100 | ''q'': quit'), 'parent', h(2).axes,'color','g','HorizontalAlignment','center','VerticalAlignment','bottom','FontSize',14,'FontWeight','bold', 'BackgroundColor','black');
end

repo: MSCKF


text:
function [y] = cellcrossproduct(x)

% [Y] = CELLSUBSELECT(X, BOOLVEC, DIM) outputs a cell-arry Y with the same dimension as X
% but for each input cell the n(n+1)/2 cross products of the rows are computed.
% the cross-products are ordered along the lower-triangle of the pairwise
% matrix.

nx = size(x);
if ~iscell(x) || length(nx)>2 || all(nx>1),
  error('incorrect input for cellcrossproduct');
end

n = cellfun('size', x, 1);
if ~all(n==n(1))
  error('each cell should have the same number of rows');
end
%[ix{1}(:,1),ix{1}(:,2)] = find(tril(ones(n(1))));

%y = cellfun(@crossprod, x, repmat(ix, nx), 'UniformOutput', 0);
y = cellfun(@crossprod, x, 'UniformOutput', 0);

function [y] = crossprod(x)

% FIXME only works in case of dim=2
n    = size(x);
indx = tril(ones(n(1)))==1;
y = zeros(0.5*n(1)*(n(1)+1), n(2));
for k = 1:n(2)
  tmp = tril(x(:,k)*x(:,k)');
  y(:,k) = tmp(indx);
end

% function [y] = crossprod(x, ix)
%
% %FIXME works only in case of dim=2
% n = size(x);
% y = zeros(0.5*n(1)*(n(1)+1), n(2));
% for k = 1:size(ix,1)
%   y(k,:) = x(ix(k,1),:).*x(ix(k,2),:);
% end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function pce_ode_hermite_test01 ( )

%*****************************************************************************80
%
%% PCE_ODE_HERMITE_TEST01 runs a test problem with PCE_ODE_HERMITE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PCE_ODE_HERMITE_TEST01:\n' );
  fprintf ( 1, '  Call PCE_ODE_HERMITE to compute a polynomial chaos expansion\n' );
  fprintf ( 1, '  for the ODE:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    u'' = - alpha * u,\n' );
  fprintf ( 1, '    u(0) = 1.\n' );

  ti = 0.0;
  tf = 2.0;
  nt = 200;
  ui = 1.0;
  np = 4;
  alpha_mu = 0.0;
  alpha_sigma = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial time         TI = %g\n', ti );
  fprintf ( 1, '  Final time           TF = %g\n', tf );
  fprintf ( 1, '  Number of time steps NT = %d\n', nt );
  fprintf ( 1, '  Initial condition    UI = %g\n', ui );
  fprintf ( 1, '  Expansion degree     NP = %d\n', np );
  fprintf ( 1, '  E(ALPHA)       ALPHA_MU = %g\n', alpha_mu );
  fprintf ( 1, '  STD(ALPHA)  ALPHA_SIGMA = %g\n', alpha_sigma );

  [ t, u ] = pce_ode_hermite ( ti, tf, nt, ui, np, alpha_mu, alpha_sigma );
%
%  Evaluate the exact expected value function.
%
  uex = ui * exp ( t.^2 / 2 );
%
%  Compare the first computed component against the exact expected value.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, ' i  T(i)  E(U(T(i)))    U(T(i),0)\n' );
  fprintf ( 1, '\n' );
  for i = 0 : 10 : nt
    fprintf ( 1, '  %4d  %6.3f  %14.6g  %14.6g  %14.6g\n', ...
      i, t(i+1), uex(i+1), u(i+1,1), abs ( uex(i+1) - u(i+1,1) ) );
  end
%
%  Plot the first computed component against the exact expected value function.
%
  figure ( 1 )
  plot ( t, u(:,1), t, uex )
  grid on
  xlabel ( '<--- T --->' )
  ylabel ( '<--- U --->' )
  title ( 'PCE Expected Value Versus Exact Expected Value' )

  filename = 'expected_value_plot.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Expected value comparison plotted in "%s".\n', filename );

  return
end


text:
function tinv = t_inv (x, v, mu, sigma)
% T_INV   Inverse of Student's T cumulative distribution function (cdf).
%    X=TINV(P,V,MU,SIGMA) returns the inverse of Student's T cdf with V degrees
%    of freedom, at the values in P, with mean MU and standard deviation
%    SIGMA.
%
%    The size of X is the common size of P and V. A scalar input
%    functions as a constant matrix of the same size as the other input.

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

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.
if nargin < 4
  sigma = 1;
end
if nargin < 3
  mu = 0;
end
if (~isscalar(mu) || ~isscalar(sigma))
  if ~((size(x,1) == size(mu,1)) && (size(mu,1) == size(sigma,1)))
    error ('norm_inv: x, mu and sigma must be of common size or scalars');
  end
end

x = (x-mu)./sigma;

if (nargin < 2)
  error('must give atleast 2 input arguments')
end

if (~isscalar(v))

  if (size(x,1) ~= size(v,1))
    error ('tinv: x and v must be of common size or scalar');
  end
end

tinv = zeros(size(x));

k = find( (x<0) || (x > 1) || isnan(x) || (v<0));
if (any(k))
  tinv(k) = NaN;
end

k = find ((x == 0) & (v > 0));
if (any (k))
  tinv(k) = -Inf;
end

k = find ((x == 1) & (v > 0));
if (any (k))
  tinv(k) = Inf;
end

k = find ((x > 0) & (x < 1) & (v > 0) & (v < 10000));
if (any (k))
  if (isscalar (v))
    tinv(k) = (sign (x(k) - 1/2) ...
    .* sqrt (v .* (1 ./ beta_inv (2*min(x(k),1-x(k)), ...
    v/2,1/2)-1)));
  else
    tinv(k) = (sign(x(k)-1/2) ...
    .* sqrt (v(k).*(1./beta_inv(2*min(x(k), 1-x(k)), ...
    v(k)/2,1/2)-1)));
  end
end

% For large v, use the quantiles of the standard normal
k = find ((x > 0) & (x < 1) & (v >= 10000));
if (any (k))
  tinv(k) = sqrt (2)*erfinv(2*x(k)-1);
end

end

% inv = sqrt (2) * erfinv (2 * x - 1);


text:
%%******************************************************************
%% infeaspt: generate an initial point for sdp.m
%%
%%  [X0,y0,Z0] = infeaspt(blk,At,C,b,options,scalefac);
%%
%%  options = 1  if want X0,Z0 to be scaled identity matrices
%%          = 2  if want X0,Z0 to be scalefac*(identity matrices).
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function [X0,y0,Z0] = infeaspt(blk,At,C,b,options,scalefac);
%%
   if (nargin < 5); options = 1; end;
   if (options == 1); scalefac = []; end;
   if (options == 2) & (nargin < 6); scalefac = 1000; end;
   if (scalefac <= 0); error('scalefac must a positive number'); end;
%%
   if ~iscell(At); At = {At}; end;
   if ~iscell(C);  C = {C}; end;
   m = length(b);
   if all(size(At) == [size(blk,1) m]);
      convertyes = zeros(size(blk,1),1);
      for p = 1:size(blk,1)
         if strcmp(blk{p,1},'s') & all(size(At{p,1}) == sum(blk{p,2}))
            convertyes(p) = 1;
         end
      end
      if any(convertyes)
         At = svec(blk,At,ones(size(blk,1),1));
      end
   end;
%%
   %%[blk,At,C,b] = validate(blk,At,C,b);
%%
   X0 = cell(size(C)); Z0 = cell(size(C));
   m = length(b);
   for p = 1:size(blk,1);
      pblk = blk(p,:);
      blktmp = pblk{2};
      n = length(C{p});
      y0 = zeros(m,1);
      b2 = 1 + abs(b');
      if (options == 1);
         if strcmp(pblk{1},'s');
            normAni = [];
            X0{p} = sparse(n,n); Z0{p} = sparse(n,n);
            ss = [0, cumsum(blktmp)];
            tt = [0, cumsum(blktmp.*(blktmp+1)/2)];
            for i = 1:length(pblk{2})
               if ~isempty(At{p,1})
                  pos = [tt(i)+1 : tt(i+1)];
                  Ai = At{p,1}(pos,:);
                  normAni = 1+sqrt(sum(Ai.*Ai));
               end
               if (length(At(p,:)) >= 2) %% for low rank constraints
                  dd = At{p,3};
                  qq = [0, cumsum(pblk{3})]; normtmp = ones(1,length(pblk{3}));
                  idxD = [0; find(diff(dd(:,1))); size(dd,1)];
                  for k = 1:length(pblk{3})
                      idx = [qq(k)+1 : qq(k+1)];
                      idx2 = [idxD(k)+1: idxD(k+1)];
                      Ak = At{p,2}(:,idx);
                      ii = dd(idx2,2)-qq(k); %% undo cumulative indexing
                      jj = dd(idx2,3)-qq(k);
                      len = pblk{3}(k);
                      Dk = spconvert([ii,jj,dd(idx2,4); len,len,0]);
                      tmp = Ak'*Ak*Dk;
                      normtmp(1,k) = 1+sqrt(sum(sum(tmp.*tmp')));
                  end
                  normAni = [normAni, normtmp];
               end
               pos = [ss(i)+1 : ss(i+1)];  ni = length(pos);
               tmp = C{p}(pos,pos);
               normCni = 1+sqrt(sum(sum(tmp.*tmp)));
               const  = 10; %%--- old: const = 1;
               constX = max([const,sqrt(ni),ni*(b2./normAni)]);
               constZ = max([const,sqrt(ni),normAni,normCni]);
               X0{p}(pos,pos) = constX*spdiags(1+1e-10*randmat(ni,1,0,'u'),0,ni,ni);
               Z0{p}(pos,pos) = constZ*spdiags(1+1e-10*randmat(ni,1,0,'u'),0,ni,ni);
            end
         elseif strcmp(pblk{1},'q');
            s = 1+[0, cumsum(blktmp)];
            len = length(blktmp);
            normC = 1+norm(C{p});
            normA = 1+sqrt(sum(At{p,1}.*At{p,1}));
            idenqX = zeros(sum(blktmp),1);
            idenqZ = zeros(sum(blktmp),1);
            idenqX(s(1:len)) = max([1,b2./normA])*sqrt(blktmp') ;
            idenqZ(s(1:len)) = max([sqrt(blktmp); max([normA,normC])*ones(1,len)])';
            idenqX(s(1:len)) = idenqX(s(1:len)).*(1+1e-10*randmat(len,1,0,'u'));
            idenqZ(s(1:len)) = idenqZ(s(1:len)).*(1+1e-10*randmat(len,1,0,'u'));
            X0{p} = idenqX;
            Z0{p} = idenqZ;
         elseif strcmp(pblk{1},'l');
            normC = 1+norm(C{p});
            normA = 1+sqrt(sum(At{p,1}.*At{p,1}));
            const = 10; %%--- old: const =1;
            constX = max([const,sqrt(n),sqrt(n)*b2./normA]);
            constZ = max([const,sqrt(n),normA,normC]);
            X0{p} = constX*(1+1e-10*randmat(n,1,0,'u'));
            Z0{p} = constZ*(1+1e-10*randmat(n,1,0,'u'));
         elseif strcmp(pblk{1},'u');
            X0{p} = sparse(n,1);
            Z0{p} = sparse(n,1);
         else
            error(' blk: some fields not specified correctly');
         end;
      elseif (options == 2);
         if strcmp(pblk{1},'s');
            n = sum(blktmp);
            X0{p} = scalefac*spdiags(1+1e-10*randmat(n,1,0,'u'),0,n,n);
            Z0{p} = scalefac*spdiags(1+1e-10*randmat(n,1,0,'u'),0,n,n);
         elseif strcmp(pblk{1},'q');
            s = 1+[0, cumsum(blktmp)];
            len = length(blktmp);
            idenq = zeros(sum(blktmp),1);
            idenq(s(1:len)) = 1+1e-10*randmat(len,1,0,'u');
            X0{p} = scalefac*idenq;
            Z0{p} = scalefac*idenq;
         elseif strcmp(pblk{1},'l');
            X0{p} = scalefac*(1+1e-10*randmat(n,1,0,'u'));
            Z0{p} = scalefac*(1+1e-10*randmat(n,1,0,'u'));
         elseif strcmp(pblk{1},'u');
            X0{p} = sparse(n,1);
            Z0{p} = sparse(n,1);
         else
            error(' blk: some fields not specified correctly');
         end
      end
   end
%%********************************************************************

author: yu-jiang

repo: radpbook

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


text:
function g = hsvargplvmLogLikeGradients(model)

%%%%%TEMP: to test the parallel function
   %g=hsvargplvmLogLikeGradientsPar(model);
   %return

%%%%%%

try
    pool_open = matlabpool('size')>0;
catch e
    pool_open = 0;
end

if pool_open && (isfield(model,'parallel') && model.parallel)
    % Run the parallel version of the current function.
    g=hsvargplvmLogLikeGradientsPar(model);
    % A 'dirty' trick to fix some parameters
    if isfield(model, 'fixParamIndices') && ~isempty(model.fixParamIndices)
        g(model.fixParamIndices) = 0;
    end
    return
end

if (model.H == 1 && isfield(model.layer{1}, 'dynamics') && ~isempty(model.layer{1}.dynamics))
    model.layer{1}.numModels = model.layer{1}.M;
    for m = 1:length(model.layer{1}.comp)
        model.layer{1}.comp{m}.vardist = model.layer{1}.vardist;
        model.layer{1}.comp{m}.dynamics = model.layer{1}.dynamics;
    end
    g = svargplvmLogLikeGradients(model.layer{1});
    %-----TMP
%     mm = model.layer{1}.comp{1};
%     mm.dynamics = model.layer{1}.dynamics;
%     mm.vardist = model.layer{1}.vardist;
%     pp = vargplvmExtractParam(mm);
%     mm = vargplvmExpandParam(mm, pp);
%     g = vargplvmLogLikeGradients(mm);
    %---
    % A 'dirty' trick to fix some parameters
    if isfield(model, 'fixParamIndices') && ~isempty(model.fixParamIndices)
        g(model.fixParamIndices) = 0;
    end
    return
end

g_leaves = hsvargplvmLogLikeGradientsLeaves(model.layer{1});

[g_nodes g_sharedLeaves] = hsvargplvmLogLikeGradientsNodes(model);
% Amend the vardistr. derivatives of the leaves only with those of the
% higher layer.
%%% TEMP %%%% !!
%g_leaves(1:model.layer{1}.vardist.nParams) = g_leaves(1:model.layer{1}.vardist.nParams).*0;
%%%%%%%%%%
g_leaves(1:model.layer{1}.vardist.nParams) = g_leaves(1:model.layer{1}.vardist.nParams) + g_sharedLeaves;

if model.H > 1 % NEW
    g_entropies = hsvargplvmLogLikeGradientsEntropies(model.layer{1});
else
    g_entropies = zeros(1,model.layer{1}.N*model.layer{1}.q);
end

% This is the gradient of the entropies (only affects covars). It's just -0.5*I
N = model.layer{1}.N;
Q1 = model.layer{1}.q;
%g_entropies = -0.5.*g_leaves(N*Q1+1:model.layer{1}.vardist.nParams);
g_leaves(N*Q1+1:model.layer{1}.vardist.nParams) = g_leaves(N*Q1+1:model.layer{1}.vardist.nParams) + g_entropies;

% Augment the derivatives of the vardistr. covars of each layer with the derivs.
% of the entropies (the leaves are done separately)
% Remember g_nodes has derivs for h=2 and upwards!!!
offset=1;
for h=2:model.H-1
    % g_entropies is the same for all models (constant) or changes only if
    % each layer has different Q.
    g_entropiesCur = hsvargplvmLogLikeGradientsEntropies(model.layer{h});
    % Find the indices for the var. covariances %%%%%%%%%% THIS IS WRONG
    NQ=N*model.layer{h}.q; % size of var. means (and var covars)
    % Inds of the vardistr.covars of the prev. node
    indStart = offset+NQ;
    indEnd = indStart + NQ - 1;
    g_nodes(indStart:indEnd)=g_nodes(indStart:indEnd)+g_entropiesCur;
    offset = offset + model.layer{h}.nParams;
end
g = [g_leaves  g_nodes];

%%%-- TEMP
%params = hsvargplvmExtractParam(model);
%g=zeros(1, length(params));
%g(1:length([g_leaves g_nodes])) = [g_leaves g_nodes];
%return
%%--

dynUsed = (isfield(model.layer{end},'dynamics') && ~isempty(model.layer{end}.dynamics));
% We only have to amend the parent's var. distr. derivatives if there are
% no dynamics, because then the derivs. are only calculated from the
% likelihood term. If there are dynamics, the derivs are calculated from
% the svargplvm function including the KL part.
if ~dynUsed
    % Amend the derivatives of the parent's var.distr. with the terms coming
    % from the KL
    g_parent = hsvargplvmLogLikeGradientsParent(model.layer{model.H});
    startInd = model.nParams-model.layer{model.H}.nParams+1;
    endInd = startInd + model.layer{model.H}.vardist.nParams-1;
    g(startInd:endInd) = g(startInd:endInd) + g_parent;
end

%--- NEW!! TEST
% If there are priors on parameters, add the contribution here
% %{
for h=1:model.H
    priorFound{h} = 0;
    for mm=1:model.layer{h}.M
        if isfield(model.layer{h}.comp{mm}, 'paramPriors') && ~isempty(model.layer{h}.comp{mm}.paramPriors)
            priorFound{h} = 1;
            break
        end
    end
end
c=1;
for h=1:model.H
    % TODO!!! Perhaps this M>1 should be forbidden if ANY layer has a
    % prior, check!!
    if priorFound{h}
        ind_cur = c:c+model.layer{h}.nParams-1;
        c=c+length(ind_cur);
        if model.layer{h}.M > 1
            error('Not implemented M > 1 yet!');
        end
        g(ind_cur) = g(ind_cur) + vargplvmParamPriorGradients(model.layer{h}.comp{1}, length(g(ind_cur)));
    end
end
% %}
%----

% A 'dirty' trick to fix some parameters
if isfield(model, 'fixParamIndices') && ~isempty(model.fixParamIndices)
    g(model.fixParamIndices) = 0;
end

end

%---------------  HELPER LOCAL FUNCTIONS ------------------------------%

% The derivatives for the leaves, ie only the likelihood part (no KL) of
% the X^1 -> Y bayesian gplvm model (svargplvm, actually) of layer 1.
% The gradients of the vardist part of this layer are incomplete, because
% they need to be amended with the derivatives of the higher layer. This
% extra derivatives are found in hsvargplvmLogLikeGradientNodes and the
% summation is done in the higher level of this function.
function g = hsvargplvmLogLikeGradientsLeaves(model)

g = [];
gShared = zeros(1, model.vardist.nParams);

for i=1:model.M
    model.comp{i}.vardist = model.vardist;
    model.comp{i}.onlyLikelihood = true;
    g_i = vargplvmLogLikeGradients(model.comp{i});
    % Now add the derivatives for the shared parameters (vardist)
    gShared = gShared + g_i(1:model.vardist.nParams);
    g_i = g_i((model.vardist.nParams+1):end);
    g = [g g_i];
end
g = [gShared g];
end

% g: The gradients for the "likelihood"+vardist part (ie no entropies, no
% KL) of all nodes
% gSharedLeaves: The gradients for the vardist part of the variational
% distribution of the first layer only, which has to be amended with the
% gradient obtained by hsvargplvmLogLikeGradientsLeaves. The vardist.
% gradients of the variational distributions of the other layers (apart from
% the parent) is handled internally (no need if H == 2)
function [g gSharedLeaves] = hsvargplvmLogLikeGradientsNodes(model)
g=[]; % Final derivative
gSharedLeavesMeans = zeros(model.layer{1}.N, model.layer{1}.q);
gSharedLeavesCovars = zeros(model.layer{1}.N, model.layer{1}.q);
for h=2:model.H
    dynUsed = false;
    if (h==model.H && (isfield(model.layer{h},'dynamics') && ~isempty(model.layer{h}.dynamics)));
        dynUsed = true;
    end

    g_sharedPrevNode = zeros(1,model.layer{h-1}.vardist.nParams);
    gSharedNodeMeans = zeros(model.layer{h-1}.N, model.layer{h-1}.q);
    gSharedNodeCovars = zeros(model.layer{h-1}.N, model.layer{h-1}.q);

    g_h = []; % Derivative wrt params of layer h

    % If we don't have dynamics, then we compute the gradients only of the
    % likelihood part, and then g_parent (in the main function of this
    % file) is calculated and amends with the KL part. However, if dynamics
    % kernel is used, then the gradients of the var. params are using the
    % likelihood and KL part in a mixed way and have to be computed
    % together. This is taken care of in the svargplvmLogLikeGradients
    % function, since the upper layer is, basically, an svargplvm. In that
    % case, we won't have to add the g_parent, as we will do everything
    % here at once.
    if ~dynUsed
        gShared = zeros(1, model.layer{h}.vardist.nParams); % Derivatives of the var. distr. of layer h
        for i=1:model.layer{h}.M % Derivative of the likelihood term of the m-th model of layer h
            model.layer{h}.comp{i}.vardist = model.layer{h}.vardist;
            model.layer{h}.comp{i}.onlyLikelihood = true;

            % Calculates the derivatives for the vardist. params and the other
            % params (apart from the KL) of the layer h > 1
            g_i = vargplvmLogLikeGradients(model.layer{h}.comp{i});

            % Now add the derivatives for the shared parameters (vardist)
            gShared = gShared + g_i(1:model.layer{h}.vardist.nParams);

            % g_i will now hold all the non-vardist. derivatives (the rest are
            % taken care of in gShared)
            g_i = g_i((model.layer{h}.vardist.nParams+1):end);

            g_h = [g_h g_i];
        end % for m...
        g_h = [gShared g_h];
    else
        model.layer{h}.numModels = model.layer{h}.M;
        for i=1:model.layer{h}.M
            model.layer{h}.comp{i}.vardist = model.layer{h}.vardist;
            model.layer{h}.comp{i}.dynamics = model.layer{h}.dynamics;
        end
        g_h = svargplvmLogLikeGradients(model.layer{h});
    end

    % Find the quantity that needs to be added to the PREVIOUS gShared.
    % This quantity is either the 'gSharedLeaves' which is returned and
    % added outside of this function (in the main body) or the
    % gSharedPrevNode , which is added in this function, after this for
    % loop finishes.
    for i=1:model.layer{h}.M
        means = model.layer{h-1}.vardist.means;
        covars = model.layer{h-1}.vardist.covars;

        %if model.centerMeans
        % TODO
        %end

        if ~isempty(model.layer{h}.comp{i}.latentIndices)
            latInd = model.layer{h}.comp{i}.latentIndices;
        else
            % All indices (M == 1)
            latInd = 1:model.layer{h-1}.vardist.latentDimension;
        end
        % In this layer h the "means" ie the X of layer h-1 are
        % here outputs. If we also have multOutput option, i.e. only
        % the full output space "means" will be grouped into
        % smaller subpsaces as defined in latentIndices. If we
        % don't have the multOutput>1 option, then means/covars will be the
        % original ones.
        means = means(:, latInd);
        covars = covars(:, latInd);

        %-- Amend the gShared OF THE PREVIOUS LAYER with the new terms due to the expectation
        % If multOutput(h) is set, then the expectation of <\tilde{F}>_{q(X_{h-1}} will
        % be split into Q expectations wrt q(X_{h-1}(:,q)), q=1,...,Q. In
        % general (not tested yet!!!) it will not be grouped by q, but
        % it could be defined as arbitrary subsets of X_{h-1}.
        beta = model.layer{h}.comp{i}.beta;
        if h == 2
            %-- Previous layer is leaves
            % Amend for the F3 term of the bound
            gSharedLeavesMeans(:, latInd) = beta * ...
                model.layer{h}.comp{i}.Z' * means;
            gSharedLeavesCovars(:, latInd) = 0.5*beta * ...
                repmat(diag(model.layer{h}.comp{i}.Z), 1, size(means,2));

            % Amend for the F0 term of the bound
            gSharedLeavesMeans(:, latInd) = gSharedLeavesMeans(:, latInd) - beta * means;
            gSharedLeavesCovars(:, latInd) = gSharedLeavesCovars(:, latInd) - 0.5*beta * ones(size(gSharedLeavesCovars(:, latInd)));

            % Reparametrization: from dF/dS -> dF/d log(S)
            gSharedLeavesCovars(:, latInd) = covars.*gSharedLeavesCovars(:, latInd);
        else % h > 2
           % error('H > 2 not implemented yet!!')
            % TODO!!!! For > 2 layers:
            % Previous layer is intermediate nodes
            % TODO: Here, we don't return the result but we added directly
            % to the previous iterations derivative (since everything
            % happens in the same function).

            % These computations are as above for the leaves case...

            % Amend for the F3 term of the bound
            gSharedNodeMeans(:, latInd) = beta * ...
                model.layer{h}.comp{i}.Z' * means;
            gSharedNodeCovars(:, latInd) = 0.5*beta * ...
                repmat(diag(model.layer{h}.comp{i}.Z), 1, size(means,2));

            % Amend for the F0 term of the bound
            gSharedNodeMeans(:, latInd) = gSharedNodeMeans(:, latInd) - beta * means;
            gSharedNodeCovars(:, latInd) = gSharedNodeCovars(:, latInd) - 0.5*beta * ones(size(gSharedNodeCovars(:, latInd)));

            % Reparametrization: from dF/dS -> dF/d log(S)
            gSharedNodeCovars(:, latInd) = covars.*gSharedNodeCovars(:, latInd);
            gSharedNodeVardist = [gSharedNodeMeans(:)' gSharedNodeCovars(:)'];

            g_sharedPrevNode(1:model.layer{h-1}.vardist.nParams) = ...
                    g_sharedPrevNode(1:model.layer{h-1}.vardist.nParams) + gSharedNodeVardist;

        end
        %--
    end

    % (for H > 2)
    if h > 2
        % The gShared of the PREVIOUS iteration needs to be augmented
           % Inds of the vardistr of the prev. node
        indStart = length(g) - model.layer{h-1}.nParams+1;
        indEnd = indStart + model.layer{h-1}.vardist.nParams - 1;
        g(indStart:indEnd) = g(indStart:indEnd) + g_sharedPrevNode;
    end

    % Now we can add the g_h of the current iteration
    g = [g g_h];
end % for h...
gSharedLeaves = [gSharedLeavesMeans(:)' gSharedLeavesCovars(:)'];
end

% The entropies of all variational distributions q(X^h) (this does not
% include the parent var. distr q(Z), for which we need to instead
% calculate a KL quantity.
function g = hsvargplvmLogLikeGradientsEntropies(model)
g=0.5*ones(1,model.N * model.q);
if isfield(model, 'DEBUG_entropy') && model.DEBUG_entropy
    g = -g;
end
end

% This calculates the derivatives of the KL quantity of the parent node
% It can be sped up by NOT calculating gKern, gInd etc, which are not used.
function gVar = hsvargplvmLogLikeGradientsParent(modelParent)

if isfield(modelParent, 'learnInducing')
    learnInducing = modelParent.learnInducing;
else
    learnInducing = true;
end

model = modelParent.comp{1};
model.vardist = modelParent.vardist;

%--- Now this is "copied" from vargplvmLogLikeGradients
gVarmeans = - model.vardist.means(:)';
% !!! the covars are optimized in the log space (otherwise the *
% becomes a / and the signs change)
gVarcovs = 0.5 - 0.5*model.vardist.covars(:)';

%{
% The following is not needed... parent node is not connected to any
% inducing points...!!!
if isfield(model, 'fixInducing') & model.fixInducing
    % TODO!!!
    warning('Implementation for fixing inducing points is not complete yet...')
    % Likelihood terms (coefficients)
    [gK_uu, gPsi0, gPsi1, gPsi2, g_Lambda, gBeta, tmpV] = vargpCovGrads(model);

    % Get (in three steps because the formula has three terms) the gradients of
    % the likelihood part w.r.t the data kernel parameters, variational means
    % and covariances (original ones). From the field model.vardist, only
    % vardist.means and vardist.covars and vardist.lantentDimension are used.
    [gKern1, gVarmeans1, gVarcovs1, gInd1] = kernVardistPsi1Gradient(model.kern, model.vardist, model.X_u, gPsi1',learnInducing);
    [gKern2, gVarmeans2, gVarcovs2, gInd2] = kernVardistPsi2Gradient(model.kern, model.vardist, model.X_u, gPsi2,learnInducing);
    [gKern0, gVarmeans0, gVarcovs0] = kernVardistPsi0Gradient(model.kern, model.vardist, gPsi0);

    %%% Compute Gradients with respect to X_u %%%
    gKX = kernGradX(model.kern, model.X_u, model.X_u);
    % The 2 accounts for the fact that covGrad is symmetric
    gKX = gKX*2;
    dgKX = kernDiagGradX(model.kern, model.X_u);
    for i = 1:model.k
        gKX(i, :, i) = dgKX(i, :);
    end
    % Allocate space for gX_u
    gX_u = zeros(model.k, model.q);
    % Compute portion associated with gK_u
    for i = 1:model.k
        for j = 1:model.q
            gX_u(i, j) = gKX(:, j, i)'*gK_uu(:, i);
        end
    end
    % This should work much faster
    %gX_u2 = kernKuuXuGradient(model.kern, model.X_u, gK_uu);
    gInd = gInd1 + gInd2 + gX_u(:)';
    %gVarmeans(model.inducingIndices, :) = gVarmeans(model.inducingIndices,
    %:) + gInd; % This should work AFTER reshaping the matrices...but here
    %we use all the indices anyway.
    gVarmeans = gVarmeans + gInd;
end
%}

gVar = [gVarmeans gVarcovs];
if isfield(model.vardist,'paramGroups')
    gVar = gVar*model.vardist.paramGroups;
end

end

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
% SP_PRECOMPUTE_PARAM: precompute all the fields, as in the space structure of the technical report, before mapping to the physical domain.
%
%     space = sp_precompute_param (space, msh);
%     space = sp_precompute_param (space, msh, 'option');
%
% INPUT:
%
%    space: object representing the discrete function space (see sp_scalar).
%    msh: mesh object containing the quadrature information (see msh_cartesian)
%   'option', value: additional optional parameters, currently available options are:
%
%              Name     |   Default value |  Meaning
%           ------------+-----------------+----------------------------------
%            value      |      true       |  compute shape_functions
%            gradient   |      false      |  compute shape_function_gradients
%            hessian    |      false      |  compute shape_function_hessians
%
% OUTPUT:
%
%    space: object containing the information of the input object, plus the
%            fields of the old structure, that are listed below. If no option
%            is given all the fields are computed. If an option is given,
%            only the selected fields will be computed.
%
%    FIELD_NAME      (SIZE)                             DESCRIPTION
%    ncomp           (scalar)                           number of components of the functions of the space (actually, 1)
%    ndof            (scalar)                           total number of degrees of freedom
%    ndof_dir        (1 x ndim vector)                  degrees of freedom along each direction
%    nsh_max         (scalar)                           maximum number of shape functions per element
%    nsh             (1 x msh.nel vector)               actual number of shape functions per each element
%    connectivity    (nsh_max x msh.nel vector)         indices of basis functions that do not vanish in each element
%    shape_functions (msh.nqn x nsh_max x msh.nel)      basis functions evaluated at each quadrature node in each element
%    shape_function_gradients
%          (ndim x msh.nqn x nsh_max x msh.nel)         basis function gradients evaluated at each quadrature node in each element
%    shape_function_hessians
%          (ndim x ndim x msh.nqn x nsh_max x msh.nel)  basis function hessians evaluated at each quadrature node in each element
%
% Copyright (C) 2009, 2010 Carlo de Falco
% Copyright (C) 2011, 2015, 2019 Rafael Vazquez
%
%    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.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function sp = sp_precompute_param (sp, msh, varargin)

  if (isempty (varargin))
    value = true;
    gradient = true;
    hessian = true;
  else
    if (~rem (length (varargin), 2) == 0)
      error ('sp_precompute: options must be passed in the [option, value] format');
    end
    value = false;
    gradient = false;
    hessian = false;
    for ii=1:2:length(varargin)-1
      if (strcmpi (varargin{ii}, 'value'))
        value = varargin{ii+1};
      elseif (strcmpi (varargin{ii}, 'gradient'))
        gradient = varargin{ii+1};
      elseif (strcmpi (varargin{ii}, 'hessian'))
        hessian = varargin{ii+1};
      else
        warning ('Ignoring unknown option %s', varargin {ii});
      end
    end
  end

  sp = struct (sp);

  sp_univ = sp.sp_univ;

  for idim = 1:msh.ndim
    elem_list{idim} = 1:msh.nel_dir(idim);
  end

  for idim = 1:msh.ndim
    nsh_dim{idim} = sp_univ(idim).nsh(elem_list{idim});
  end

  [nsh_grid{1:msh.ndim}] = ndgrid (nsh_dim{:});
  nsh = 1;
  for idim = 1:msh.ndim
    nsh = nsh .* nsh_grid{idim};
  end
  sp.nsh = nsh(:)';
  clear nsh

  for idim = 1:msh.ndim
    csize = ones (1, 2*msh.ndim);
    csize([idim, msh.ndim+idim]) = [sp_univ(idim).nsh_max, msh.nel_dir(idim)];
    crep = [sp_univ.nsh_max, msh.nel_dir];
    crep([idim, msh.ndim+idim]) = 1;

    conn{idim} = reshape (sp_univ(idim).connectivity(:,elem_list{idim}), csize);
    conn{idim} = repmat (conn{idim}, crep);
    conn{idim} = reshape (conn{idim}, [], msh.nel);
  end

  connectivity = zeros (sp.nsh_max, msh.nel);
  indices = ones (size (conn{1}));
  for idim = 1:msh.ndim
    indices = indices & conn{idim} ~= 0;
  end
  for idim = 1:msh.ndim
    conn{idim} = conn{idim}(indices);
  end
  connectivity(indices) = sub2ind ([sp.ndof_dir, 1], conn{:}); % The extra 1 makes things work in any dimension
  sp.connectivity = reshape (connectivity, sp.nsh_max, msh.nel);
  clear conn csize crep indices connectivity

  if (value || gradient || hessian)
    shp = cell(1,msh.ndim); shg = cell(1,msh.ndim); shh = cell(1,msh.ndim);
    for idim = 1:msh.ndim
      ssize = ones (1, 3*msh.ndim);
      ssize([idim, msh.ndim+idim, 2*msh.ndim+idim]) = [msh.nqn_dir(idim), sp_univ(idim).nsh_max, msh.nel_dir(idim)];
      srep = [msh.nqn_dir, sp_univ.nsh_max, msh.nel_dir];
      srep([idim, msh.ndim+idim, 2*msh.ndim+idim]) = 1;
      shp{idim} = reshape (sp_univ(idim).shape_functions(:,:,elem_list{idim}), ssize);
      shp{idim} = repmat (shp{idim}, srep);
      shp{idim} = reshape (shp{idim}, msh.nqn, sp.nsh_max, msh.nel);
      shg{idim} = reshape (sp_univ(idim).shape_function_gradients(:,:,elem_list{idim}), ssize);
      shg{idim} = repmat (shg{idim}, srep);
      shg{idim} = reshape (shg{idim}, msh.nqn, sp.nsh_max, msh.nel);
      shh{idim} = reshape (sp_univ(idim).shape_function_hessians(:,:,elem_list{idim}), ssize);
      shh{idim} = repmat (shh{idim}, srep);
      shh{idim} = reshape (shh{idim}, msh.nqn, sp.nsh_max, msh.nel);
    end

    if (value)
      sp.shape_functions = 1;
      for idim = 1:msh.ndim
        sp.shape_functions = sp.shape_functions .* shp{idim};
      end
    end

    if (gradient)
      for idim = 1:msh.ndim
        shape_fun_grad = shg{idim};
        for jdim = setdiff (1:msh.ndim, idim)
          shape_fun_grad = shape_fun_grad .* shp{jdim};
        end
        sp.shape_function_gradients(idim,:,:,:) = shape_fun_grad;
      end
      sp.shape_function_gradients = reshape (sp.shape_function_gradients, ...
                                    msh.ndim, msh.nqn, sp.nsh_max, msh.nel);
    end

% To be changed with isprop, as soon as classdef is working
    if (hessian && isfield (struct (msh), 'geo_map_der2'))
      for idim = 1:msh.ndim
        shape_fun_hess = shh{idim};
        for jdim = setdiff (1:msh.ndim, idim)
          shape_fun_hess = shape_fun_hess .* shp{jdim};
        end
        sp.shape_function_hessians(idim,idim,:,:,:) = shape_fun_hess;

        for jdim = setdiff (1:msh.ndim, idim)
          shape_fun_hess = shg{idim} .* shg{jdim};
          for kdim = setdiff (1:msh.ndim, [idim, jdim])
            shape_fun_hess = shape_fun_hess .* shp{kdim};
          end
          sp.shape_function_hessians(idim,jdim,:,:,:) = shape_fun_hess;
        end
      end
    end

    if (strcmpi (sp.space_type, 'NURBS'))
      sp = bsp_2_nrb__ (sp, msh, sp.weights);
    end
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
%GM_PHD_Predict_Birth
%Last modified 2nd September 2013
%Matlab code by Bryan Clarke b.clarke@acfr.usyd.edu.au

%This file performs prediction for newly birthed and spawned targets
%This is necessary since their positions were initialised in a previous timestep and they may
%have moved between then and now
%We iterate through j = 1..J_b,k where J_b,k is number of birthed landmarks
%at time k (but the landmarks correspond to the measurement from time k-1)
%For this implementation, birthed and spawned targets are identical except
%they have weights that are calculated from different functions, and
%different starting covariances. A target is spawned or birthed depending
%on which weighting function will give it a higher starting weight.

%The means of these will be the position in Cartesian space where they are detected
%The covariances & weights are calculated according to Vo&Ma
s = sprintf('Step 1: Prediction for birthed and spawned targets.');
disp(s);

m_birth_before_prediction = [m_birth, m_spawn]; %Need to store these BEFORE prediction for use in the update step.

%Perform prediction for birthed targets using birthed velocities.
for j = 1:numBirthedTargets
    i = i + 1;
    %w_birth was already instantiated in GM_PHD_Create_Birth
    m_birth(:,j) = F * m_birth(:,j);
    P_range = calculateDataRange4(j);
    P_birth(:,P_range) = Q + F * P_birth(:,P_range) * F';
end
%Perform prediction for birthed targets using birthed velocities.
for j = 1:numSpawnedTargets
    i = i + 1;
    %w_birth was already instantiated in GM_PHD_Create_Birth
    m_spawn(:,j) = F * m_spawn(:,j);
    P_range = calculateDataRange4(j);
    P_spawn(:,P_range) = Q + F * P_spawn(:,P_range) * F';
end

if(VERBOSE == 1)
  for j = 1:numBirthedTargets
        thisM = m_birth(:,j);
        s = sprintf('Birthed target %d: %3.4f %3.4f %3.4f %3.4f Weight %3.9f', j, thisM(1), thisM(2), thisM(3), thisM(4), w_birth(j));
        disp(s);
  end
  for j = 1:numSpawnedTargets
        thisM = m_spawn(:,j);
        s = sprintf('Spawned target %d: %3.4f %3.4f %3.4f %3.4f Weight %3.9f', j, thisM(1), thisM(2), thisM(3), thisM(4), w_birth(j));
        disp(s);
   end
end


text:
function v = repmat(v,varargin)
% overloads repmat

v.x = repmat(v.x,varargin{:});
v.y = repmat(v.y,varargin{:});
v.z = repmat(v.z,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
clear all
close all
clc

%% get sie data

sie2mat('example.sie','video',1 );

%% save video data

%convert2vid( SIE_VIDEO );


text:
% [INPUT]
% data = A float t-by-n matrix (-Inf,Inf) representing the model input.
% tail = A string representing the target tail:
%   - 'L' for lower tail;
%   - 'U' for upper tail.
% bw = An integer [21,252] representing the dimension of each rolling window (optional, default=252).
% f = A float [0.05,0.20] representing the percentage of observations to be included in tails (optional, default=0.10).
% pt = A float [0,1) representing the initial penantly term for underrepresented samples with respect to the bandwidth (optional, default=0.5).
%
% [OUTPUT]
% chi = A float n-by-n-by-t matrix [0,1] representing the Chi coefficients.
% chi_bar = A float n-by-n-by-t matrix [-1,1] representing the Chi Bar coefficients.
%
% [NOTES]
% The bandwidth is automatically expanded as much as possible following an optimality criteria.

function [chi,chi_bar] = asymptotic_tail_dependence(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('data',@(x)validateattributes(x,{'double'},{'real' '2d' 'nonempty'}));
        ip.addRequired('tail',@(x)any(validatestring(x,{'L' 'U'})));
        ip.addOptional('bw',252,@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 21 '<=' 252 'scalar'}));
        ip.addOptional('f',0.10,@(x)validateattributes(x,{'double'},{'real' 'finite' '>=' 0.05 '<=' 0.2 'scalar'}));
        ip.addOptional('pt',0.5,@(x)validateattributes(x,{'double'},{'real' 'finite' '>=' 0 '<' 1 'scalar'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    [data,bw] = validate_input(ipr.data,ipr.bw);
    tail = ipr.tail;
    f = ipr.f;
    pt = ipr.pt;

    nargoutchk(1,2);

    [chi,chi_bar] = asymptotic_tail_dependence_internal(data,bw,tail,f,pt);

end

function [chi,chi_bar] = asymptotic_tail_dependence_internal(data,bw,tail,f,pt)

    up = isempty(getCurrentTask());

    [t,n] = size(data);

    c = nchoosek(1:n,2);
    c_len = size(c,1);

    dc = cell(c_len,2);

    for i = 1:c_len
        c_i = c(i,:);
        dc(i,:) = {c_i data(:,c_i)};
    end

    pt = [linspace(pt,1,bw).'; ones(t-bw,1)];

    dc_results = cell(c_len,2);

    if (up)
        parfor k = 1:c_len
            windows = extract_rolling_windows(dc{k,2},bw);

            chi_k = zeros(t,1);
            chibar_k = zeros(t,1);

            for w = 1:t
                dc_kw = windows{w};

                [chi_kw,chibar_kw] = calculate_coefficients(dc_kw,tail,f);
                chi_k(w) = chi_kw;
                chibar_k(w) = chibar_kw;
            end

            dc_results(k,:) = {(chi_k .* pt) (chibar_k .* pt)};
        end
    else
        for k = 1:c_len
            windows = extract_rolling_windows(dc{k,2},bw);

            chi_k = zeros(t,1);
            chibar_k = zeros(t,1);

            for w = 1:t
                dc_kw = windows{w};

                [chi_kw,chibar_kw] = calculate_coefficients(dc_kw,tail,f);
                chi_k(w) = chi_kw;
                chibar_k(w) = chibar_kw;
            end

            dc_results(k,:) = {(chi_k .* pt) (chibar_k .* pt)};
        end
    end

    en = eye(n);
    chi = repmat(en,1,1,t);
    chi_bar = repmat(en,1,1,t);

    for k = 1:c_len
        dc_k = dc{k,1};
        i = dc_k(1);
        j = dc_k(2);

        [chi_k,chibar_k] = deal(dc_results{k,:});
        chi(i,j,:) = chi_k;
        chi(j,i,:) = chi_k;
        chi_bar(i,j,:) = chibar_k;
        chi_bar(j,i,:) = chibar_k;
    end

    for i = 1:t
        nan_indices = isnan(data(i,:));

        if (any(nan_indices))
            chi(nan_indices,:) = NaN;
            chi(:,nan_indices) = NaN;

            chi_bar(nan_indices,:) = NaN;
            chi_bar(:,nan_indices) = NaN;
        end
    end

end

function [chi,chi_bar] = calculate_coefficients(data,tail,f)

    if (any(any(isnan(data),1)))
        chi = NaN;
        chi_bar = NaN;
        return;
    end

    t = size(data,1);
    t1 = t + 1;

    nu = max(round(t * f),1);
    nu1 = nu + 1;

    if (strcmp(tail,'L'))
        u1 = 1 - (tiedrank(data(:,1)) ./ t1);
        u2 = 1 - (tiedrank(data(:,2)) ./ t1);
    else
        u1 = tiedrank(data(:,1)) ./ t1;
        u2 = tiedrank(data(:,2)) ./ t1;
    end

    zs = -1 ./ log(u1);
    zt = -1 ./ log(u2);
    z =  sort(min(zs,zt),'descend');

    if (nu == 1)
        eta = 0;
    else
        eta = sum(log(z(1:nu) ./ z(nu1)));
    end

    chi_bar = ((2 / nu1) * eta) - 1;
    sigma = (chi_bar + 1)^2 / nu1;
    ci = chi_bar + (norminv(0.975) * sigma^0.5);

    if (ci >= 1)
        chi = z(nu1) * (nu1 / t);
    else
       chi = 0;
    end

end

function [data,bw] = validate_input(data,bw)

    [t,n] = size(data);

    if (n < 2)
        error('The value of ''data'' is invalid. Expected input to be a matrix with at least 2 columns.');
    end

    nan_counts = sum(isnan(data),1);
    nan_threshold = round(t * 0.70,0);

    if (any(nan_counts > nan_threshold))
        error(['The value of ''data'' is invalid. Expected input to contain no more than 70% of NaN values (' num2str(nan_threshold) ') for each time series.']);
    end

    for i = 2:ceil(t / bw)
        bw_i = bw * i;

        if ((bw_i / t) >= 0.3)
            bw = bw_i;
            break;
        end
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function D = addoccludedtag(D, words)
%
% Adds the field object.crop (it also considers an occluded object as a
% cropped object).
%
% Removes the words 'occluded' from the object name.

if nargin < 2
    words = {'occluded'};
end

Nwords = length(words);
changes = 0;

for i = 1:length(D)
    if isfield(D(i).annotation, 'object')
        Nobjects = length(D(i).annotation.object);
        for j = 1:Nobjects

%             if isfield(D(i).annotation.object(j), 'occluded')
%                 D(i).annotation.object(j).occluded == 'yes'
%             end

            w = lower(getwords(D(i).annotation.object(j).name));
            k = ismember(w, words);

            if sum(k)>0
                D(i).annotation.object(j).occluded = 'yes';
                changes = changes + 1;
                % remove words
                for n = 1:Nwords
                    D(i).annotation.object(j).name = strrep(D(i).annotation.object(j).name, words{n}, '');
                end
            end
        end
    end
end

% Make sure there are no new empty object names
D = LMvalidobjects(D);

fprintf('occluded field: %d changes to ''yes''\n', changes')

author: CSAILVision

repo: LabelMeToolbox


text:
function res = isrelaxfeasible(F)

% Check if solution avaliable
currsol = evalin('caller','sdpvar(''getSolution'')');
if isempty(currsol)
    disp('No solution available.')
    return
end

F = flatten(F);

nlmi = length(F.LMIid);
spaces = ['                                    '];
if (nlmi == 0) & (neq == 0)
    disp('empty LMI')
    return
end

lmiinfo{1} = 'LMI';
lmiinfo{2} = 'Element-wise';
lmiinfo{3} = 'Equality constraint';
lmiinfo{4} = 'Second order cone constraint';
lmiinfo{5} = 'Rotated Lorentz constraint';

header = {'ID','Constraint','Type','Residual (should be > 0)','Tag'};

if nlmi>0
    for j = 1:nlmi
        F0 = relaxdouble(F.clauses{j}.data);
        if any(isnan(F0(:)))
            res = NaN;
        else
            switch F.clauses{j}.type
                case 1
                    res = min(eig(F0));
                case 2
                    res = min(min(F0));
                case 3
                    res = -max(max(abs(F0)));
                case 4
                    res = F0(1)-norm(F0(2:end));
                case 5
                    res = 2*F0(1)*F0(2)-norm(F0(3:end))^2
            end
        end
        data{j,1} = ['#' num2str(j)];
        data{j,2} = F.clauses{j}.symbolic;
        data{j,3} = lmiinfo{F.clauses{j}.type};
        data{j,4} = res;
        data{j,5} = F.clauses{j}.handle;
    end
end

res = [data{:,4}];

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [ y, symm ] = cvx_s_sparse( m, n, symm, i, j )

%CVX_S_SPARSE Matrices with a fixed sparsity pattern.

if nargin < 5,
    error( 'Sparsity structure missing.' );
elseif ~isnumeric( i ) || ~isnumeric( j ),
    error( 'Sparsity arguments must be vectors of nonnegative integers.' );
elseif any( i <= 0 ) || any( j <= 0 ) || any( i ~= floor( i ) ) || any( j ~= floor( j ) ),
    error( 'Sparsity arguments must be vectors nonnegative integers.' );
elseif numel( i ) ~= 1 && numel( j ) ~= 1 && numel( i ) ~= numel( j ),
    error( 'Sparsity arguments have incompatible size.' );
elseif any( i > m ) || any( j > n ),
    error( 'One or more indices are out of range.' );
elseif symm && m ~= n,
    error( 'Symmetric structure requires a square matrix.' );
end
i = i(:);
j = j(:);
nz = max( numel(i), numel(j) );
if symm,
    t = max(i,j);
    j = min(i,j);
    i = t;
end
[ c, cndxs ] = sort( i + ( j - 1 ) * m );
tt = [true;c(2:end)~=c(1:end-1)];
c = c(tt);
r = 1 : length(c);
if symm,
    if numel(i) > 1,
        i = i(cndxs(tt));
    end
    if numel(j) > 1,
        j = j(cndxs(tt));
    end
    tt = i ~= j;
    r = [ r , r(tt) ];
    c = [ c ; j(tt) + ( i(tt) - 1 ) * m ];
end
y = min( sparse( r, c, 1, nz, m * n ), 1 );
symm = false;

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [mPerc, mProb] =calc_zlta(mCat0,mCat1,mCat2,fTstart,fT,fTw,nTbin,nN)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example: mLTA=calc_zlta(mCat00,mCat20,params.fTstart, fT,fTw,nTbin, nN);
%
% This function calculates the rate changes (z-value) of earthquake
% occurrencs between two periods. This function calculates rate changes for
% all the grid nodes together. Input is either a single vector or a whole
% matrix columnswise only with dates (not the whole catalog). Output is the
% z(lta)- and its probability value for each grid point.
%
% Author: van Stiphout, Thomas
% Email: vanstiphout@sed.ethz.ch
% Created: 7. Aug. 2007
% Changed: 14. Aug.2007
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Variables:
% mCat0          Catalog complete (only origin time)
% mCat1          Catalog period 1 (vector or matrix with yrs in column)
% mCat2          Catalog period 2 (vector or matrix with yrs in column)
% fTstart        Begin of time period 1
% fT             Date for which rate change is calculated
% fTw            Length of Time window of second period
% nTbin          Length of Time steps for histogram
% nN             Sampling volume
%
% Output:
% mLTA           Scalar or vector with beta values for each input column
% mProb          Scalar or vector with probability for beta values for each
%                input column. The probability is either calculated based
%                on a synthetic catalog or a real-like catalog.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% disp('~/zmap/src/thomas/seismicrates/calc_zlta.m');
% probability calculation based on synthetic catalogs.
% Synthetic catalogs either based catalog with uniform rates (o) or according to
% complete catalog (1)
bSyn=logical(0);

% calculate histogram for different time periods
% nSteps=floor((fT-fTstart)/(nTbin/365));
% vR1=histc(mCat1,  fTstart:(nTbin/365):fTstart+nSteps*(nTbin/365)  );
vR1=histc(mCat1,fTstart : nTbin/365 : fT-fTw,1);
vR1=vR1(1:end-1,:);

% vR1=histc(mCat1,fTimeStart:fTimeSteps/365:fTimeCut+fTimeWindow);
vR2=histc(mCat2,fT-fTw : nTbin/365:fT,1);
vR2=vR2(1:end-1,:);

% calculate the mean rate for different periods
mean1=mean(vR1);
mean2=mean(vR2);

var1 = var(vR1);
var2 = var(vR2);

if isempty(vR1)
    disp('Warning - Time Period 1 is without any event');
elseif isempty(vR2)
    disp('Warning - Time Period 2 is without any event');
end

% create synthetic catalogs to extimate significance level
% reset random number generator
rand('state',sum(100*clock));
if bSyn
    vPos=ceil(rand(nN,1000).*size(mCat0,1));
    mSyn1=mCat0(vPos);
 else
    mSyn1=rand(nN,1000)*(fT-fTstart)+fTstart;
end
% apply histogram to synthetic catalog
vS1=histc(mSyn1,fTstart : nTbin/365 : fT-fTw,1);
vS1=vS1(1:end-1,:);
vS2=histc(mSyn1,fT-fTw : nTbin/365:fT,1);
vS2=vS2(1:end-1,:);
% calculate the mean rate for different periods in synthetic catalog
mSynMean1=mean(vS1);
mSynMean2=mean(vS2);
% calculate z(lta) values for synthetic catalog
mSynPerc=(mSynMean2./mSynMean1.*100)-100;
% calculate values for normal distribution
% [jbt(i), jbp(i)]=jbtest(mSynLTA);
% llt(i)=lillietest(mSynLTA);
[mu,s] = normfit(mSynPerc);

% z(lta)
mPerc=(mean2./mean1.*100)-100;
% calculate the probability of z(lta)-values
[mProb] = 1-normpdf(mPerc,mu,s);
% figure;plot(mLTA,mProb,'.');

author: CelsoReyes


text:
function [N,x0] = affine_null_space(A,b,varargin)
  % Given a system Ax = b, determine a matrix N spanning the right null space
  % of A and a feasible solution x0 so that:
  %
  %     A * (N * y + x0) = b  for any y
  %
  % Inputs:
  %   A  m by n (sparse) matrix.
  %   b  m by #b right-hand side
  %   Options:
  %     'Tol'  followed by tolerance for determine rank (what's considered
  %       zero?)
  %     'Method'  followed by either:
  %        {'qr'}  use QR decomposition of A' (robust, best understood, slowest)
  %        'luq'  use LUQ decomposition
  %        'rq'  use QR decomposition of A (good when m << n)
  %        'svd'  use SVD decompostion (only use for small/dense matrices)
  %        'rrlu'  simplifed version of LUQ (**broken**)
  % Outputs:
  %   N  n by #N matrix spanning null space, where #N = m - rowrank(A)
  %   x0 n by #b, so that columns are feasible solutions
  %

  tol = [];
  method = 'qr';
  % default values
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Method','Tol'},{'method','tol'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end
  if nargin<2 || isempty(b)
    b = zeros(size(A,1),1);
  end

  if isempty(tol)
    tol = max(max(size(A)) * norm(A,1) * eps,100*eps);
  end

  switch method
  case 'luq'
    %% Special sparse LUQ decomposition
    %tol
    %[L,U,Q] = luq(A,1,tol);
    %% Rank
    %nc = find(any(abs(U)>tol,2),1,'last');
    %nc
    %if isempty(nc)
    %  nc = 0;
    %  m = size(A,1)-nc;
    %  if nargout>=2
    %    x0 = ones(m,1);
    %  end
    %  N = speye(size(A,1),m);
    %else
    %  if nargout>=2
    %    y0 = U(1:nc,1:nc)\(speye(nc,size(L,1))*(L\b));
    %    x0 = Q\[y0;zeros(size(Q,1)-nc,size(b,2))];
    %  end
    %  QQ = Q^-1;
    %  N = QQ(:,nc+1:end);
    %end
    %%big = max(abs(N));
    %%if big > 0
    %%  N = N/big;
    %%end

    % FAR MORE ACCURATE THAN luq(A,...); see doc of spspaces (L is more accurate
    % than Q)
    [L,U,Q] = luq(A',0,tol);
    % this looks gross but somehow doesn't destroy sparsity
    LL = L^-1;
    S = max(abs(U),[],2);
    if ~isempty(S)
      J = find(S<=tol);
    else
      J = (1:size(S,1))';
    end
    N = LL(J,:)';
    nc = find(any(abs(U)>tol,2),1,'last');
    Unc = U(1:nc,1:nc);
    Lnc = L(:,1:nc);
    %x0 = A'*(((A')'*A')\b);
    %x0 = L*U*Q*((Q'*U'*L'*L*U*Q)\b);
    %x0 = L*U*(((U'*L'*L*U)\(Q'\b)));
    %x0 = Lnc*Unc*(((Unc'*Lnc'*Lnc*Unc)\(Q'\b)));
    %x0 = Lnc*(((Lnc'*Lnc)\(Unc'\(Q'\b))));
    x0 = LL'*(speye(size(L,2),size(Unc,1))*(Unc'\(Q'\b)));
  case 'rrlu'
    % "Strong rank revealing LU factorizations" [Miranian & Gu 2002]
    % https://math.berkeley.edu/~luiza/RRLU.pdf
    %
    % Seems there is a typo. They write: "Nr = [-A11\A12;I_{m-k,n-k}]" But so
    % that A*Nr makes sense, I think it should be "Nr = [-A11\A12;I_{n-k,n-k}]"
    %
    %
    m = size(A,1);
    n = size(A,2);
    [L,U,P,Q] = lu(sparse(A),tol);

    NZ = find((any(abs(U)>tol,2)));
    % LUQ just checks the diagonal:
    NZluq = find(abs(diag(U))>tol);
    if ~isempty(setxor(NZ,NZluq))
      size(NZ)
      size(NZluq)
      warning('Not handling non-zero bottom right corner (use luq)');
    end

    Z = find(~(any(abs(U)>tol,2)));
    R = sparse((1:size(U,1))',[NZ;Z],1);

    C = blkdiag(R',speye(size(U,2)-size(R,2)));
    % L U               = P A Q
    % L R' R U C C'     = P A Q
    % L R' R U C        = P A Q C
    % L R' [U11 U12; 0] = P A Q C
    %
    % let M = [-U11\U12;I]
    %
    % [U11 U12] M = 0
    %
    % L R' [U11 U12] M = P A Q C M
    % L R' 0           = P A Q C M
    %
    % Implies
    %
    % N = Q C M --> A N = 0
    %
    UU = R*U*C;
    nc = sum((any(abs(UU)>tol,2)));
    U11 = UU(1:nc,1:nc);
    U12 = UU(1:nc,nc+1:end);
    M = [-U11\U12;speye(size(U12,2),size(U12,2))];
    N = Q*C*M;

    % We have:
    %
    % L R' [U11 U12; 0] = P A Q C
    %
    % A x = b
    % Let x = Q C y
    % P A Q C y = P b
    % L R' [U11 U12; 0] y = P b
    % Assume: (L R') is invertible:
    % [U11 U12; 0] y = (L R')\(P b)
    % [U11 U12; 0] y = [b1;b2];
    % Invertible iff b2 = 0
    % [U11 U12; 0] y = [b1;0];
    % y = [U11\b1;y2] for any y2
    % And might as well set y2 to 0
    % y = [U11\b1; 0]
    %
    LL = L*R';
    bb = LL\(P*b);
    y = [U11\(bb(1:nc,:));zeros(size(U,2)-nc,size(bb,2))];
    x0 = Q*C*y;

  case 'qr'
    [Q,R,E] = qr(A');
    % Rank of A
    nc = find(any(abs(R)>tol,2),1,'last');
    % Q = [Q₁,N]
    Q1 = Q(:,1:nc);
    % A possibly non-unique solution
    if nargout>=2
      x0 = Q1*(R(1:nc,1:nc)'\(E(:,1:nc)'*(b)));
    end
    N = Q(:,nc+1:end);
  case 'rq'
    % http://mathoverflow.net/a/253997/23064
    [Q,R,E] = qr(A);
    nc = find(any(abs(R)>tol,2),1,'last');
    R1 = R(1:nc,1:nc);
    R2 = R(1:nc,nc+1:end);
    n = size(A,2);
    N = E*[-(R1\R2);speye(n-nc,n-nc)];
    %assert(nargout <= 1 && 'x0 not supported for rq');
    b1 = Q(:,1:nc)'*b;
    x0 = E*[R1\b1;zeros(size(E,2)-size(R1,1),size(b1,2))];
  case 'svd'
    [U,S,V] = svd(full(A));
    % Carefully extract diagonal of S
    Sdiag = S(sub2ind(size(S),1:min(size(S)),1:min(size(S))));
    Z = abs(Sdiag)<tol;
    N = V(:,setdiff(1:end,find(~Z)));
    Sdiag(Z) = 0;
    Sdiag(~Z) = 1./Sdiag(~Z);
    % Place back into S without changing size of S
    S(sub2ind(size(S),1:min(size(S)),1:min(size(S)))) = Sdiag;
    Apinv = V*S'*U';
    x0 = Apinv * b;
  end
  % Zap anything below tolerance
  %N(abs(N) < tol) = 0
  [NI,NJ,NV] = find(N);
  N = sparse(NI,NJ,(abs(NV)>tol).*NV,size(N,1),size(N,2));

  %assert(max(abs(A*(N*rand(size(N,2),size(b,2)) + x0) - b)) < 1e-10, ...
  %  'Should span solutions to A x = b');
  if nargout>1 && ~(all(all(abs(A*x0-b)<tol)))
    % Should check that constraint right-hand sides are compatible:
    %   [Q,R,E] = qr(A');
    %   rank_A = find(any(abs(R)>tol,2),1,'last');
    %   [Q,R,E] = qr([A b]');
    %   rank_Ab = find(any(abs(R)>tol,2),1,'last');
    %   assert(rank_Ab <= rank_A);
    %
    warning('MATLAB:singularMatrix', ...
      'A*x0 ~= b, this may mean that  A*x = b is impossible to satisfy');
  end
end

author: alecjacobson

repo: gptoolbox


text:
function [kg] = ug2kg(ug)
% Convert mass from micrograms to kilograms.
% Chad Greene 2012
kg = ug*1E-9 ;

author: Sable


text:
function [Rot,v,x,PosAmers,P] = iekfPropagation(dt,Rot,v,x,omega_b,a_b,...
    PosAmers,P,omega,acc,Q,g)
NbAmers = size(PosAmers,2);

% state propagation
Rot = Rot*expSO3((omega-omega_b)*dt);
v = v+(Rot*(acc-a_b)+g)*dt;
x = x+v*dt;

% covariance propagation
F = eye(size(P));
F(4:6,1:3) = vecto(g)*dt;
F(7:9,1:3) = vecto(g)*dt*dt;
F(7:9,4:6) = eye(3)*dt;
F(1:3,10:12) = -Rot*dt;
F(4:6,10:12) = -vecto(v)*Rot*dt;
F(7:9,10:12) = -vecto(x)*Rot*dt;
for i = 1:NbAmers
    posAmers_i = PosAmers(:,i);
    F(13+3*i:15+3*i,10:12) = -vecto(posAmers_i)*Rot*dt;
end
F(4:6,13:15) = -Rot*dt;
F(7:9,13:15) = -Rot*dt*dt;

G = zeros(size(P,1),size(Q,1));
G(1:3,1:3) = Rot;
G(4:6,1:3) = vecto(v)*Rot;
G(7:9,1:3) = vecto(x)*Rot;
G(4:6,4:6) = Rot;
G(7:9,4:6) = Rot*dt;
G(10:15,7:12) = eye(6);
G(1:3,7:9) = Rot*dt;
G(4:6,7:9) = vecto(v)*Rot*dt*dt;
G(7:9,7:9) = vecto(x)*Rot*dt*dt*dt;
G(4:6,10:12) = Rot*dt*dt;
G(7:9,10:12) = Rot*dt*dt*dt;

P = F*P*F' + G*(Q*dt)*G'*dt;
end

author: mbrossar


text:
function ShowCmap(M, h)
% ShowCmap(M, h)
% a function to display and examine a colormap
% M (nx3): color map
% h: figure handle
%
% see also readXcol, ROIcmap, rgbdectohex, and ScaleToMap
%
   figure(h); clf
   subplot (211);
   colormap(M);
   image ([1:1:size(M,1)]);
   str = sprintf('%d colors color map. Pick colors with mouse\nHit "enter" to quit', size(M,1));
   title(str, 'fontsize',14);
   drawnow;

   i = 0;
   subplot (269);cla

   x1 = 1;
   while (~isempty(x1)),
      [x1,y] = ginput (1);
      if (~isempty(x1)),
         x1 = floor(x1(1)-0.5)+1;
         subplot (269);
         addsquare([0 i], [2.5 1+i], M(x1,:)); hold on
         plot (-0.2, i+0.5, 'k*');
         axis ([-1 3 -1 11]);
         str = sprintf ('Col %d: %.3g %.3g %.3g', x1, M(x1,1), M(x1,2), M(x1,3));
         ht = text (3, 0.3+i, 0, str, 'fontsize',14, 'color', M(x1,:));
         title (str);
         i = rem(i +1, 10);
      end
   end
return;