text:
function [features, gparams, feature_info] = init_features(features, gparams, is_color_image, img_sample_sz, size_mode)

if nargin < 3
    size_mode = 'same';
end

% Set missing global parameters to default values
if ~isfield(gparams, 'normalize_power')
    gparams.normalize_power = [];
end
if ~isfield(gparams, 'normalize_size')
    gparams.normalize_size = true;
end
if ~isfield(gparams, 'normalize_dim')
    gparams.normalize_dim = false;
end
if ~isfield(gparams, 'square_root_normalization')
    gparams.square_root_normalization = false;
end
if ~isfield(gparams, 'use_gpu')
    gparams.use_gpu = false;
end

% find which features to keep
feat_ind = false(length(features),1);
for n = 1:length(features)

    if ~isfield(features{n}.fparams,'useForColor')
        features{n}.fparams.useForColor = true;
    end

    if ~isfield(features{n}.fparams,'useForGray')
        features{n}.fparams.useForGray = true;
    end

    if (features{n}.fparams.useForColor && is_color_image) || (features{n}.fparams.useForGray && ~is_color_image)
        % keep feature
        feat_ind(n) = true;
    end
end

% remove features that are not used
features = features(feat_ind);

num_features = length(features);

feature_info.min_cell_size = zeros(num_features,1);

% Initialize features by
% - setting the dimension (nDim)
% - specifying if a cell array is returned (is_cell)
% - setting default values of missing feature-specific parameters
% - loading and initializing necessary data (e.g. the lookup table or the network)
for k = 1:length(features)
    if isequal(features{k}.getFeature, @get_fhog)
        if ~isfield(features{k}.fparams, 'nOrients')
            features{k}.fparams.nOrients = 9;
        end
        features{k}.fparams.nDim = 3*features{k}.fparams.nOrients+5-1;
        features{k}.is_cell = false;
        features{k}.is_cnn = false;

    elseif isequal(features{k}.getFeature, @get_table_feature)
        table = load(['lookup_tables/' features{k}.fparams.tablename]);
        features{k}.fparams.nDim = size(table.(features{k}.fparams.tablename),2);
        features{k}.is_cell = false;
        features{k}.is_cnn = false;

    elseif isequal(features{k}.getFeature, @get_colorspace)
        features{k}.fparams.nDim = 1;
        features{k}.is_cell = false;
        features{k}.is_cnn = false;

    elseif isequal(features{k}.getFeature, @get_cnn_layers) || isequal(features{k}.getFeature, @get_OFcnn_layers)
        % make sure the layers are correcly sorted
        features{k}.fparams.output_layer = sort(features{k}.fparams.output_layer);

        % Set default parameters
        if ~isfield(features{k}.fparams, 'input_size_mode')
            features{k}.fparams.input_size_mode = 'adaptive';
        end
        if ~isfield(features{k}.fparams, 'input_size_scale')
            features{k}.fparams.input_size_scale = 1;
        end
        if ~isfield(features{k}.fparams, 'downsample_factor')
            features{k}.fparams.downsample_factor = ones(1, length(features{k}.fparams.output_layer));
        end

        % load the network
        net = load_cnn(features{k}.fparams, img_sample_sz);

        % find the dimensionality of each layer
        features{k}.fparams.nDim = net.info.dataSize(3, features{k}.fparams.output_layer+1)';

        % find the stride of the layers
        if isfield(net.info, 'receptiveFieldStride')
            net_info_stride = cat(2, [1; 1], net.info.receptiveFieldStride);
        else
            net_info_stride = [1; 1];
        end

        % compute the cell size of the layers (takes down-sampling factor
        % into account)
        features{k}.fparams.cell_size = net_info_stride(1, features{k}.fparams.output_layer+1)' .* features{k}.fparams.downsample_factor';

        % this feature will always return a cell array
        features{k}.is_cell = true;
        features{k}.is_cnn = true;
    elseif isequal(features{k}.getFeature,@get_eitel_cnn)
        features{k}.fparams = make_eitel_feature(features{k}.fparams);
    else
        error('Unknown feature type');
    end

    % Set default cell size
    if ~isfield(features{k}.fparams, 'cell_size')
        features{k}.fparams.cell_size = 1;
    end

    % Set default penalty
    if ~isfield(features{k}.fparams, 'penalty')
        features{k}.fparams.penalty = zeros(length(features{k}.fparams.nDim),1);
    end

    % Find the minimum cell size of each layer
    feature_info.min_cell_size(k) = min(features{k}.fparams.cell_size);
end

% Order the features in increasing minimal cell size
[~, feat_ind] = sort(feature_info.min_cell_size);
features = features(feat_ind);
feature_info.min_cell_size = feature_info.min_cell_size(feat_ind);

% Set feature info
feature_info.dim_block = cell(num_features,1);
feature_info.penalty_block = cell(num_features,1);

for k = 1:length(features)
    % update feature info
    feature_info.dim_block{k} = features{k}.fparams.nDim;
    feature_info.penalty_block{k} = features{k}.fparams.penalty(:);
end
% Feature info for each cell block
feature_info.dim = cell2mat(feature_info.dim_block);
feature_info.penalty = cell2mat(feature_info.penalty_block);

% Find if there is any CNN feature
cnn_feature_ind = -1;
for k = 1:length(features)
    if features{k}.is_cnn
        cnn_feature_ind = k;
    end
end

% This ugly code sets the image sample size to be used for extracting the
% features. It then computes the data size (size of the features) and the
% image support size (the corresponding size in the image).
if cnn_feature_ind > 0
    scale = features{cnn_feature_ind}.fparams.input_size_scale;

    new_img_sample_sz = img_sample_sz;

    % First try decrease one
    net_info = net.info;

    if ~strcmpi(size_mode, 'same') && strcmpi(features{cnn_feature_ind}.fparams.input_size_mode, 'adaptive')
        orig_sz = net.info.dataSize(1:2,end)' / features{cnn_feature_ind}.fparams.downsample_factor(end);

        if strcmpi(size_mode, 'exact')
            desired_sz = orig_sz + 1;
        elseif strcmpi(size_mode, 'odd_cells')
            desired_sz = orig_sz + 1 + mod(orig_sz,2);
        end

        while desired_sz(1) > net_info.dataSize(1,end)
            new_img_sample_sz = new_img_sample_sz + [1, 0];
            net_info = vl_simplenn_display(net, 'inputSize', [round(scale * new_img_sample_sz), 3 1]);
        end
        while desired_sz(2) > net_info.dataSize(2,end)
            new_img_sample_sz = new_img_sample_sz + [0, 1];
            net_info = vl_simplenn_display(net, 'inputSize', [round(scale * new_img_sample_sz), 3 1]);
        end
    end

    feature_info.img_sample_sz = round(new_img_sample_sz);

    if strcmpi(features{cnn_feature_ind}.fparams.input_size_mode, 'adaptive')
        features{cnn_feature_ind}.img_input_sz = feature_info.img_sample_sz;
    else
        features{cnn_feature_ind}.img_input_sz = net.meta.normalization.imageSize(1:2);
    end

    % Sample size to be input to the net
    scaled_sample_sz = round(scale * features{cnn_feature_ind}.img_input_sz);

    if isfield(net_info, 'receptiveFieldStride')
        net_info_stride = cat(2, [1; 1], net_info.receptiveFieldStride);
    else
        net_info_stride = [1; 1];
    end

    net_stride = net_info_stride(:, features{cnn_feature_ind}.fparams.output_layer+1)';
    total_feat_sz = net_info.dataSize(1:2, features{cnn_feature_ind}.fparams.output_layer+1)';

    shrink_number = max(2 * ceil((net_stride(end,:) .* total_feat_sz(end,:) - scaled_sample_sz) ./ (2 * net_stride(end,:))), 0);

    deepest_layer_sz = total_feat_sz(end,:) - shrink_number;
    scaled_support_sz = net_stride(end,:) .* deepest_layer_sz;

    % Calculate output size for each layer
    cnn_output_sz = round(bsxfun(@rdivide, scaled_support_sz, net_stride));
    features{cnn_feature_ind}.fparams.start_ind = floor((total_feat_sz - cnn_output_sz)/2) + 1;
    features{cnn_feature_ind}.fparams.end_ind = features{cnn_feature_ind}.fparams.start_ind + cnn_output_sz - 1;

    feature_info.img_support_sz = round(scaled_support_sz .* feature_info.img_sample_sz ./ scaled_sample_sz);

    % Set the input size
    features{cnn_feature_ind}.fparams.net = set_cnn_input_size(net, feature_info.img_sample_sz);

    if gparams.use_gpu
        if isempty(gparams.gpu_id)
            gpuDevice();
        elseif gparams.gpu_id > 0
            gpuDevice(gparams.gpu_id);
        end
        features{cnn_feature_ind}.fparams.net = vl_simplenn_move(features{cnn_feature_ind}.fparams.net, 'gpu');
    end
else
    max_cell_size = max(feature_info.min_cell_size);

    if strcmpi(size_mode, 'same')
        feature_info.img_sample_sz = round(img_sample_sz);
    elseif strcmpi(size_mode, 'exact')
        feature_info.img_sample_sz = round(img_sample_sz / max_cell_size) * max_cell_size;
    elseif strcmpi(size_mode, 'odd_cells')
        new_img_sample_sz = (1 + 2*round(img_sample_sz / (2*max_cell_size))) * max_cell_size;

        % Check the size with the largest number of odd dimensions (choices in the
        % third dimension)
        feature_sz_choices = floor(bsxfun(@rdivide, bsxfun(@plus, new_img_sample_sz, reshape(0:max_cell_size-1, 1, 1, [])), feature_info.min_cell_size));
        num_odd_dimensions = sum(sum(mod(feature_sz_choices, 2) == 1, 1), 2);
        [~, best_choice] = max(num_odd_dimensions(:));
        pixels_added = best_choice - 1;
        feature_info.img_sample_sz = round(new_img_sample_sz + pixels_added);
    else
        error('Unknown size_mode');
    end

    % Setting the feature size and support size
    %     feature_info.data_sz = floor(bsxfun(@rdivide, feature_info.img_sample_sz, feature_info.min_cell_size));
    feature_info.img_support_sz = feature_info.img_sample_sz;
end

% Set the sample size and data size for each feature
feature_info.data_sz_block = cell(num_features,1);
for k = 1:length(features)
    if features{k}.is_cnn
        % CNN features have a different sample size, since the receptive
        % field is often larger than the support size
        features{k}.img_sample_sz = feature_info.img_sample_sz(:)';

        % Set the data size based on the computed output size
        feature_info.data_sz_block{k} = floor(bsxfun(@rdivide, cnn_output_sz, features{k}.fparams.downsample_factor'));
    else
        % implemented classic features always have the same sample and
        % support size
        features{k}.img_sample_sz = feature_info.img_support_sz(:)';
        features{k}.img_input_sz = features{k}.img_sample_sz;

        % Set data size based on cell size
        feature_info.data_sz_block{k} = floor(bsxfun(@rdivide, features{k}.img_sample_sz, features{k}.fparams.cell_size));
    end
end

feature_info.data_sz = cell2mat(feature_info.data_sz_block);

author: martin-danelljan

repo: ECO


text:
%% check Clebsch Gordan Tensor

%% the reference

% some arbitrary rotation
g = rotation.byEuler(-72*degree,-88*degree,-134*degree);

% we want to express the product of two wigner D functions
D1 = WignerD(g,'order',1);
D2 = WignerD(g,'order',2);
D1D2_ref = D1(:) * D2(:).'

%% expansion into Wigner functions of lower order

% zero order component
D0 = WignerD(g,'order',0);
CG0 = ClebschGordanTensor(1,2,0);
C0 = EinsteinSum(CG0,[1 3 -1],D0,[-1,-2],CG0,[2 4 -2])

% first order component
D1 = WignerD(g,'order',1);
CG1 = ClebschGordanTensor(1,2,1);
C1 = EinsteinSum(CG1,[1 3 -1],D1,[-1 -2],CG1,[2 4 -2])

% second order component
D2 = WignerD(g,'order',2);
CG2 = ClebschGordanTensor(1,2,2);
C2 = EinsteinSum(CG2,[1 3 -1],D2,[-1 -2],CG2,[2 4 -2])

% third order component
D3 = WignerD(g,'order',3);
CG3 = ClebschGordanTensor(1,2,3);
C3 = EinsteinSum(CG3,[1 3 -1],D3,[-1 -2],CG3,[2 4 -2])

a = reshape(matrix(C0 + C1 + C2 + C3),[9,25]) ./ D1D2_ref;

imagesc(real(a))

author: mtex-toolbox

repo: mtex

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


text:
function [indexSelect, infoChange] = ivmSelectPoint(model, add);

% IVMSELECTPOINT Choose a point for inclusion or removal.
% FORMAT
% DESC identifies the next point for inclusion or removal.
% ARG model : IVM structure for which the next point is being
% selected.
% ARG add : flag which indicates whether or not we are adding a
% point. If we are not adding we are assumed to be removing a point
% (default is true).
%
% SEEALSO : ivmOptimiseIvm, ivmSelectPoints, ivmCreate
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2007

% IVM

if nargin < 2
  % If add is 1, then we are including a point.
  add = true;
end

switch model.selectionCriterion
 case 'random'
  if add
    indexSelect = ceil(rand(1)*length(model.J));
    infoChange = -.5*sum(log2(1-model.varSigma(indexSelect, :).* ...
                              model.nu(indexSelect, :)), 2);
  else
    indexSelect = ceil(rand(1)*length(model.I));
    infoChange = -.5*sum(log2(1-model.varSigma(indexSelect, :).* ...
                          model.beta(indexSelect, :)+1e-300), 2);
  end
 case 'entropy'
  delta = ivmComputeInfoChange(model, add);
  [infoChange, indexSelect] = max(delta);
  numSelect = sum(delta==infoChange);
  if numSelect>1
    index1 = find(delta==infoChange);
    index1Select = ceil(rand(1)*numSelect);
    indexSelect = index1(index1Select);
  end
 case 'rentropy'
  % entropy with first point random
  if length(model.I)
    % if point is already selected select another.
    delta = ivmComputeInfoChange(model, add);
    [infoChange, indexSelect] = max(delta);
  else
    % otherwise select one randomly
    if add
      indexSelect = ceil(rand(1)*length(model.J));
      infoChange = -.5*sum(log2(1-model.varSigma(indexSelect, :)* ...
                                model.nu(indexSelect, :)), 2);
    else
      indexSelect = ceil(rand(1)*length(model.I));
      infoChange = -.5*sum(log2(1-model.varSigma(indexSelect, :)* ...
                                model.beta(indexSelect, :)+1e-300), 2);
    end
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [params, names] = cmpndKernExtractParam(kern)

% CMPNDKERNEXTRACTPARAM Extract parameters from the CMPND kernel structure.
% FORMAT
% DESC Extract parameters from the compound kernel matrix into a vector of
% parameters for optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the
% kernel. The vector of 'transforms' is assumed to be empty
% here. Any transformations of parameters should be done in
% component kernels.
%
% SEEALSO cmpndKernParamInit, cmpndKernExpandParam, kernExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

params = zeros(1, kern.nParams);
if nargout > 1
  namesTemp = cell(1, kern.nParams);
end
startVal = 1;
endVal = 0;
storedTypes = cell(0);
for i = 1:length(kern.comp)
  endVal = endVal + kern.comp{i}.nParams;
  if nargout > 1
    [params(1, startVal:endVal), namesTemp(startVal:endVal)] = kernExtractParam(kern.comp{i});
    instNum = sum(strcmp(kern.comp{i}.type, storedTypes)) + 1;
    for j = startVal:endVal
      namesTemp{1, j} = [kern.comp{i}.type ' ' num2str(instNum) ' ' ...
                         namesTemp{1, j}];
    end
    storedTypes{end+1} = kern.comp{i}.type;
  else
    params(1, startVal:endVal) = kernExtractParam(kern.comp{i});
  end
    startVal = endVal + 1;
end

% If any parameters are 'tied together' deal with them.
paramGroups = kern.paramGroups;
for i = 1:size(paramGroups, 2)
  ind = find(paramGroups(:, i));
  if nargout > 1
    names{i} = namesTemp{ind(1)};
    if length(ind) > 1
      for j = 2:length(ind)
        names{i} = [names{i} ', ' namesTemp{ind(j)}];
      end
    end
  end
  paramGroups(ind(2:end), i) = 0;
end
params = params*paramGroups;


text:
function [recall, precision, info] = vl_pr(labels, scores, varargin)
%VL_PR   Precision-recall curve.
%   [RECALL, PRECISION] = VL_PR(LABELS, SCORES) computes the
%   precision-recall (PR) curve. LABELS are the ground truth labels,
%   greather than zero for a positive sample and smaller than zero for
%   a negative one. SCORES are the scores of the samples obtained from
%   a classifier, where lager scores should correspond to positive
%   samples.
%
%   Samples are ranked by decreasing scores, starting from rank 1.
%   PRECISION(K) and RECALL(K) are the precison and recall when
%   samples of rank smaller or equal to K-1 are predicted to be
%   positive and the remaining to be negative. So for example
%   PRECISION(3) is the percentage of positive samples among the two
%   samples with largest score. PRECISION(1) is the precision when no
%   samples are predicted to be positive and is conventionally set to
%   the value 1.
%
%   Set the zero the lables of samples that should be ignored in the
%   evaluation. Set to -INF the scores of samples which are not
%   retrieved. If there are samples with -INF score, then the PR curve
%   may have maximum recall smaller than 1, unless the INCLUDEINF
%   option is used (see below). The options NUMNEGATIVES and
%   NUMPOSITIVES can be used to specify additional samples with -INF
%   score (see below).
%
%   [RECALL, PRECISION, INFO] = VL_PR(...) returns an additional
%   structure INFO with the following fields:
%
%   info.auc::
%     The area under the precision-recall curve. If the INTERPOLATE
%     option is set to FALSE, then trapezoidal interpolation is used
%     to integrate the PR curve. If the INTERPOLATE option is set to
%     TRUE, then the curve is piecewise constant and no other
%     approximation is introduced in the calculation of the area. In
%     the latter case, INFO.AUC is the same as INFO.AP.
%
%   info.ap::
%     Average precision as defined by TREC. This is the average of the
%     precision observed each time a new positive sample is
%     recalled. In this calculation, any sample with -INF score
%     (unless INCLUDEINF is used) and any additional positive induced
%     by NUMPOSITIVES has precision equal to zero. If the INTERPOLATE
%     option is set to true, the AP is computed from the interpolated
%     precision and the result is the same as INFO.AUC. Note that AP
%     as defined by TREC normally does not use interpolation [1].
%
%   info.ap_interp_11::
%     11-points interpolated average precision as defined by TREC.
%     This is the average of the maximum precision for recall levels
%     greather than 0.0, 0.1, 0.2, ..., 1.0. This measure was used in
%     the PASCAL VOC challenge up to the 2008 edition.
%
%   info.auc_pa08::
%     Deprecated. It is the same of INFO.AP_INTERP_11.
%
%   VL_PR(...) with no output arguments plots the PR curve in the
%   current axis.
%
%   VL_PR() accepts the following options:
%
%   Interpolate:: false
%     If set to true, use interpolated precision. The interpolated
%     precision is defined as the maximum precision for a given recall
%     level and onwards. Here it is implemented as the culumative
%     maximum from low to high scores of the precision.
%
%   NumPositives:: []
%   NumNegatives:: []
%     If set to a number, pretend that LABELS contains this may
%     positive/negative labels. NUMPOSITIVES/NUMNEGATIVES cannot be
%     smaller than the actual number of positive/negative entrires in
%     LABELS. The additional positive/negative labels are appended to
%     the end of the sequence, as if they had -INF scores (not
%     retrieved). This is useful to evaluate large retrieval systems
%     for which one stores ony a handful of top results for efficiency
%     reasons.
%
%   IncludeInf:: false
%     If set to true, data with -INF score SCORES is included in the
%     evaluation and the maximum recall is 1 even if -INF scores are
%     present. This option does not include any additional positive or
%     negative data introduced by specifying NUMPOSITIVES and
%     NUMNEGATIVES.
%
%   Stable:: false
%     If set to true, RECALL and PRECISION are returned the same order
%     of LABELS and SCORES rather than being sorted by decreasing
%     score (increasing recall). Samples with -INF scores are assigned
%     RECALL and PRECISION equal to NaN.
%
%   About the PR curve::
%     This section uses the same symbols used in the documentation of
%     the VL_ROC() function. In addition to those quantities, define:
%
%       PRECISION(S) = TP(S) / (TP(S) + FP(S))
%       RECALL(S) = TPR(S) = TP(S) / P
%
%     The precision is the fraction of positivie predictions which are
%     correct, and the recall is the fraction of positive labels that
%     have been correctly classified (recalled). Notice that the recall
%     is also equal to the true positive rate for the ROC curve (see
%     VL_ROC()).
%
%   REFERENCES:
%   [1] C. D. Manning, P. Raghavan, and H. Schutze. An Introduction to
%   Information Retrieval. Cambridge University Press, 2008.
%
%   See also VL_ROC(), VL_HELP().

% Author: Andrea Vedaldi

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% 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).

[tp, fp, p, n, perm, varargin] = vl_tpfp(labels, scores, varargin{:}) ;
opts.stable = false ;
opts.interpolate = false ;
opts = vl_argparse(opts,varargin) ;

% compute precision and recall
small = 1e-10 ;
recall = tp / max(p, small) ;
precision = max(tp, small) ./ max(tp + fp, small) ;

% interpolate precision if needed
if opts.interpolate
  precision = fliplr(vl_cummax(fliplr(precision))) ;
end

% --------------------------------------------------------------------
%                                                      Additional info
% --------------------------------------------------------------------

if nargout > 2 || nargout == 0

  % area under the curve using trapezoid interpolation
  if ~opts.interpolate
    info.auc = 0.5 * sum((precision(1:end-1) + precision(2:end)) .* diff(recall)) ;
  end

  % average precision (for each recalled positive sample)
  sel = find(diff(recall)) + 1 ;
  info.ap = sum(precision(sel)) / p ;
  if opts.interpolate
    info.auc = info.ap ;
  end

  % TREC 11 points average interpolated precision
  info.ap_interp_11 = 0.0 ;
  for rc = linspace(0,1,11)
    pr = max([0, precision(recall >= rc)]) ;
    info.ap_interp_11 = info.ap_interp_11 + pr / 11 ;
  end

  % legacy definition
  info.auc_pa08 = info.ap_interp_11 ;
end

% --------------------------------------------------------------------
%                                                                 Plot
% --------------------------------------------------------------------

if nargout == 0
  cla ; hold on ;
  plot(recall,precision,'linewidth',2) ;
  spline([0 1], [1 1] * p / length(labels), 'r--', 'linewidth', 2) ;
  axis square ; grid on ;
  xlim([0 1]) ; xlabel('recall') ;
  ylim([0 1]) ; ylabel('precision') ;
  title(sprintf('PR (AUC: %.2f%%, AP: %.2f%%, AP11: %.2f%%)', ...
                info.auc * 100, ...
                info.ap * 100, ...
                info.ap_interp_11 * 100)) ;
  if opts.interpolate
    legend('PR interp.', 'PR rand.', 'Location', 'SouthEast') ;
  else
    legend('PR', 'PR rand.', 'Location', 'SouthEast') ;
  end
  clear recall precision info ;
end

% --------------------------------------------------------------------
%                                                        Stable output
% --------------------------------------------------------------------

if opts.stable
  precision(1) = [] ;
  recall(1) = [] ;
  precision_ = precision ;
  recall_ = recall ;
  precision = NaN(size(precision)) ;
  recall = NaN(size(recall)) ;
  precision(perm) = precision_ ;
  recall(perm) = recall_ ;
end

% --------------------------------------------------------------------
function h = spline(x,y,spec,varargin)
% --------------------------------------------------------------------
prop = vl_linespec2prop(spec) ;
h = line(x,y,prop{:},varargin{:}) ;

author: Cloud-CV

repo: object-proposals


text:
function fd=pderiv(f,dim,difforder)
%PDERIV   Derivative of smooth periodic function
%   Usage:  fd=pderiv(f);
%           fd=pderiv(f,dim);
%           fd=pderiv(f,dim,difforder);
%
%   `pderiv(f)` will compute the derivative of *f* using a using a 4th order
%   centered finite difference scheme. *f* must have been obtained by a
%   regular sampling. If *f* is a matrix, the derivative along the columns
%   will be found.
%
%   `pderiv(f,dim)` will do the same along dimension *dim*.
%
%   `pderiv(f,dim,difforder)` uses a centered finite difference scheme of
%   order *difforder* instead of the default.
%
%   `pderiv(f,dim,Inf)` will compute the spectral derivative using a DFT.
%
%   `pderiv` assumes that *f* is a regular sampling of a function on the
%   torus $[0,1)$. The derivative of a function on a general torus $[0,T)$
%   can be found by scaling the output by $1/T$.

% Assert correct input.

complainif_argnonotinrange(nargin,1,3,mfilename);

if nargin==1
  dim=[];
end;

if nargin<3
  difforder=4;
end;

[f,L,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,[],dim,'PDERIV');

switch(difforder)
 case 2
  fd = L*(circshift(f,-1)-circshift(f,1))/2;
 case 4
  fd = L*(-circshift(f,-2)+8*circshift(f,-1)-8*circshift(f,1)+ ...
          circshift(f,2))/12;
 case Inf
  n=fftindex(L,0);
  n=repmat(n,1,W);

  fd=2*pi*ifft(i*n.*fft(f));

  if isreal(f)
    fd=real(fd);
  end;

 otherwise
  error('The specified differentation order is not implemented.');
end;

fd=assert_sigreshape_post(fd,dim,permutedsize,order);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function v = asgn(a1, a2) %#ok
% Author Bruno Luong <brunoluong@yahoo.com>
% Last update: 05/Apil/2009
v = a2;
end

author: Sable

repo: mcbench-benchmarks


text:
function [V,S] = alphavol(X,R,fig)
% ALPHAVOL Alpha shape of 2D or 3D point set.
%
%   V = ALPHAVOL(X,R) gives the area or volume V of the basic alpha shape
%   for a 2D or 3D point set. X is a coordinate matrix of size Nx2 or Nx3.
%
%   R is the probe radius with default value R = Inf. In the default case
%   the basic alpha shape (or alpha hull) is the convex hull.
%
%   [V,S] = ALPHAVOL(X,R) outputs a structure S with fields:
%    S.tri - Triangulation of the alpha shape (Mx3 or Mx4)
%    S.vol - Area or volume of simplices in triangulation (Mx1)
%    S.rcc - Circumradius of simplices in triangulation (Mx1)
%    S.bnd - Boundary facets (Px2 or Px3)
%
%   ALPHAVOL(X,R,1) plots the alpha shape.
%
%   % 2D Example - C shape
%   t = linspace(0.6,5.7,500)';
%   X = 2*[cos(t),sin(t)] + rand(500,2);
%   subplot(221), alphavol(X,inf,1);
%   subplot(222), alphavol(X,1,1);
%   subplot(223), alphavol(X,0.5,1);
%   subplot(224), alphavol(X,0.2,1);
%
%   % 3D Example - Sphere
%   [x,y,z] = sphere;
%   [V,S] = alphavol([x(:),y(:),z(:)]);
%   trisurf(S.bnd,x,y,z,'FaceColor','blue','FaceAlpha',1)
%   axis equal
%
%   % 3D Example - Ring
%   [x,y,z] = sphere;
%   ii = abs(z) < 0.4;
%   X = [x(ii),y(ii),z(ii)];
%   X = [X; 0.8*X];
%   subplot(211), alphavol(X,inf,1);
%   subplot(212), alphavol(X,0.5,1);
%
%   See also DELAUNAY, TRIREP, TRISURF

%   Author: Jonas Lundgren <splinefit@gmail.com> 2010

%   2010-09-27  First version of ALPHAVOL.
%   2010-10-05  DelaunayTri replaced by DELAUNAYN. 3D plots added.
%   2012-03-08  More output added. DELAUNAYN replaced by DELAUNAY.

if nargin < 2 || isempty(R), R = inf; end
if nargin < 3, fig = 0; end

% Check coordinates
dim = size(X,2);
if dim < 2 || dim > 3
    error('alphavol:dimension','X must have 2 or 3 columns.')
end

% Check probe radius
if ~isscalar(R) || ~isreal(R) || isnan(R)
    error('alphavol:radius','R must be a real number.')
end

% Unique points
[X,imap] = unique(X,'rows');

% Delaunay triangulation
T = delaunay(X);

% Remove zero volume tetrahedra since
% these can be of arbitrary large circumradius
if dim == 3
    n = size(T,1);
    vol = volumes(T,X);
    epsvol = 1e-12*sum(vol)/n;
    T = T(vol > epsvol,:);
    holes = size(T,1) < n;
end

% Limit circumradius of simplices
[~,rcc] = circumcenters(TriRep(T,X));
T = T(rcc < R,:);
rcc = rcc(rcc < R);

% Volume/Area of alpha shape
vol = volumes(T,X);
V = sum(vol);

% Return?
if nargout < 2 && ~fig
    return
end

% Turn off TriRep warning
warning('off','MATLAB:TriRep:PtsNotInTriWarnId')

% Alpha shape boundary
if ~isempty(T)
    % Facets referenced by only one simplex
    B = freeBoundary(TriRep(T,X));
    if dim == 3 && holes
        % The removal of zero volume tetrahedra causes false boundary
        % faces in the interior of the volume. Take care of these.
        B = trueboundary(B,X);
    end
else
    B = zeros(0,dim);
end

% Plot alpha shape
if fig
    if dim == 2
        % Plot boundary edges and point set
        x = X(:,1);
        y = X(:,2);
        plot(x(B)',y(B)','r','linewidth',2), hold on
        plot(x,y,'k.'), hold off
        str = 'Area';
    elseif ~isempty(B)
        % Plot boundary faces
        trisurf(TriRep(B,X),'FaceColor','red','FaceAlpha',1/3);
        str = 'Volume';
    else
        cla
        str = 'Volume';
    end
    axis equal
    str = sprintf('Radius = %g,   %s = %g',R,str,V);
    title(str,'fontsize',12)
end

% Turn on TriRep warning
warning('on','MATLAB:TriRep:PtsNotInTriWarnId')

% Return structure
if nargout == 2
    S = struct('tri',imap(T),'vol',vol,'rcc',rcc,'bnd',imap(B));
end

%--------------------------------------------------------------------------
function vol = volumes(T,X)
%VOLUMES Volumes/areas of tetrahedra/triangles

% Empty case
if isempty(T)
    vol = zeros(0,1);
    return
end

% Local coordinates
A = X(T(:,1),:);
B = X(T(:,2),:) - A;
C = X(T(:,3),:) - A;

if size(X,2) == 3
    % 3D Volume
    D = X(T(:,4),:) - A;
    BxC = cross(B,C,2);
    vol = dot(BxC,D,2);
    vol = abs(vol)/6;
else
    % 2D Area
    vol = B(:,1).*C(:,2) - B(:,2).*C(:,1);
    vol = abs(vol)/2;
end

%--------------------------------------------------------------------------
function B = trueboundary(B,X)
%TRUEBOUNDARY True boundary faces
%   Remove false boundary caused by the removal of zero volume
%   tetrahedra. The input B is the output of TriRep/freeBoundary.

% Surface triangulation
facerep = TriRep(B,X);

% Find edges attached to two coplanar faces
E0 = edges(facerep);
E1 = featureEdges(facerep, 1e-6);
E2 = setdiff(E0,E1,'rows');

% Nothing found
if isempty(E2)
    return
end

% Get face pairs attached to these edges
% The edges connects faces into planar patches
facelist = edgeAttachments(facerep,E2);
pairs = cell2mat(facelist);

% Compute planar patches (connected regions of faces)
n = size(B,1);
C = sparse(pairs(:,1),pairs(:,2),1,n,n);
C = C + C' + speye(n);
[~,p,r] = dmperm(C);

% Count planar patches
iface = diff(r);
num = numel(iface);

% List faces and vertices in patches
facelist = cell(num,1);
vertlist = cell(num,1);
for k = 1:num

    % Neglect singel face patches, they are true boundary
    if iface(k) > 1

        % List of faces in patch k
        facelist{k} = p(r(k):r(k+1)-1);

        % List of unique vertices in patch k
        vk = B(facelist{k},:);
        vk = sort(vk(:))';
        ik = [true,diff(vk)>0];
        vertlist{k} = vk(ik);

    end
end

% Sort patches by number of vertices
ivert = cellfun(@numel,vertlist);
[ivert,isort] = sort(ivert);
facelist = facelist(isort);
vertlist = vertlist(isort);

% Group patches by number of vertices
p = [0;find(diff(ivert));num] + 1;
ipatch = diff(p);

% Initiate true boundary list
ibound = true(n,1);

% Loop over groups
for k = 1:numel(ipatch)

    % Treat groups with at least two patches and four vertices
    if ipatch(k) > 1 && ivert(p(k)) > 3

        % Find double patches (identical vertex rows)
        V = cell2mat(vertlist(p(k):p(k+1)-1));
        [V,isort] = sortrows(V);
        id = ~any(diff(V),2);
        id = [id;0] | [0;id];
        id(isort) = id;

        % Deactivate faces in boundary list
        for j = find(id')
            ibound(facelist{j-1+p(k)}) = 0;
        end

    end
end

% Remove false faces
B = B(ibound,:);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function VoigtSpectrGUI
%this is a GUI program i wrote about simulating absorption lineshape and
%spectrum from HITRAN output data
% this program simulate the spectrum using the HITRAN data and
% Calculating the Voigt porfile According to the paper:
%   Applied spectro. V58, 468(2004);
%   Phi= A*K(x,y)
% The lineshape depends on Guassian (vD) width and Lorentzian width(vL).
% Absorption depends on number density and absorption length, I/I0 = Exp(-aNL);
% Notice the Half Width at Half Maxium for the vD and vL inputs;

    fh = figure('Visible','off','Name','Spectrum Simulation',...
             'position',[50,50, 850, 600]);

    hax = axes('Units','pixels','Position',[60,60,500,500]);

    %construct the components
    hp = uipanel('Units','pixels','ShadowColor',[0.1,1,0.01],...
              'Position',[628 160 200 430]);

    hpb1 = uicontrol('Style','pushbutton','String','Open HITRFile',...
          'Foregroundcolor',[1 0.1 0.2],'FontWeight','bold',...
           'Position',[650 540, 110, 30],'Callback',{@fileopen_Callback});

    het1 = uicontrol('Style','edit','String','file name',...
           'Foregroundcolor',[0.2 0.4 0.1],'FontWeight','bold',...
           'Position',[650,510,110,30],'Callback',{@edittext_Callback});

    hst1 = uicontrol('Style','text', 'String','HWHM[vD vL]/cm-1',...
            'Foregroundcolor',[0.2 0.3 0.3],'FontWeight','bold',...
             'Backgroundcolor',[0.7 0.98 0.98],...
             'Position', [650 470 160 15]);

    het2 = uicontrol('Style','edit', 'String','0.003',...
            'Foregroundcolor',[0.30 0.1 1],'FontWeight','bold',...
             'Backgroundcolor',[0.97842 0.98 1],...
             'Position', [650 445 75 20]);

    het3 = uicontrol('Style','edit', 'String','0.006',...
            'Foregroundcolor',[0.30 0.1 1],'FontWeight','bold',...
             'Backgroundcolor',[1 0.99 1],...
             'Position', [735 445 75 20]);

    hst2 = uicontrol('Style','text', 'String','NumDensity/Lightpath/cm',...
            'Foregroundcolor',[0.2 0.3 0.3],'FontWeight','bold',...
             'Backgroundcolor',[0.7 0.98 1],...
             'Position', [650 410 160 15]);

    het4 = uicontrol('Style','edit', 'String','1E+10',...
            'Foregroundcolor',[0.30 0.1 1],'FontWeight','bold',...
             'Backgroundcolor',[0.97842 0.98 1],...
             'Position', [650 385 75 20]);

    het5 = uicontrol('Style','edit', 'String','1E+5',...
            'Foregroundcolor',[0.30 0.1 1],'FontWeight','bold',...
             'Backgroundcolor',[1 0.99 1],...
             'Position', [735 385 75 20]);

    hpb2 = uicontrol('Style','pushbutton','String','Simu-LineShape',...
        'Foregroundcolor',[1 0.1 0.2],'FontWeight','bold',...
           'Position',[650 340, 110, 30],'Callback',{@SimuLS_Callback});

    het6 = uicontrol('Style','edit', 'String','VoigtWidth',...
        'Foregroundcolor',[0.30 0.1 1],'FontWeight','bold',...
         'Backgroundcolor',[0.4 0.9 0.7],...
         'Position', [735 315 75 20]);

    hpb3 = uicontrol('Style','pushbutton','String','Simu-Absorption',...
        'Foregroundcolor',[1 0.1 0.2],'FontWeight','bold',...
           'Position',[650 280, 110, 30],'Callback',{@SimuAp_Callback});

    hpb4 = uicontrol('Style','pushbutton','String','Open ExpSpec',...
        'Foregroundcolor',[1 0.1 0.2],'FontWeight','bold',...
           'Position',[650 240, 110, 30],'Callback',{@ExpSpectr_Callback});

    hst4 = uicontrol('Style','text', 'String','Display Options',...
          'Backgroundcolor',[0.678402 0.98 1],...
            'Foregroundcolor',[0.2 0.3 0.3],'FontWeight','bold',...
             'Position',[650 210 130 15]);

    hpm = uicontrol('Style','popupmenu',...
                'String',{'LineShape','Absorption','CombinedPlot'},...
                'Foregroundcolor',[0.99 0.01 0.01],'FontWeight','bold',...
                'Value',1,'Position',[690 170 110 30],'Callback',{@popup_Callback});

    hp2 = uipanel('Units','pixels','ShadowColor',[0.1 0.1 1],...
               'Position',[628 40 200 100]);

    het7 = uicontrol('Style','edit','String','FileName.dat',...
             'Foregroundcolor',[0.2 0.4 0.1],'FontWeight','bold',...
             'Position',[640,105,90,25]);

    hcb1 = uicontrol('Style','checkbox','String','Line Shape',...
                'Position',[655 75 85 25],'Foregroundcolor',[0.2 0.3 0.9]);
    hcb2 = uicontrol('Style','checkbox','String','Absorption',...
                'Position',[655 50 85 25],'Foregroundcolor',[0.2 0.3 0.9]);

    hpb5 = uicontrol('Style','pushbutton','String','Save',...
              'Foregroundcolor',[0.6 0.2 0.3],'FontWeight','bold',...
              'Position',[750 66, 60, 50],'Callback',{@save_Callback});

      % to align the controls
    align([hp,het1,het6,hpb1,hpb2,hpb3,hpb4,hst1,hst2,hst4,hpm,hp2],'Center','None');

    % Initialize the GUI.
    % Change units to normalized so components resize automatically.
    set([hp,het1,het2,het3,het4,het5,het6,het7,hpb1,hpb2,hax,hst1,hst2,...
       hst4,hpm,hpb3,hpb4,hp2,hcb1,hcb2,hpb5],'Units','normalized');

    %make figure visible.
    set(fh, 'Visible','on');
    set(fh, 'NumberTitle','off')%% supress the NumverTitle

    % define some variables
    format long e ; %% define numerical format
    stp=0.0002; %% calculation(convolution) stepsize
    Dd = 0.5; %% define the edge extended beyond the spectrum region
    Nd = 1E+10; %% Number density
    L  = 1E+5; %% Optical path length cm

    gL = 0.001; %% in cm -1, pressure broadening cofffient at 50 torr.
    gD = 0.003; %% 0.0032 doppler broadening in cm-1
    SgmvTot= zeros(1,10);  %% initiate the arrays
    expx = zeros(1,10);
    expy = zeros(1,10);
    A = zeros(1,10);
    dtfl = zeros(1,10);
    grd = zeros(1,10); %% default grid

    %%file open callback
    function fileopen_Callback(source,eventdata)
        % **********uiputfile for open the save file dialog*********

        [flnm, flpth] = uigetfile({'*.out','All Files' },...
                         'Select HITRAN outputfile','fl.out');
        %specify application data, so it can be used by other objects;
        % setappdata(h,'name',value)
        setappdata(hpb1,'fname',flnm);

        set(het1, 'String',flnm);
        % set the edit text box value.

     end

  % read in the HITRAN data and simulate Line Shape spectrum
    function SimuLS_Callback(source, eventdata)%callback is just a regular func

        format long e ;
        % get the value of specified data
        % value = getappdata(h,name)
        % values = getappdata(h)
        flnm = getappdata(hpb1,'fname');

        fid = fopen(flnm);
        C = textscan(fid,...
        '%d%f%f%u%f%f%f%f %d%d%d%d%d %d%d%d%d %d%d%d%d %d%d%d%d%d%d%d%d%d%d%d%d%d ');
        fclose(fid);  %% returned C is a cell

        N = length(C{2}); %% Number of lines/ Cell length
        %% define a grid
        rangeL= -Dd + min(C{2});  %% spectrum region: First line minus 1 cm-1
        rangeH = Dd + max(C{2}); %% Last line position + 1 cm-1
        grd=[rangeL:stp:rangeH];  %% define the grid

        v = grd;
        gD = str2num(get(het2, 'string'));
        gL = str2num(get(het3, 'string'));
        gV = 0.5346*gL + sqrt(0.2166*gL^2 + gD^2); %% Voigt profile half width
        x = gL/gV;
        SgmvTot = 0;

        for i=1:N  %% calculate the line shape for each peak

            v0(i) = C{2}(i);
            S(i) = C{3}(i);
            y = abs(v-v0(i))/gV;
            Sgmv0(i) = S(i)/(2*gV*(1.065 + 0.447*x + 0.058*x^2));
            Sgmv = Sgmv0(i)*((1-x)*exp(-0.693.*y.^2) + (x./(1+y.^2)) + ...
              0.016*(1-x)*x*(exp(-0.0841.*y.^2.25)-1./(1 + 0.021.*y.^2.25)));

         SgmvTot = SgmvTot + Sgmv;

        end

        plot(grd,SgmvTot,'-b.');
        xlabel('WaveLength in cm-1');
        ylabel('Effective absorption Cross Section in cm2');

        set(hpm, 'Value',1)% change the popup menu settings;

        %% display the Voigt width;
        strg = sprintf('gV= %1.4f',gV);
        set(het6,'string',strg);
        set(het6,'Backgroundcolor',[0.9 0.9 0.9])

        %%figure, plot(v,A);
        %%ylabel('Absorption')
    end

    function SimuAp_Callback(source, eventdata)
        %% calculate absorption

        Nd = str2num(get(het4, 'string')); %% Number density in molecule/cm3
        L = str2num(get(het5, 'string')); %% effective length path in cm
        A = 1 - exp(-SgmvTot.*Nd*L); %% I/I0 = Exp(-aNL);

        plot(grd,A,'-m');
        xlabel('WaveLength in cm-1');
        set(hpm, 'Value',2)% change the popup menu settings;

    end

    %slice X callback
    function ExpSpectr_Callback(source, eventdata)

        format long e ;
        [flnm, flpth] = uigetfile({'*.dat','All Files' },'Select Exp. file','ExpSpectr.dat');
        fid = fopen(flnm, 'r');
        a = fscanf(fid, '%f %f', [2 inf]);    % It has two rows now.
        a = a';
        fclose(fid)
        expx = a(:,1);
        expy = a(:,2);
        hold on;
        plot3 = plot(expx,expy,'-.b');
        ylabel('Absorption')

        hold off;

    end

   %popup menu callback
    function popup_Callback(source,eventdata)
        str = get(source, 'String');
        val = get(hpm, 'Value'); %use either function handle 'hpm' or 'source';
        % the string of popup menu is a cell array {}, cell{num} returns num th
        % element of the cell.
        switch val
        case 1  %'showImage'
             %disp(str{val});
                hold off;
                plot(grd,SgmvTot,'-b.');
                xlabel('WaveLength in cm-1')
                ylabel('Effective absorption Cross Section in cm2')

        case 2
                hold off;
                plot(grd,A,'-m');

        case 3 %combined plot
                hold on;
                plot(expx,expy,'-.b');
                ylabel('Absorption')

                hold off;

        end
    end

    function save_Callback(source, eventdata)

        if get(hcb1,'Value')> 0.9 % it should be 1.0
            dtfl = [grd;SgmvTot]; %% create a data matrix for data saving
            fid = fopen(get(het7, 'string'), 'wt');
            fprintf(fid, '%6.4f %e\n', dtfl); %% set the data format
            fclose(fid);

        end

        if get(hcb2,'Value')> 0.9 % it should be 1.0

           dtfl = [grd;A]; %% create a data matrix for data saving
           fid = fopen(get(het7, 'string'), 'wt');
           fprintf(fid, '%6.4f %5.6f\n', dtfl); %% set the data format
           fclose(fid);
        end

        if get(hcb1,'Value')> 0.9 && get(hcb2,'Value')> 0.9 % it should be 1.0

            dtfl = [grd;SgmvTot;A]; %% create a data matrix for data saving

            fid = fopen(get(het7, 'string'), 'wt');
            fprintf(fid, '%6.4f %e %5.6f\n', dtfl); %% set the data format
            fclose(fid);

        end
    end

%%% read in the HITRAN file to obtain line positions and intensities.
end


text:
% SCRIPT TEST THE DIRECT DYNAMICS OF A 2 DOF PLANAR ROBOT ROBOT

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.

close all

fprintf('\nTHE SIMULATION PRESENTS THE ROBOT AT AN INITIAL POSITION WHEN NO TORQUES ARE APPLIED\n')

%load robot parameters
robot=load_robot('example', '2dofplanar');

%simulate for 10 seconds, change this depending on your computer speed and
%the total time that you want to simulate
total_simulation_time = 10;

%initial position and joint speed
q0 = [0 0]';
qd0 = [0 0]';

drawrobot3d(robot, q0);
adjust_view(robot);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The student should try different combinations of the following parameters:
%   g: the direction of the gravity vector. In this case, if we select g=[0  0 9.81]';
%       the movement of the arm is not affected by the gravity, since it
%       moves in a plane perpendicular to the gravity vector.
%   tau: the torques applied to each joint. The student should observe the effects of selecting
%       tau = [0 0 0]' or different values in combination with the
%       direction of the vector g.
%   robot.dynamics.friction = 0 selects no friction at the joints, whereas
%       robot.dynamics.friction = 1 considers that there exists friction. This
%       friction is modelled by robot.motors.Viscous (viscous friction) and
%       robot.motors.Coulomb (Coulomb friction). The student should observe
%       that selecting g=[0  9.81 0]' and tau = [0 0 0]' and
%       robot.dynamics.friction = 0 turns into an infinite triple pendulum
%       movement. In addition, selecting selecting g=[0  9.81 0]' and tau = [0 0 0]' and
%       robot.dynamics.friction = 1 simulates the case in which the triple
%       pendulum converges to a steady solution with the three links
%       hanging along the Y direction.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%you may redefine the gravity vector
%in this case you may one of the next two lines, that define
%the gravity acting along the Y axis or the Z axis, respectively.
g=[0  -9.81 0]'; %y0 axis

tau = [20 20]';%no torques applied
%select friction or not
robot.dynamics.friction = 1;

fprintf('\nCOMPUTING FORWARD DYNAMICS (this may take a while)')

%this may take a while, since it requires integration
%of the acceleration at each time step
[t, q, qd] = forwarddynamic(robot, total_simulation_time, q0, qd0, tau, g, []);

figure, plot(t, q), grid, title('Position vs. time')
xlabel('time (s)'), ylabel('Position (rad)')
legend('q_1', 'q_2', 'q_3', 'q_4', 'q_5', 'q_6');

figure, plot(t, qd), grid, title('Speed vs. time')
xlabel('time (s)'), ylabel('Speed (rad/s)')
legend('qd_1', 'qd_2', 'qd_3', 'qd_4', 'qd_5', 'qd_6');

%plot results faster in animate!
qq = q(:,1:2:end);

%animate it!!
animate(robot, qq)

repo: ARTE


text:
function camAlignImageToDTI(dti_ref_filename,high_res_filename,output_filename,pixdim_tosave)
%
%   camAlignImageToDTI(dti_ref_filename,high_res_filename,output_filename,pixdim_tosave)
%
% EXAMPLE:
%           camAlignImageToDTI('b0.nii.gz','mask_t1.nii.gz','xformed_mask.nii.gz');
%
% (c) Stanford VISTA, 2010

% Experiment with xform
ni_src = niftiRead(high_res_filename);
ni_ref = niftiRead(dti_ref_filename);

ref_dim = ni_ref.dim(1:3);
ref_pixdim = ni_ref.pixdim(1:3);

% If no desired pixel size given, choose the source pixel size
if notDefined('pixdim_tosave')
    pixdim_tosave = ni_src.pixdim(1:3);
end

if pixdim_tosave ~= ni_src.pixdim(1:3)
    % Get bounding box to match the world space of the dti reference
    bound_box = mrAnatXformCoords(ni_ref.qto_xyz,[1 1 1;ref_dim]);
    % Must take into account the pixel size difference that slightly shifts the
    % origin and therefore we must shift the bounding box
    % also (see mrAnatResliceSpm)
    bound_box = bound_box + repmat(rem((pixdim_tosave - ref_pixdim)/2,ref_pixdim),2,1);

    % reslice
    [new_data xform] = mrAnatResliceSpm(ni_src.data, ni_src.qto_ijk, bound_box, pixdim_tosave, [0 0 0 0 0 0]);
else
    new_data = ni_src.data;
    xform = ni_src.qto_xyz;
end

% write
dtiWriteNiftiWrapper(new_data, xform, output_filename);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = integral2(varargin)
%INTEGRAL2  Double integral of a SPHEREFUN over its domain.
%   I = INTEGRAL2(F) returns a value representing the double integral of a
%   SPHEREFUN.
%
% See also SPHEREFUN/INTEGRAL, SPHEREFUN/SUM2, SPHEREFUN/QUAD2D.

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

[varargout{1:nargout}] = sum2(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function varargout = log(varargin)

switch class(varargin{1})

    case 'sdpvar'
         % Try to detect logsumexp construction etc
        varargout{1} = check_for_special_cases(varargin{:});
        % Nope, then just define this logarithm
        if isempty(varargout{1})
            varargout{1} = InstantiateElementWise(mfilename,varargin{:});
        end

    case 'char'

        operator = CreateBasicOperator('concave','increasing','callback');
        operator.derivative = @(x)(1./(abs(x)+eps));
        operator.inverse = @(x)(exp(x));
        operator.domain = [0 inf];
        operator.singularity = [0 -inf -inf] ;

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

function f = check_for_special_cases(x)
f = [];
% Check for log(1+x)
base = getbase(x);
if all(base(:,1)==1)
    f = slog(x-1);
    return;
end
% Check if user is constructing log(sum(exp(x)))
if base(1)~=0
    return
end
if ~all(base(2:end)==1)
    return
end
modelst = yalmip('extstruct',getvariables(x));
if isempty(modelst)
    return;
end
if length(modelst)==1
    models{1} = modelst;
else
    models = modelst;
end
% LOG(DET(X))
if length(models)==1
    if strcmp(models{1}.fcn,'det_internal')
        n = length(models{1}.arg{1});
        try
            f = logdet(reshape(models{1}.arg{1},sqrt(n),sqrt(n)));
        catch
        end
        return
    end
end
% LOG(EXP(x1)+...+EXP(xn))
for i = 1:length(models)
    if ~strcmp(models{i}.fcn,'exp')
        return
    end
end
p = [];
for i = 1:length(models)
    p = [p;models{i}.arg{1}];
end
f = logsumexp(p);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function r = arange(start, stop, step)
    r = [start:step:stop]';

author: petercorke

repo: robotics-toolbox-matlab


text:
function [ r, z, rho, c, s, info ] = dchdd ( r, ldr, p, x, z, ldz, nz, y, rho )

%*****************************************************************************80
%
%% DCHDD downdates an augmented Cholesky decomposition.
%
%  Discussion:
%
%    DCHDD can also downdate the triangular factor of an augmented QR
%    decomposition.
%
%    Specifically, given an upper triangular matrix R of order P, a
%    row vector X, a column vector Z, and a scalar Y, DCHDD
%    determines an orthogonal matrix U and a scalar ZETA such that
%
%          (R   Z )     (RR  ZZ)
%      U * (      )  =  (      ),
%          (0 ZETA)     ( X   Y)
%
%    where RR is upper triangular.
%
%    If R and Z have been obtained from the factorization of a least squares
%    problem, then RR and ZZ are the factors corresponding to the problem
%    with the observation (X,Y) removed.  In this case, if RHO
%    is the norm of the residual vector, then the norm of
%    the residual vector of the downdated problem is
%    sqrt ( RHO * RHO - ZETA * ZETA ). DCHDD will simultaneously downdate
%    several triplets (Z, Y, RHO) along with R.
%
%    For a less terse description of what DCHDD does and how
%    it may be applied, see the LINPACK guide.
%
%    The matrix U is determined as the product U(1)*...*U(P)
%    where U(I) is a rotation in the (P+1,I)-plane of the form
%
%      ( C(I)      -S(I)    )
%      (                    ).
%      ( S(I)       C(I)    )
%
%    The rotations are chosen so that C(I) is real.
%
%    The user is warned that a given downdating problem may be impossible
%    to accomplish or may produce inaccurate results.  For example, this
%    can happen if X is near a vector whose removal will reduce the
%    rank of R.  Beware.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 May 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real R(LDR,P), the upper triangular matrix that
%    is to be downdated.  The part of R below the diagonal is not referenced.
%
%    Input, integer LDR, the leading dimension of the array R.
%    LDR must be at least P.
%
%    Input, integer P, the order of the matrix R.
%
%    Input, real X(P), the row vector that is to be removed from R.
%
%    Input, real Z(LDZ,NZ), an array of NZ P-vectors
%    which are to be downdated along with R.
%
%    Input, integer LDZ, the leading dimension of the array Z.
%    LDZ must be at least P.
%
%    Input, integer NZ, the number of vectors to be downdated.
%    NZ may be zero, in which case Z, Y, and RHO are not referenced.
%
%    Input, real Y(NZ), the scalars for the downdating of
%    the vectors Z.
%
%    Input, real RHO(NZ), the norms of the residual vectors.
%
%    Output, real R(LDR,P), the downdated upper triangular matrix.
%
%    Output, real Z(LDZ,NZ), the downdated vectors.
%
%    Output, real RHO(NZ), the norms of the residual vectors,
%    which have been changed along with R and Z.
%
%    Output, real C(P), S(P), the cosines and sines of the
%    transforming rotations.
%
%    Output, integer INFO, return flag.
%     0, the entire downdating was successful.
%    -1, if R could not be downdated.  In this case, all quantities
%        are left unaltered.
%     1, if some RHO could not be downdated.  The offending RHO's are
%        set to -1.
%

%
%  Solve R' * A = X, placing the result in the array S.
%
  info = 0;
  s(1) = x(1) / r(1,1);

  for j = 2 : p
    s(j) = x(j) - ddot ( j-1, r(1:j-1,j), 1, s(1:j-1), 1 );
    s(j) = s(j) / r(j,j);
  end

  norm = dnrm2 ( p, s(1:p), 1 );

  if ( 1.0 <= norm )
    info = -1;
    return
  end

  alpha = sqrt ( 1.0 - norm * norm );
%
%  Determine the transformations.
%
  for ii = 1 : p
    i = p - ii + 1;
    scale = alpha + abs ( s(i) );
    a = alpha / scale;
    b = s(i) / scale;
    norm = sqrt ( a * a + b * b );
    c(i) = a / norm;
    s(i) = b / norm;
    alpha = scale * norm;
  end
%
%  Apply the transformations to R.
%
  for j = 1 : p
    xx = 0.0;
    for ii = 1 : j
      i = j - ii + 1;
      t = c(i) * xx + s(i) * r(i,j);
      r(i,j) = c(i) * r(i,j) - s(i) * xx;
      xx = t;
    end
  end
%
%  If required, downdate Z and RHO.
%
  for j = 1 : nz

    zeta = y(j);
    for i = 1 : p
      z(i,j) = ( z(i,j) - s(i) * zeta ) / c(i);
      zeta = c(i) * zeta - s(i) * z(i,j);
    end

    azeta = abs ( zeta );

    if ( rho(j) < azeta )
      info = 1;
      rho(j) = -1.0;
    else
      rho(j) = rho(j) * sqrt ( 1.0 - ( azeta / rho(j) )^2 );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%*********************************************************************
%% sortA: sort columns of At{p} in ascending order according to the
%%        number of nonzero elements.
%%
%% [At,C,b,X0,Z0,permA,permZ] = sortA(blk,At,C,b,X0,Z0);
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*********************************************************************

   function [At,C,X0,Z0,permA,permZ] = sortA(blk,At,C,b,X0,Z0);

   global spdensity smallblkdim
%%
   if isempty(spdensity); spdensity = 0.4; end
   if isempty(smallblkdim); smallblkdim = 50; end
%%
   numblk = size(blk,1);
   m = length(b);
   nnzA = zeros(numblk,m);
   permA = kron(ones(numblk,1),[1:m]);
   permZ = cell(size(blk,1),1);
%%
   for p=1:size(blk,1)
      pblk = blk(p,:);
      n = sum(pblk{2});
      numblk = length(pblk{2});
      if strcmp(pblk{1},'s') & (max(pblk{2}) > smallblkdim)
         n2 = sum(pblk{2}.*pblk{2});  n22 = sum(pblk{2}.*(pblk{2}+1))/2;
         m1 = size(At{p,1},2);
         if (length(pblk{2}) == 1)
            tmp = abs(C{p}) + abs(Z0{p});
            if  (~isempty(At{p,1}))
                tmp = tmp + smat(blk(p,:),abs(At{p,1})*ones(m1,1),1);
            end
            if (nnz(tmp) < spdensity*n22);
               per = symamd(tmp);
               invper = zeros(n,1); invper(per) = [1:n];
               permZ{p} = invper;
               if (~isempty(At{p,1}))
                  isspAt = issparse(At{p,1});
                  for k = 1:m1
                     Ak = smat(pblk,At{p,1}(:,k),1);
                     At{p,1}(:,k) = svec(pblk,Ak(per,per),isspAt);
                  end
               end
               C{p}  = C{p}(per,per);
               Z0{p} = Z0{p}(per,per);
               X0{p} = X0{p}(per,per);
            else
               per = [];
            end
            if (length(pblk) > 2) & (~isempty(per))
               m2 = length(pblk{3});
               P = spconvert([(1:n)', per', ones(n,1)]);
               At{p,2} = P*At{p,2};
            end
         end
         if ~isempty(At{p,1}) & (mexnnz(At{p,1}) < m*n22/2)
            for k = 1:m1
                Ak = At{p,1}(:,k);
                nnzA(p,k) = length(find(abs(Ak) > eps));
            end
            [dummy,permAp] = sort(nnzA(p,1:m1));
            At{p,1}  = At{p,1}(:,permAp);
            permA(p,1:m1) = permAp;
         end
      elseif strcmp(pblk{1},'q') | strcmp(pblk{1},'l') | strcmp(pblk{1},'u');
         if ~issparse(At{p,1});
            At{p,1} = sparse(At{p,1});
         end
      end
   end
%%*********************************************************************

repo: REST


text:
function [tissueModel, cRes] = pruningModel(model, rankNonCore, coreRxn, zeroExpRxns, precursorMets, eta, tol)
% Prunes the reactions of the model based on their expression,
% connectivity to core and confidence score.
%
% USAGE:
%   [tissueModel, cRes] = pruningModel(model, rankNonCore, coreRxn, zeroExpRxns, precursorMets, eta, tol)
%
% INPUTS:
%   model:               input model (COBRA model structure)
%   rankNonCore:         order for reaction removal
%	coreRxn:             core reactions in model
%	zeroExpRxns:         reactions with zero expression (i.e., measured zero, not just
%                        missing from expression data)
%   precursorMets:       list of metabolites involved in the
%                        reactions defined in protected reactions
%   eta:                 tradeoff between removing core and zero-expression
%                        reactions (default value: 1/3)
%   tol:                 minimum flux threshold for "expressed" reactions
%                        (default 1e-8)
%
% OUTPUTS:
%	tissueModel:         pruned, context-specific model
%	cRes:                result of model checks (consistency/function)
%                           - vs. +: reaction r removed from generic model or
%                           not
%                           1 vs. 2: reaction r had zero or non-zero expression evidence
%                           -x.y: removal of reaction r corresponded with removal of y (num.) total
%                           core reactions
%                           +x.1 vs. x.0: precursor production possible after removal of
%                           reaction r or not
%                           3: removal of reaction r by itself prevented production of required
%                           metabolites (therefore was not removed)
%
%
% Authors: - This script is an adapted version of the implementation from
%            https://github.com/jaeddy/mcadre.
%          - Modified and commented by S. Opdam and A. Richelle,May 2017

    tissueModel = model;
    R_P = model.rxns;
    nonCoreRemoved = 0;
    coreRemoved = 0;
    cRes = zeros(3000, 1);
    count = 1;

    paramConsistency.epsilon=tol;
    paramConsistency.modeFlag=0;
    paramConsistency.method='fastcc';

    while numel(rankNonCore) > 0

        display(['Reaction no. ', num2str(count)])
        r = rankNonCore(1);
        display(['Attempting to remove reaction ', r{:}, '...'])
        modelR = removeRxns(tissueModel, r);

        % First check precursor production; if this test fails, no need to
        % check model consistency with FVA (time-saving step)
        rStatus = checkModelFunction(modelR, precursorMets);

        if rStatus

            % Check for inactive reactions after removal of r
            if numel(r)
                % Remove reaction r from the model
                model_rem = removeRxns(tissueModel, r);
            end
            % Check for inactive reactions after removal of r
            try
                [fluxConsistentMetBool,fluxConsistentRxnBool] = findFluxConsistentSubset(model_rem,paramConsistency);
                rStatus_and_not_error = true;
            catch
                rStatus_and_not_error = false;
            end
        else
            rStatus_and_not_error = false;
        end

        if rStatus_and_not_error
            inactive_G= [ r; model_rem.rxns(fluxConsistentRxnBool==0)];

            inactiveCore = intersect(inactive_G, coreRxn);
            inactiveNonCore = setdiff(inactive_G, inactiveCore);

            % Remove reactions with zero expression (previously penalized in
            % rank_reactions) and corresponding inactive core reactions, only if
            % sufficiently more non-core reactions are removed
            if ismember(r, zeroExpRxns)

                display('Zero-expression evidence for reaction...')

                % Check model function with all inactive reactions removed
                modelTmp = removeRxns(tissueModel, inactive_G);
                tmpStatus = checkModelFunction(modelTmp, precursorMets);

                if (numel(inactiveCore) / numel(inactiveNonCore) <= eta) && tmpStatus

                    R_P = setdiff(R_P, inactive_G);
                    tissueModel = removeRxns(tissueModel, inactive_G);

                    rankNonCore(ismember(rankNonCore, inactive_G)) = [];

                    nonCoreRemoved = nonCoreRemoved + numel(inactiveNonCore);
                    coreRemoved = coreRemoved + numel(inactiveCore);
                    num_removed = nonCoreRemoved + coreRemoved;

                    display('Removed all inactive reactions')

                    % result = -1.x indicates that reaction r had zero
                    % expression evidence and was removed along with any
                    % consequently inactivated reactions; x indicates the number of
                    % core reactions removed
                    if numel(inactiveCore) > 100
                        removed_C_indicator = numel(inactiveCore) / 100;
                    else
                        removed_C_indicator = numel(inactiveCore) / 10;
                    end
                    result = -1 - removed_C_indicator;

                else
                    % Note: no reactions (core or otherwise) are actually
                    % removed in this step, but it is necessary to update the
                    % total number of removed reactions to avoid errors below
                    num_removed = nonCoreRemoved + coreRemoved;
                    rankNonCore(1) = [];

                    display('No reactions removed')

                    % result = 1.x indicates that no reactions were removed
                    % because removal of r either led to a ratio of inactivated
                    % core vs. non-core reactions above the specified threshold
                    % eta (x = 1) or the removal of r and consequently
                    % inactivated reactions prevented production of required
                    % metabolites (x = 0)
                    result = 1 + tmpStatus / 10;
                end

            % If reaction has expression evidence, only attempt to remove
            % inactive non-core reactions
            else
                % Check model function with non-core inactive reactions removed
                modelTmp = removeRxns(tissueModel, inactiveNonCore);
                tmpStatus = checkModelFunction(modelTmp, precursorMets);

                if numel(inactiveCore) == 0 && tmpStatus
                    R_P = setdiff(R_P, inactiveNonCore);
                    tissueModel = removeRxns(tissueModel, inactiveNonCore);

                    rankNonCore(ismember(rankNonCore, inactiveNonCore)) = [];

                    nonCoreRemoved = nonCoreRemoved + numel(inactiveNonCore);
                    num_removed = nonCoreRemoved + coreRemoved;
                    display('Removed non-core inactive reactions')

                    % result = -2 indicates that reaction r had expression.
                    % evidence and was removed along with (only) non-core
                    % inactivated reactions; x indicates the number of
                    % core reactions removed (should be zero!)
                    if numel(inactiveCore) > 100
                        removed_C_indicator = numel(inactiveCore) / 100;
                    else
                        removed_C_indicator = numel(inactiveCore) / 10;
                    end
                    result = -2 - removed_C_indicator;

                else
                    num_removed = nonCoreRemoved + coreRemoved;
                    rankNonCore(1) = [];

                    display('No reactions removed')

                    % result = 2.x indicates that no reactions were removed
                    % because removal of r either led to inactivated core
                    % reactions (x = 1) or the removal of r and consequently
                    % inactivated reactions prevented production of required
                    % metabolites (x = 0)
                    result = 2 + tmpStatus / 10;
                end
            end
        else
            num_removed = nonCoreRemoved + coreRemoved;
            rankNonCore(1) = [];

            % result = 3 indicates that no reactions were removed because
            % removal of r by itself prevented production of required
            % metabolites
            result = 3;
        end

        cRes(count) = result;
        count = count + 1;
        display(sprintf(['Num. removed: ', num2str(num_removed), ...
            ' (', num2str(coreRemoved), ' core, ', ...
            num2str(nonCoreRemoved), ' non-core); ', ...
            'Num. remaining: ', num2str(numel(rankNonCore)), '\n']))
    end
    cRes(count:end) = [];
end

author: opencobra

repo: cobratoolbox


text:
function [dim] = pos2dim3d(pos,dimold)

% POS2DIM3D reconstructs the volumetric dimensions from an ordered list of
% positions. optionally, the original dim can be provided, and the (2:end)
% elements are appended to the output.
%
% Use as
%   [dim] = pos2dim3d(pos, dimold)
% where pos is an ordered list of positions and where the (optional)
% dimold is a vector with the original dimensionality of the anatomical
% or functional data.
%
% The output dim is a 1x3 or 1xN vector of which the first three elements
% correspond to the 3D volumetric dimensions.
%
% See also POS2DIM, POS2TRANSFORM

% Copyright (C) 2009, Jan-Mathijs Schoffelen

if nargin==1 && ~isstruct(pos),
  dimold = zeros(0,2);
elseif isstruct(pos),
  % the input is a FieldTrip data structure
  dimord = pos.dimord;
  dimtok = tokenize(dimord, '_');
  for i = 1:length(dimtok)
    if strcmp(dimtok{i},'pos'),
      dimold(i,1) = size(pos.pos,1);
    else
      dimold(i,1) = numel(getfield(pos, dimtok{i}));
    end
  end
  pos = pos.pos;
else
  if size(pos,1)~=dimold(1),
    ft_error('the first element in the second input should be equal to the number of positions');
  end
end

% extract the dim now that the bookkeeping is done
dim = pos2dim(pos);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Gamma = simgamma_ltd(T,P,Pi,rate,L,refractory_period,grouping)
%
% Simulate state time courses with longer-than-order-1-markov time dependencies
%
% INPUTS:
%
% T                     Number of time points for each time series
% P                     Transition probability matrix (K by K)
% Pi                    Initial probabilities (K by 1)
% rate                  The weights that model the contribution of the
%                       latest L points are modelled by a Gamma
%                       distribution, whose shape is 1 - and rate is
%                       specified here
% L                     The length of the history (in number of time points)
%                       that influences the state at time t
% refractory_period     to prevent bursts of quick changes, refractory_period
%                       can be set so that, after a change, you cannot change again
%                       after 'refractory_period' number of iterations
%
% OUTPUTS
%
% Gamma         simulated  p(state | data)
%
% Author: Diego Vidaurre, OHBA, University of Oxford

N = length(T); K = length(Pi);

if nargin<4, rate = 2; end
if nargin<5, L = 10; end
if nargin<6, refractory_period = 2; end
if nargin<7, grouping = []; end

Gamma = zeros(sum(T),K);
weights = gampdf(0:L-1,1,rate)';
weights = weights / sum(weights);
weights = weights(end:-1:1);

for n = 1:N
    if ~isempty(grouping)
        i = grouping(n);
        Pn = P(:,:,i); Pin = Pi(:,i)';
    else
        Pn = P; Pin = Pi;
    end
    Gammai = zeros(T(n),K);
    if any(Pin==1)
        Gammai(1,Pin==1) = 1;
    else
        Gammai(1,:) = mnrnd(1,Pin);
    end
    last_ch = Inf;
    for t = 2:L
        if last_ch < refractory_period
            Gammai(t,:) = Gammai(t-1,:);
            last_ch = last_ch + 1;
        else
            if t==2
                g = repmat(weights(end-t+2:end,1),1,K) .* Gammai(1:t-1,:);
            else
                g = sum(repmat(weights(end-t+2:end,1),1,K) .* Gammai(1:t-1,:));
            end
            g = g / sum(g);
            Gammai(t,:) = mnrnd(1,g*Pn);
            if any(Gammai(t,:)~=Gammai(t-1,:)), last_ch = 1; end
            Gammai(t,:) = Gammai(t,:) / sum(Gammai(t,:));
        end
    end
    for t=L+1:T(n)
        if last_ch < refractory_period
            Gammai(t,:) = Gammai(t-1,:);
            last_ch = last_ch + 1;
        else
            g = sum(repmat(weights,1,K) .* Gammai(t-L:t-1,:));
            Gammai(t,:) = mnrnd(1,g*Pn);
            if any(Gammai(t,:)~=Gammai(t-1,:)), last_ch = 1; end
        end
    end
    t = (1:T(n)) + sum(T(1:n-1));
    Gamma(t,:) = Gammai;
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function [ x, fx ] = local_min ( a, b, epsi, t, f, x )

%*****************************************************************************80
%
%% LOCAL_MIN seeks a local minimum of a function F(X) in an interval [A,B].
%
%  Discussion:
%
%    The method used is a combination of golden section search and
%    successive parabolic interpolation.  Convergence is never much slower
%    than that for a Fibonacci search.  If F has a continuous second
%    derivative which is positive at the minimum (which is not at A or
%    B), then convergence is superlinear, and usually of the order of
%    about 1.324....
%
%    The values EPSI and T define a tolerance TOL = EPSI * abs ( X ) + T.
%    F is never evaluated at two points closer than TOL.
%
%    If F is a unimodal function and the computed values of F are always
%    unimodal when separated by at least SQEPS * abs ( X ) + (T/3), then
%    LOCAL_MIN approximates the abscissa of the global minimum of F on the
%    interval [A,B] with an error less than 3*SQEPS*abs(LOCAL_MIN)+T.
%
%    If F is not unimodal, then LOCAL_MIN may approximate a local, but
%    perhaps non-global, minimum to the same accuracy.
%
%    Thanks to Jonathan Eggleston for pointing out a correction to the
%    golden section step, 01 July 2013.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2008
%
%  Author:
%
%    Original FORTRAN77 version by Richard Brent.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Richard Brent,
%    Algorithms for Minimization Without Derivatives,
%    Dover, 2002,
%    ISBN: 0-486-41998-3,
%    LC: QA402.5.B74.
%
%  Parameters:
%
%    Input, real A, B, the endpoints of the interval.
%
%    Input, real EPSI, a positive relative error tolerance.
%    EPSI should be no smaller than twice the relative machine precision,
%    and preferably not much less than the square root of the relative
%    machine precision.
%
%    Input, real T, a positive absolute error tolerance.
%
%    Input, function value = F ( x ), the name of a user-supplied
%    function whose local minimum is being sought.
%
%    Output, real X, the estimated value of an abscissa
%    for which F attains a local minimum value in [A,B].
%
%    Output, real FX, the value F(X).
%

%
%  C is the square of the inverse of the golden ratio.
%
  c = 0.5 * ( 3.0 - sqrt ( 5.0 ) );

  sa = a;
  sb = b;
  x = sa + c * ( b - a );
  w = x;
  v = w;
  e = 0.0;
  fx = f ( x );
  fw = fx;
  fv = fw;

  while ( 1 )

    m = 0.5 * ( sa + sb );
    tol = epsi * abs ( x ) + t;
    t2 = 2.0 * tol;
%
%  Check the stopping criterion.
%
    if ( abs ( x - m ) <= t2 - 0.5 * ( sb - sa ) )
      break
    end
%
%  Fit a parabola.
%
    r = 0.0;
    q = r;
    p = q;

    if ( tol < abs ( e ) )

      r = ( x - w ) * ( fx - fv );
      q = ( x - v ) * ( fx - fw );
      p = ( x - v ) * q - ( x - w ) * r;
      q = 2.0 * ( q - r );

      if ( 0.0 < q )
        p = - p;
      end

      q = abs ( q );

      r = e;
      e = d;

    end

    if ( abs ( p ) < abs ( 0.5 * q * r ) & ...
         q * ( sa - x ) < p & ...
         p < q * ( sb - x ) )
%
%  Take the parabolic interpolation step.
%
      d = p / q;
      u = x + d;
%
%  F must not be evaluated too close to A or B.
%
      if ( ( u - sa ) < t2 | ( sb - u ) < t2 )

        if ( x < m )
          d = tol;
        else
          d = - tol;
        end

      end
%
%  A golden-section step.
%
    else

      if ( x < m )
        e = sb - x;
      else
        e = sa - x;
      end

      d = c * e;

    end
%
%  F must not be evaluated too close to X.
%
    if ( tol <= abs ( d ) )
      u = x + d;
    elseif ( 0.0 < d )
      u = x + tol;
    else
      u = x - tol;
    end

    fu = f ( u );
%
%  Update A, B, V, W, and X.
%
    if ( fu <= fx )

      if ( u < x )
        sb = x;
      else
        sa = x;
      end

      v = w;
      fv = fw;
      w = x;
      fw = fx;
      x = u;
      fx = fu;

    else

      if ( u < x )
        sa = u;
      else
        sb = u;
      end

      if ( fu <= fw | w == x )
        v = w;
        fv = fw;
        w = u;
        fw = fu;
      elseif ( fu <= fv | v == x | v == w )
        v = u;
        fv = fu;
      end

    end

  end

  return
end


text:
function [xx, yy, zz] = chebpts3(nx, ny, nz, dom, kind)
%CHEBPTS3   3D tensor product Chebyshev grid.
%   [XX YY ZZ] = CHEBPTS3(N) constructs an N by N by N grid of Chebyshev
%   tensor points on [-1 1]^3.
%
%   [XX YY ZZ] = CHEBPTS3(NX,NY,NZ) constructs an NX by NY by NZ grid of
%   Chebyshev tensor points on [-1 1]^3.
%
%   [XX YY ZZ] = CHEBPTS3(NX,NY,NZ,DOM) constructs an NX by NY by NZ grid
%   of Chebyshev tensor points on the cube [a b] x [c d] x [e g], where
%   DOM = [a b c d e g].
%
%   [XX YY ZZ] = CHEBPTS3(NX,NY,NZ,D,KIND) constructor Chebyshev tensor
%   grid of the kind KIND. KIND = 2 is default.
%
% See also CHEBPTS and CHEBPTS2.

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

if ( nargin > 3 )
   % Third argument should be a domain.
   dom = dom(:).';  % make a row vector.
   if ( ~all(size(dom) == [1 6]) )
        error('CHEBFUN:chebpts3:domain', 'Unrecognised domain.');
   end
else  % Default to the canoncial domain.
    dom = [-1, 1, -1, 1, -1, 1];
end

if ( nargin == 1 )
    % Make it a square Chebyshev grid if only one input.
    ny = nx;
    nz = nx;
end
if ( nargin < 5 )
    % default to Chebyshev points of the 2nd kind.
    kind = 2;
end

% Get points:
x = chebpts(nx, dom(1:2), kind);
y = chebpts(ny, dom(3:4), kind);
z = chebpts(nz, dom(5:6), kind);

% Tensor product.
[xx, yy, zz] = ndgrid(x, y, z);

end


text:
function pass = test_jumps_manual(pref)

if ( nargin == 0 )
    pref = cheboppref();
end

% BVP with nonhomogeneous jump condition
% Toby Driscoll, June 2014
% Based on issue #665

tol = 1e1*pref.bvpTol;

%%
N = chebop(@(s, V) diff(V, 2), [-1, 1]);
N.lbc = @(V) V + 1;
N.rbc = @(V) V - 1;
N.bc = @(s,V) [jump(V, 0) - 1 ; jump(diff(V),0)];
y = N\0;

%%
err(1) = norm(N(y), inf);
jump1  = feval(y, 0, 'right') - feval(y, 0, 'left');
err(2) = abs(jump1 - 1);
jump2  = feval(diff(y), 0, 'right') - feval(diff(y), 0, 'left');
err(3) = abs(jump2);
err(4) = abs(y(-1) + 1);
err(5) = abs(y(1) - 1);

%%
pass = err < tol;

end


text:
function [pca_dim, k_var, U, S, V] = PCA_MP_fit(M, pvar, flag_plot)
% assume M is already preprocessed
if ~exist('flag_plot', 'var')
    flag_plot = 1;
end
if ~exist('pvar', 'var')
    pvar = 0.9;
end

[U, S, V] = svd(M, 'econ');
s = diag(S).^2;

ep = 1e-4;
splot = s(abs(s)>ep);
ns = length(splot);

gap = 1; % only for the moment estimator
pout = 0.2;
niter = 10;
nx_target = 3.2;
[lb_e, sigma_e, tfconv] = estimate_MP2(splot, pout, gap, niter, nx_target, 'L2');
lb1_e = (1 + sqrt(lb_e))^2 * sigma_e^2;
ns_bulk_e = sum(splot < lb1_e);
pca_dim = sum(splot > lb1_e);
s_out = splot(splot > lb1_e);
if isempty(s_out)
    s_out = max(splot);
    pca_dim = 1;
    k_var = 1;
else
    s_out = sort(s_out, 'descend');
    var_sum = cumsum(s_out);
    k_var = find(var_sum >= pvar * var_sum(end), 1);
end

if flag_plot ==1
    figure;
    set(gcf,'position', [88, 571, 656, 534]);
    dbin = adaptBinWidth2(splot, 4);
    nbins_all = round(max(splot)/dbin);
    ht = histogram(splot, nbins_all, 'Normalization', 'pdf');
    ns_plot = ns;
    % histogram(splot, 40000);
    % h=histogram(splot(splot< 2*lb1), 80, 'Normalization', 'pdf');
    % xlim([0,15])
    hold on;
    x0 = linspace(0, 2*lb1_e, 1000);
    y0 = MPdistr(x0, lb_e, sigma_e);
    y0 = y0 / ns_plot * ns_bulk_e;
    plot(x0, y0, 'r')
    xlim([0, 5*lb1_e])
    hold off;
    title(['Estimated signal dimension = ', num2str(pca_dim)])
    length(splot)
    % display('Dim  lb  sigma  flag_conv')
    pca_dim, lb_e, sigma_e, tfconv

    figure;
    set(gcf,'position', [88, 571, 656, 534]);
    plot(var_sum / var_sum(end),'bo-');
    title('Explained variance');
    legend([num2str(pvar*100),'% = ', num2str(k_var)])
end

end

author: xiuyechen

repo: FishExplorer


text:
%  Figure 3.4      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
% script to generate Fig. 3.4
%% fig3_04.m    Example 3.5
clf;
k=1;
num=1;                        % form numerator
den=[1 k];                    % form denominator
% sinusoidal input signal
deltaT = 0.001;
t=0:deltaT:10;                % form time vector
u=sin(10*(t));                % form input
sys=tf(num,den);              % form system
[y]=lsim(sys,u,t);            % linear simulation
% plot response
figure();
plot(t,y);
xlabel('Time (sec)');
ylabel('Output');
title('Fig. 3.4 (a): transient response');
pause;
hold on;
y1=(10/101)*exp(-t);
phi=atan(-10);
y2=(1/sqrt(101))*sin(10*t+phi);
plot(t,y1,t,y2,t,y1+y2);
% grid
nicegrid
hold off;
pause;
figure();
ii=[9001:10001];
plot(t(ii),y(ii),t(ii),u(ii));
xlabel('Time (sec)');
ylabel('Output, input');
title('Fig. 3.4 (b): Steady-state response');
text(9.4,0.65,'u(t)');
text(9.24,0.12,'y(t)');
% grid
nicegrid


text:
function d=dsphere(p,xc,yc,zc,r)

%   Copyright (C) 2004-2006 Per-Olof Persson. See COPYRIGHT.TXT for details.

d=sqrt((p(:,1)-xc).^2+(p(:,2)-yc).^2+(p(:,3)-zc).^2)-r;

repo: ifem


text:
%% string gen_file_name(string,string)
% fileName - string
% extra - string
% newFileName - string
%
function [newFileName] = gen_file_name(fileName,extra)
  %ext = get_file_extension(fileName);
  %newFileName = [fileName(1:length(fileName)-4) extra ext];
  [pathstr,name,ext] = fileparts(fileName);
  newFileName = fullfile(pathstr,[name extra ext]);
end

author: andrewssobral

repo: lrslibrary


text:
base_graph_list = 1:2;

Z_c_list = [2, 4, 8, 16, 32, 64, 128, 256, ...
3, 6, 12, 24, 48, 96, 192, 384, ...
5, 10, 20, 40, 80, 160, 320, ...
7, 14, 28, 56, 112, 224, ...
9, 18, 36, 72, 144, 288, ...
11, 22, 44, 88, 176, 352, ...
13, 26, 52, 104, 208, ...
15, 30, 60, 120, 240];

pass_or_failure = 1;
for base_graph = base_graph_list
    for Z_c = Z_c_list
        pass_flag = test_ldpc(base_graph, Z_c);
        if ~pass_flag
          pass_or_failure = 0;
          disp(['The case base_graph = ', num2str(base_graph), ' Z_c= ', num2str(Z_c), ' passed.']);
          break;
        else
          disp(['case base_graph = ', num2str(base_graph), ' Z_c= ', num2str(Z_c), ' passed.']);
        end
    end
end

if pass_or_failure == 1
   disp('test passed.');
else
   disp('test failed.');
end

author: xiaoshaoning


text:
function DCM = spm_dcm_erp_dipfit(DCM, save_vol_sens)
% Prepare structures for ECD forward model (EEG, MEG and LFP)
% FORMAT DCM = spm_dcm_erp_dipfit(DCM, save_vol_sens)
% DCM           - DCM structure
% save_vol_sens - optional argument indicating whether to perform
%                 the time consuming step required for actually using
%                 the forward model to compute lead fields (1, default)
%                 or skip it if the function is only called for
%                 verification of the input (0).
%
% Input DCM structure requires:
%       DCM.xY.Dfile
%       DCM.xY.Ic
%       DCM.Lpos
%       DCM.options.spatial - 'ERP', 'LFP' or 'IMG'
%
% fills in:
%
%       DCM.M.dipfit
%
%    dipfit.location - 0 or 1 for source location priors
%    dipfit.symmetry - 0 or 1 for symmetry constraints on sources
%    dipfit.modality - 'EEG', 'MEG', 'MEGPLANAR' or 'LFP'
%    dipfit.type     - 'ECD', 'LFP' or 'IMG''
%    dipfit.symm     - distance (mm) for symmetry constraints (ECD)
%    dipfit.Lpos     - x,y,z source positions (mm)            (ECD)
%    dipfit.Nm       - number of modes                        (Imaging)
%    dipfit.Ns       - number of sources
%    dipfit.Nc       - number of channels
%
%    dipfit.vol      - volume structure (for M/EEG)
%    dipfit.datareg  - registration structure (for M/EEG)
%__________________________________________________________________________
% Copyright (C) 2007-2015 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_erp_dipfit.m 7702 2019-11-22 11:32:26Z guillaume $

% Get data filename and good channels
%--------------------------------------------------------------------------
try
    DCM.xY.Dfile;
    DCM.M.dipfit.Ic = DCM.xY.Ic;
catch
    spm('alert*','Please specify data');
    error('Data not specified.')
end

if nargin == 1
  save_vol_sens = 0;
end

% D - SPM data structure
%--------------------------------------------------------------------------
try
    D  = spm_eeg_load(DCM.xY.Dfile);
catch
    D  = spm_eeg_load(spm_file(DCM.xY.Dfile,'filename'));
end

% set options in dipfit
%--------------------------------------------------------------------------
try, spatial  = DCM.options.spatial;  catch, spatial  = 'IMG'; end
try, location = DCM.options.location; catch, location = 0;     end
try, symmetry = DCM.options.symmetry; catch, symmetry = 0;     end

DCM.M.dipfit.type     = spatial;
DCM.M.dipfit.location = location;
DCM.M.dipfit.symmetry = symmetry;

% Get source locations if MEG or EEG
%--------------------------------------------------------------------------
switch DCM.xY.modality

    % get source priors for EEG or MEG
    %----------------------------------------------------------------------
    case{'EEG','MEG','MEGPLANAR'}

        [D, ok] = check(D, '3d');
        if ~ok
            spm('alert*',{'File not ready for source reconstruction.',...
                'Use Prepare to specify sensors and fiducials.'});
        end

        try
            DCM.M.dipfit.Lpos = DCM.Lpos;
        catch
            spm('alert*',{'Please specify source locations','in DCM.Lpos'})
        end

        DCM.M.dipfit.modality = DCM.xY.modality;
        DCM.M.dipfit.Ns       = length(DCM.Sname);
        DCM.M.dipfit.Nc       = length(DCM.xY.Ic);

        % otherwise assume LFP
        %------------------------------------------------------------------
    otherwise

        DCM.M.dipfit.modality = 'LFP';
        DCM.M.dipfit.Ns       = length(DCM.Sname);
        DCM.M.dipfit.Nc       = length(DCM.xY.Ic);
        return
end

% Detect silent sources for non-LFP type from GUI source name text
%--------------------------------------------------------------------------
DCM.M.dipfit.silent_source  = strfind(DCM.Sname,'silent');

% If not LFP, get electromagnetic forward model
%==========================================================================
if ~isfield(D, 'val'), D.val = 1; end

% detects old version of the struct
%--------------------------------------------------------------------------
if ~isfield(D, 'inv') || ~iscell(D.inv) ||...
        ~(isfield(D.inv{D.val}, 'forward') && isfield(D.inv{D.val}, 'datareg')) ||...
        ~isa(D.inv{D.val}.mesh.tess_ctx, 'char')
    D = spm_eeg_inv_mesh_ui(D, D.val);
    D = spm_eeg_inv_datareg_ui(D, D.val);
    D = spm_eeg_inv_forward_ui(D, D.val);
    save(D);
end

% fill in dipfit
%--------------------------------------------------------------------------
for m = 1:numel(D.inv{D.val}.forward)
    if strncmp(DCM.xY.modality, D.inv{D.val}.forward(m).modality, 3)
        DCM.M.dipfit.vol      = D.inv{D.val}.forward(m).vol;
        DCM.M.dipfit.datareg  = D.inv{D.val}.datareg(m);
        if isfield(D.inv{D.val}.forward(m), 'siunits') && D.inv{D.val}.forward(m).siunits
            DCM.M.dipfit.sens     = D.inv{D.val}.forward(m).sensors;
            DCM.M.dipfit.siunits  = true;
        else
            DCM.M.dipfit.sens     = DCM.M.dipfit.datareg.sensors;
            DCM.M.dipfit.siunits  = false;
        end
    end
end

% channels
%--------------------------------------------------------------------------
if save_vol_sens
    if ischar(DCM.M.dipfit.vol)
        DCM.M.dipfit.vol = ft_read_headmodel(DCM.M.dipfit.vol);
    end

    [DCM.M.dipfit.vol, DCM.M.dipfit.sens] = ft_prepare_vol_sens(DCM.M.dipfit.vol, ...
        DCM.M.dipfit.sens, 'channel', D.chanlabels(DCM.xY.Ic));
end

switch DCM.options.spatial

    % Imaging (distributed source reconstruction)
    %----------------------------------------------------------------------
    case{'IMG'}

        % Load Gain or Lead field matrix
        %------------------------------------------------------------------
        DCM.val = D.val;
        [L,D]   = spm_eeg_lgainmat(D, [], D.chanlabels(DCM.xY.Ic));

        % centers
        %------------------------------------------------------------------
        xyz = DCM.M.dipfit.Lpos;
        Np  = size(xyz,2);

        % parameters
        %==================================================================

        % defaults: Nm = 6; number of modes per region
        %------------------------------------------------------------------
        try, rad  = DCM.M.dipfit.radius; catch, rad  = 16;    end
        try, Nm   = DCM.M.dipfit.Nm;     catch, Nm   = 6;     end

        % Compute spatial basis (eigenmodes of lead field)
        %==================================================================

        % create MSP spatial basis set in source space
        %------------------------------------------------------------------
        vert   = D.inv{D.val}.mesh.tess_mni.vert;
        for i  = 1:Np
            Dp = sum([vert(:,1) - xyz(1,i), ...
                      vert(:,2) - xyz(2,i), ...
                      vert(:,3) - xyz(3,i)].^2,2);

            % nearest mesh points
            %--------------------------------------------------------------
            Ip = find(Dp < rad^2);
            if length(Ip) < Nm;
                [y,Ip] = sort(Dp);
                Ip     = Ip(1:Nm);
            end

            % left hemisphere
            %--------------------------------------------------------------
            U                  = spm_svd(L(:,Ip)',0);
            U                  = U(:,1:Nm);
            DCM.M.dipfit.G{i}  = L(:,Ip)*U;
            DCM.M.dipfit.U{i}  = U;
            DCM.M.dipfit.Ip{i} = Ip;
        end

        % Save results
        %==================================================================
        DCM.M.dipfit.radius  = rad;                   % VOI (XYZ, Radius)
        DCM.M.dipfit.Nm      = Nm;                    % modes per region
        DCM.M.dipfit.Nd      = length(vert);          % number of dipoles
        DCM.M.dipfit.gainmat = D.inv{D.val}.gainmat;  % Lead field filename

    otherwise

end

author: spm


text:
function cost = FDDL_fidelity(Y, Y_range, D, D_range, X)

    nClasses = numel(Y_range) - 1;
    cost = 0;
    for c = 1: nClasses
        Yc = get_block_col(Y, c, Y_range);
        Dc = get_block_col(D, c, D_range);
        Xc = get_block_row(X, c, D_range);
        Xcc = get_block_col(Xc, c, Y_range);
        cost = cost + normF2(Yc - Dc *Xcc);
        for j = 1:nClasses
            if j == c
                continue;
            else
                Xcj = get_block_col(Xc, j, Y_range);
                cost = cost + normF2(Dc*Xcj);
            end
        end
    end

end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL