text:
function figure2dot(filename, varargin)
%FIGURE2DOT    Save figure in Graphviz (.dot) file.
%  FIGURE2DOT(filename) saves the current figure as dot-file.
%
%  FIGURE2DOT(filename, 'object', HGOBJECT) constructs the graph representation
%  of the specified object (default: gcf)
%
%  You can visualize the constructed DOT file using:
%    - [GraphViz](http://www.graphviz.org) on your computer
%    - [WebGraphViz](http://www.webgraphviz.com) online
%    - [Gravizo](http://www.gravizo.com) for your markdown files
%    - and a lot of other software such as OmniGraffle
%

%   Copyright (c) 2008--2014, Nico Schlömer <nico.schloemer@gmail.com>
%   All rights reserved.
%
%   Redistribution and use in source and binary forms, with or without
%   modification, are permitted provided that the following conditions are
%   met:
%
%      * Redistributions of source code must retain the above copyright
%        notice, this list of conditions and the following disclaimer.
%      * Redistributions in binary form must reproduce the above copyright
%        notice, this list of conditions and the following disclaimer in
%        the documentation and/or other materials provided with the distribution
%
%   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%   AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%   IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%   ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%   LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%   CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%   SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%   INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%   CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%   ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%   POSSIBILITY OF SUCH DAMAGE.
%
% =========================================================================

    ipp = m2tInputParser();
    ipp = ipp.addRequired(ipp, 'filename', @ischar);
    ipp = ipp.addParamValue(ipp, 'object', gcf, @ishghandle);
    ipp = ipp.parse(ipp, filename, varargin{:});
    args = ipp.Results;

    filehandle = fopen(args.filename, 'w');
    finally_fclose_filehandle = onCleanup(@() fclose(filehandle));

    % start printing
    fprintf(filehandle, 'digraph simple_hierarchy {\n\n');
    fprintf(filehandle, 'node[shape=box];\n\n');

    % define the root node
    node_number = 0;
    p = get(args.object, 'Parent');
    % define root element
    type = get(p, 'Type');
    fprintf(filehandle, 'N%d [label="%s"]\n\n', node_number, type);

    % start recursion
    plot_children(filehandle, p, node_number);

    % finish off
    fprintf(filehandle, '}');

    % ----------------------------------------------------------------------------
    function plot_children(fh, h, parent_node)

        children = allchild(h);

        for h = children(:)'
            if shouldSkip(h), continue, end;
            node_number = node_number + 1;

            label = {};
            label = addHGProperty(label, h, 'Type', '');
            try
                hClass = class(handle(h));
                label = addProperty(label, 'Class', hClass);
            catch
                % don't do anything
            end
            label = addProperty(label, 'Handle', sprintf('%g', double(h)));
            label = addHGProperty(label, h, 'Title', '');
            label = addHGProperty(label, h, 'Axes', []);
            label = addHGProperty(label, h, 'String', '');
            label = addHGProperty(label, h, 'Tag', '');
            label = addHGProperty(label, h, 'DisplayName', '');
            label = addHGProperty(label, h, 'Visible', 'on');
            label = addHGProperty(label, h, 'HandleVisibility', 'on');

            % print node
            fprintf(fh, 'N%d [label="%s"]\n', ...
                    node_number, m2tstrjoin(label, '\n'));

            % connect to the child
            fprintf(fh, 'N%d -> N%d;\n\n', parent_node, node_number);

            % recurse
            plot_children(fh, h, node_number);
        end
    end
end
% ==============================================================================
function bool = shouldSkip(h)
    %  returns TRUE for objects that can be skipped
    objType = get(h, 'Type');
    bool = ismember(lower(objType), guitypes());
end
% ==============================================================================
function label = addHGProperty(label, h, propName, default)
    % get a HG property and assign it to a GraphViz node label
    if ~exist('default','var') || isempty(default)
        shouldOmit = @isempty;
    elseif isa(default, 'function_handle')
        shouldOmit = default;
    else
        shouldOmit = @(v) isequal(v,default);
    end

    if isprop(h, propName)
        propValue = get(h, propName);
        if numel(propValue) == 1 && ishghandle(propValue) && isprop(propValue, 'String')
            % dereference Titles, labels, ...
            propValue = get(propValue, 'String');
        elseif ishghandle(propValue)
            % dereference other HG objects to their raw handle value (double)
            propValue = double(propValue);
        elseif iscell(propValue)
            propValue = ['{' m2tstrjoin(propValue,',') '}'];
        end

        if ~shouldOmit(propValue)
            label = addProperty(label, propName, propValue);
        end
    end
end
function label = addProperty(label, propName, propValue)
    % add a property to a GraphViz node label
    if isnumeric(propValue)
        propValue = num2str(propValue);
    elseif iscell(propValue)
        propValue = m2tstrjoin(propValue,sprintf('\n'));
    end
    label = [label, sprintf('%s: %s', propName, propValue)];
end
% ==============================================================================

author: thomas-koehler

repo: SupER


text:
function test_issue623

% MEM 2gb
% WALLTIME 00:20:00
% DEPENCENCY

%%

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/original/eeg/brainvision/Mischa.vhdr');

ft_hastoolbox('compat/obsolete', 1); % ft_qualitycheck has moved to be obsolete

%%

cfg = [];
cfg.dataset = filename;
cfg.saveplot = 'no';
cfg.savemat = 'no';
[info, timelock, freq, summary] = ft_qualitycheck(cfg);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [label, center, bCon, sumD, D] = litekmeans(X, k, varargin)
%LITEKMEANS K-means clustering, accelerated by matlab matrix operations.
%
%   label = LITEKMEANS(X, K) partitions the points in the N-by-P data matrix
%   X into K clusters.  This partition minimizes the sum, over all
%   clusters, of the within-cluster sums of point-to-cluster-centroid
%   distances.  Rows of X correspond to points, columns correspond to
%   variables.  KMEANS returns an N-by-1 vector label containing the
%   cluster indices of each point.
%
%   [label, center] = LITEKMEANS(X, K) returns the K cluster centroid
%   locations in the K-by-P matrix center.
%
%   [label, center, bCon] = LITEKMEANS(X, K) returns the bool value bCon to
%   indicate whether the iteration is converged.
%
%   [label, center, bCon, SUMD] = LITEKMEANS(X, K) returns the
%   within-cluster sums of point-to-centroid distances in the 1-by-K vector
%   sumD.
%
%   [label, center, bCon, SUMD, D] = LITEKMEANS(X, K) returns
%   distances from each point to every centroid in the N-by-K matrix D.
%
%   [ ... ] = LITEKMEANS(..., 'PARAM1',val1, 'PARAM2',val2, ...) specifies
%   optional parameter name/value pairs to control the iterative algorithm
%   used by KMEANS.  Parameters are:
%
%   'Distance' - Distance measure, in P-dimensional space, that KMEANS
%      should minimize with respect to.  Choices are:
%            {'sqEuclidean'} - Squared Euclidean distance (the default)
%             'cosine'       - One minus the cosine of the included angle
%                              between points (treated as vectors). Each
%                              row of X SHOULD be normalized to unit. If
%                              the intial center matrix is provided, it
%                              SHOULD also be normalized.
%
%   'Start' - Method used to choose initial cluster centroid positions,
%      sometimes known as "seeds".  Choices are:
%         {'sample'}  - Select K observations from X at random (the default)
%           matrix   - A K-by-P matrix of starting locations; or a K-by-1
%                      indicate vector indicating which K points in X
%                      should be used as the initial center.  In this case,
%                      you can pass in [] for K, and KMEANS infers K from
%                      the first dimension of the matrix.
%
%   'MaxIter'  - Maximum number of iterations allowed.  Default is 100.
%
%   'Replicates' - Number of times to repeat the clustering, each with a
%                  new set of initial centroids. Default is 1. If the
%                  initial centroids are provided, the replicate will be
%                  automatically set to be 1.
%
%
%
%    Examples:
%
%       fea = rand(500,10);
%       [label, center] = litekmeans(fea, 5, 'MaxIter', 50);
%
%       fea = rand(500,10);
%       [label, center] = litekmeans(fea, 5, 'MaxIter', 50, 'Replicates', 10);
%
%       fea = rand(500,10);
%       [label, center, bCon, sumD, D] = litekmeans(fea, 5, 'MaxIter', 50);
%       TSD = sum(sumD);
%
%       fea = rand(500,10);
%       initcenter = rand(5,10);
%       [label, center] = litekmeans(fea, 5, 'MaxIter', 50, 'Start', initcenter);
%
%       fea = rand(500,10);
%       idx=randperm(500);
%       [label, center] = litekmeans(fea, 5, 'MaxIter', 50, 'Start', idx(1:5));
%
%
%   See also KMEANS
%
%   version 2.0 --December/2011
%   version 1.0 --November/2011
%
%   Written by Deng Cai (dengcai AT gmail.com)

if nargin < 2
    error('litekmeans:TooFewInputs','At least two input arguments required.');
end

[n, p] = size(X);

pnames = {   'distance' 'start'   'maxiter'  'replicates' 'onlinephase'};
dflts =  {'sqeuclidean' 'sample'       []        []      'off'  };
[eid,errmsg,distance,start,maxit,reps,online] = getargs(pnames, dflts, varargin{:});
if ~isempty(eid)
    error(sprintf('litekmeans:%s',eid),errmsg);
end

if ischar(distance)
    distNames = {'sqeuclidean','cosine'};
    j = strcmpi(distance, distNames);
    j = find(j);
    if length(j) > 1
        error('litekmeans:AmbiguousDistance', ...
            'Ambiguous ''Distance'' parameter value:  %s.', distance);
    elseif isempty(j)
        error('litekmeans:UnknownDistance', ...
            'Unknown ''Distance'' parameter value:  %s.', distance);
    end
    distance = distNames{j};
else
    error('litekmeans:InvalidDistance', ...
        'The ''Distance'' parameter value must be a string.');
end

center = [];
if ischar(start)
    startNames = {'sample','cluster'};
    j = find(strncmpi(start,startNames,length(start)));
    if length(j) > 1
        error(message('litekmeans:AmbiguousStart', start));
    elseif isempty(j)
        error(message('litekmeans:UnknownStart', start));
    elseif isempty(k)
        error('litekmeans:MissingK', ...
            'You must specify the number of clusters, K.');
    end
    if j == 2
        if floor(.1*n) < 5*k
            j = 1;
        end
    end
    start = startNames{j};
elseif isnumeric(start)
    if size(start,2) == p
        center = start;
    elseif (size(start,2) == 1 || size(start,1) == 1)
        center = X(start,:);
    else
        error('litekmeans:MisshapedStart', ...
            'The ''Start'' matrix must have the same number of columns as X.');
    end
    if isempty(k)
        k = size(center,1);
    elseif (k ~= size(center,1))
        error('litekmeans:MisshapedStart', ...
            'The ''Start'' matrix must have K rows.');
    end
    start = 'numeric';
else
    error('litekmeans:InvalidStart', ...
        'The ''Start'' parameter value must be a string or a numeric matrix or array.');
end

% The maximum iteration number is default 100
if isempty(maxit)
    maxit = 100;
end

% Assume one replicate
if isempty(reps) || ~isempty(center)
    reps = 1;
end

if ~(isscalar(k) && isnumeric(k) && isreal(k) && k > 0 && (round(k)==k))
    error('litekmeans:InvalidK', ...
        'X must be a positive integer value.');
elseif n < k
    error('litekmeans:TooManyClusters', ...
        'X must have more rows than the number of clusters.');
end

bestlabel = [];
sumD = zeros(1,k);
bCon = false;

for t=1:reps
    switch start
        case 'sample'
            center = X(randsample(n,k),:);
        case 'cluster'
            Xsubset = X(randsample(n,floor(.1*n)),:);
            [dump, center] = litekmeans(Xsubset, k, varargin{:}, 'start','sample', 'replicates',1);
        case 'numeric'
    end

    last = 0;label=1;
    it=0;

    switch distance
        case 'sqeuclidean'
            while any(label ~= last) && it<maxit
                last = label;

                bb = full(sum(center.*center,2)');
                ab = full(X*center');
                D = bb(ones(1,n),:) - 2*ab;

                [val,label] = min(D,[],2); % assign samples to the nearest centers
                ll = unique(label);
                if length(ll) < k
                    %disp([num2str(k-length(ll)),' clusters dropped at iter ',num2str(it)]);
                    missCluster = 1:k;
                    missCluster(ll) = [];
                    missNum = length(missCluster);

                    aa = sum(X.*X,2);
                    val = aa + val;
                    [dump,idx] = sort(val,1,'descend');
                    label(idx(1:missNum)) = missCluster;
                end
                E = sparse(1:n,label,1,n,k,n);  % transform label into indicator matrix
                center = full((E*spdiags(1./sum(E,1)',0,k,k))'*X);    % compute center of each cluster
                it=it+1;
            end
            if it<maxit
                bCon = true;
            end
            if isempty(bestlabel)
                bestlabel = label;
                bestcenter = center;
                if reps>1
                    if it>=maxit
                        aa = full(sum(X.*X,2));
                        bb = full(sum(center.*center,2));
                        ab = full(X*center');
                        D = bsxfun(@plus,aa,bb') - 2*ab;
                        D(D<0) = 0;
                    else
                        aa = full(sum(X.*X,2));
                        D = aa(:,ones(1,k)) + D;
                        D(D<0) = 0;
                    end
                    D = sqrt(D);
                    for j = 1:k
                        sumD(j) = sum(D(label==j,j));
                    end
                    bestsumD = sumD;
                    bestD = D;
                end
            else
                if it>=maxit
                    aa = full(sum(X.*X,2));
                    bb = full(sum(center.*center,2));
                    ab = full(X*center');
                    D = bsxfun(@plus,aa,bb') - 2*ab;
                    D(D<0) = 0;
                else
                    aa = full(sum(X.*X,2));
                    D = aa(:,ones(1,k)) + D;
                    D(D<0) = 0;
                end
                D = sqrt(D);
                for j = 1:k
                    sumD(j) = sum(D(label==j,j));
                end
                if sum(sumD) < sum(bestsumD)
                    bestlabel = label;
                    bestcenter = center;
                    bestsumD = sumD;
                    bestD = D;
                end
            end
        case 'cosine'
            while any(label ~= last) && it<maxit
                last = label;
                W=full(X*center');
                [val,label] = max(W,[],2); % assign samples to the nearest centers
                ll = unique(label);
                if length(ll) < k
                    missCluster = 1:k;
                    missCluster(ll) = [];
                    missNum = length(missCluster);
                    [dump,idx] = sort(val);
                    label(idx(1:missNum)) = missCluster;
                end
                E = sparse(1:n,label,1,n,k,n);  % transform label into indicator matrix
                center = full((E*spdiags(1./sum(E,1)',0,k,k))'*X);    % compute center of each cluster
                centernorm = sqrt(sum(center.^2, 2));
                center = center ./ centernorm(:,ones(1,p));
                it=it+1;
            end
            if it<maxit
                bCon = true;
            end
            if isempty(bestlabel)
                bestlabel = label;
                bestcenter = center;
                if reps>1
                    if any(label ~= last)
                        W=full(X*center');
                    end
                    D = 1-W;
                    for j = 1:k
                        sumD(j) = sum(D(label==j,j));
                    end
                    bestsumD = sumD;
                    bestD = D;
                end
            else
                if any(label ~= last)
                    W=full(X*center');
                end
                D = 1-W;
                for j = 1:k
                    sumD(j) = sum(D(label==j,j));
                end
                if sum(sumD) < sum(bestsumD)
                    bestlabel = label;
                    bestcenter = center;
                    bestsumD = sumD;
                    bestD = D;
                end
            end
    end
end

label = bestlabel;
center = bestcenter;
if reps>1
    sumD = bestsumD;
    D = bestD;
elseif nargout > 3
    switch distance
        case 'sqeuclidean'
            if it>=maxit
                aa = full(sum(X.*X,2));
                bb = full(sum(center.*center,2));
                ab = full(X*center');
                D = bsxfun(@plus,aa,bb') - 2*ab;
                D(D<0) = 0;
            else
                aa = full(sum(X.*X,2));
                D = aa(:,ones(1,k)) + D;
                D(D<0) = 0;
            end
            D = sqrt(D);
        case 'cosine'
            if it>=maxit
                W=full(X*center');
            end
            D = 1-W;
    end
    for j = 1:k
        sumD(j) = sum(D(label==j,j));
    end
end

function [eid,emsg,varargout]=getargs(pnames,dflts,varargin)
%GETARGS Process parameter name/value pairs
%   [EID,EMSG,A,B,...]=GETARGS(PNAMES,DFLTS,'NAME1',VAL1,'NAME2',VAL2,...)
%   accepts a cell array PNAMES of valid parameter names, a cell array
%   DFLTS of default values for the parameters named in PNAMES, and
%   additional parameter name/value pairs.  Returns parameter values A,B,...
%   in the same order as the names in PNAMES.  Outputs corresponding to
%   entries in PNAMES that are not specified in the name/value pairs are
%   set to the corresponding value from DFLTS.  If nargout is equal to
%   length(PNAMES)+1, then unrecognized name/value pairs are an error.  If
%   nargout is equal to length(PNAMES)+2, then all unrecognized name/value
%   pairs are returned in a single cell array following any other outputs.
%
%   EID and EMSG are empty if the arguments are valid.  If an error occurs,
%   EMSG is the text of an error message and EID is the final component
%   of an error message id.  GETARGS does not actually throw any errors,
%   but rather returns EID and EMSG so that the caller may throw the error.
%   Outputs will be partially processed after an error occurs.
%
%   This utility can be used for processing name/value pair arguments.
%
%   Example:
%       pnames = {'color' 'linestyle', 'linewidth'}
%       dflts  = {    'r'         '_'          '1'}
%       varargin = {{'linew' 2 'nonesuch' [1 2 3] 'linestyle' ':'}
%       [eid,emsg,c,ls,lw] = statgetargs(pnames,dflts,varargin{:})    % error
%       [eid,emsg,c,ls,lw,ur] = statgetargs(pnames,dflts,varargin{:}) % ok

% We always create (nparams+2) outputs:
%    one each for emsg and eid
%    nparams varargs for values corresponding to names in pnames
% If they ask for one more (nargout == nparams+3), it's for unrecognized
% names/values

%   Original Copyright 1993-2008 The MathWorks, Inc.
%   Modified by Deng Cai (dengcai@gmail.com) 2011.11.27

% Initialize some variables
emsg = '';
eid = '';
nparams = length(pnames);
varargout = dflts;
unrecog = {};
nargs = length(varargin);

% Must have name/value pairs
if mod(nargs,2)~=0
    eid = 'WrongNumberArgs';
    emsg = 'Wrong number of arguments.';
else
    % Process name/value pairs
    for j=1:2:nargs
        pname = varargin{j};
        if ~ischar(pname)
            eid = 'BadParamName';
            emsg = 'Parameter name must be text.';
            break;
        end
        i = strcmpi(pname,pnames);
        i = find(i);
        if isempty(i)
            % if they've asked to get back unrecognized names/values, add this
            % one to the list
            if nargout > nparams+2
                unrecog((end+1):(end+2)) = {varargin{j} varargin{j+1}};
                % otherwise, it's an error
            else
                eid = 'BadParamName';
                emsg = sprintf('Invalid parameter name:  %s.',pname);
                break;
            end
        elseif length(i)>1
            eid = 'BadParamName';
            emsg = sprintf('Ambiguous parameter name:  %s.',pname);
            break;
        else
            varargout{i} = varargin{j+1};
        end
    end
end

varargout{nparams+1} = unrecog;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ U, V, W ] = vals2coeffs( X,Y,Z )
% VALS2COEFFS  componentwise conversion of matrices of values to
% matrices of 2D Fourier coefficients.
%
% U, V, W = VALS2COEFFS( X,Y, Z ) converts matrices X, Y and Z of values
% sampled from doubly periodic functions on equally-spaced tensor grids of
% the domain [-pi, pi) x [-pi, pi) to matrices U, V and W,
% containing 2D Fourier coefficients for the corresponding interpolants.
%
%
% See also COEFFS2VALS

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

U = spherefun.vals2coeffs(X);
V = spherefun.vals2coeffs(Y);
W = spherefun.vals2coeffs(Z);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%PLOTF Plot feature distribution, special version
%
%   h = PLOTF(A,N)
%
% Produces 1-D density plots for all the features in dataset A. The
% densities are estimated using PARZENML. N is the number of
% feature density plots on a row.
%
% See also DATASETS, PARZENML

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: plotf.m,v 1.6 2009/11/13 08:54:18 davidt Exp $

function h_out = plotf(a,n,z)
  %DXD make a standard setting for n, I'm getting crazy!
  if nargin<2
    n = 1;
  end

  if ~isdataset(a)
    a = prdataset(a,1); % solves a lot of problems
  end
  [m,k,c] = getsize(a);

	% Define the color for each of the classes:
  if c == 1
    clrmap = [0 0 1];
	elseif c == 2
		clrmap = [0 0 1; 1 0 0];
	else
		clrmap = hsv(c);
	end

	% Make subplots for each feature, so a grid of p x q subplots is
	% defined
	h = [];
	if k >= n
		p = ceil(k/n); q = n;
	else
		p = k; q = 1;
  end

  if isempty(getfeatlab(a))
    a = setfeatlab(a,[1:k]');
  end
  % Get the feature names
  feats = getfeatlab(a,'string');
	%DXD what happens here?!
  %RD If feature labels are scalars to single characters it might be
  %nicer to put 'Feature ' in front of it.
	if size(feats,2) == 1
		feats = [repmat('Feature ',size(feats,1),1) feats];
	end
	if isempty(feats)
		feats = num2str((1:k)');
	end

	% Make the plot for each of the features:
	for j = 1:k
		b = a(:,j);
		s = zeros(1,c);
		d = zeros(121,c);
		bb = [-0.10:0.01:1.10]' * (max(b)-min(b)) + min(b);
		ex = 0;
		% Make a density estimate of each of the classes:
		for i = 1:c
			I = findnlab(a,i);
			D = +distm(bb,b(I,:));
      if nargin < 3
        s(i) = parzenml(b(I,:));
      else
        s(i) = z(i);
      end
			% Compute the density function
			d(:,i) = sum(exp(-D/(s(i).^2)),2)./(length(I)*s(i));;
		end
		% Create the subplots with the correct sizes:
    if p==1 && q==1
      % avoid subplots in case of single plot
      plot(bb,zeros(size(bb)),'w.');
      hold on;
    else
      subplot(p,q,j)
      plot(bb,zeros(size(bb)),'w.');
      hold on
    end
		h = [];
		% Scatter the data and plot the density functions for each of the
		% classes:
		for i = 1:c
			I = findnlab(a,i);
			hh = plot(b(I),zeros(size(b(I))),'x',bb,+d(:,i));
			set(hh,'color',clrmap(i,:));
			h = [h;hh];
		end
		legend(h(1:2:end)',num2str(getlablist(a))); %does not work properly
		title([getname(a) ': ' feats(j,:)]);
		V = axis;
		axis([bb(1) bb(end) V(3) V(4)]);
		set(gca,'xtick',[]);
		set(gca,'ytick',[]);
		xlabel(feats(j,:));
		hold off
	end

	% The last details to take care of:
	if k == 1, title(''); end
	if nargout > 0
		h_out = h;
	end

	return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function varargout=fxpefac(varargin)
%
% For calling details please see v_fxpefac.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_fxpefac.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: fxpefac.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_fxpefac(varargin{:});
else
    v_fxpefac(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function Z=conv(X,Y)
%CONV (overloaded)

if isnumeric(X)
    temp = X;
    X = Y;
    Y = temp;
end
if isa(X,'sdpvar') & isa(Y,'sdpvar')
    error('nonlinear CONV not supported yet. Make feature request');
end
x_lmi_variables = X.lmi_variables;
n = X.dim(1);
m = X.dim(2);
Z=X;
ii = [];
jj = [];
ss = [];
for i = 1:length(x_lmi_variables)+1
    x=reshape(X.basis(:,i),n,m);
    z=conv(full(x),full(Y));
    [iiz,jjz,ssz] = find(z(:));
    ii = [ii iiz(:)'];
    jj = [jj jjz(:)'+i-1];
    ss = [ss ssz(:)'];
end
Z.basis = sparse(ii,jj,ss);
Z.dim(1) = size(z,1);
Z.dim(2) = size(z,2);
Z = clean(Z);
% Reset info about conic terms
if isa(Z,'sdpvar')
    Z.conicinfo = [0 0];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [GAmp,GTime]=GzAreaTrapezoid2(p) % unit m-1
% Trapezoid with prescribed area & maximum slew rate & maximum gradient amp
% note: don't support ramp sampling

global VCtl
global VObj

tStart=p.tStart; % start time
Area=abs(p.Area);   % trapezoid area m-1
Duplicates=max(1,p.Duplicates);
DupSpacing=max(0,p.DupSpacing);

Area=Area/(VObj.Gyro/(2*pi));
GzAmp=VCtl.MaxGrad;   % Gz amplitude
if Area/GzAmp < VCtl.MinUpdRate
    tEnd=VCtl.MinUpdRate + tStart;
    GzAmp=Area/VCtl.MinUpdRate;   % Gz amplitude
else
    tEnd=Area/GzAmp + tStart;
end

tRamp=max(VCtl.MinUpdRate,GzAmp/VCtl.MaxSlewRate);   % ramp time

[GAmp,GTime]=StdTrap(tStart-tRamp, ...
                     tEnd+tRamp,   ...
                     tStart,               ...
                     tEnd,                 ...
                     GzAmp*sign(p.Area),2,2,2);

[GTime,m,n]=unique(GTime);
GAmp=GAmp(m);

% Create Duplicates
if Duplicates~=1 & DupSpacing ~=0
    GAmp=repmat(GAmp,[1 Duplicates]);
    TimeOffset = repmat(0:DupSpacing:(Duplicates-1)*DupSpacing,[length(GTime) 1]);
    GTime=repmat(GTime,[1 Duplicates]) + (TimeOffset(:))';
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [rho,rhou,rhov,Ener] = ChannelIC2D(x, y, time);

%  function [Q] = ChannelIC2D(x, y, time)
%  Purpose: Impose uniform plane flow

mu = 1e-2; pbar = 10; gamma = 1.5;

rho  = 1;
rhou = y.^2;
rhov = 0;
Ener = (2*mu*x + pbar)/(gamma-1) + .5*(y.^4);
return

repo: Aero-matlab


text:
function gamma=ref_gabdualns_3(g,V);
%REF_GABDUALNS_3  GABDUALNS by multiwindow method.
%   Usage:  gamma=ref_gabdualns_3(g,V);
%

[gm,a,M]=ref_nonsep2multiwin(g,V);

gmd=gabdual(gm,a,M);

gamma=gmd(:,1);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% RES = reconLpyr(PYR, INDICES, LEVS, FILT2, EDGES)
%
% Reconstruct image from Laplacian pyramid, as created by buildLpyr.
%
% PYR is a vector containing the N pyramid subbands, ordered from fine
% to coarse.  INDICES is an Nx2 matrix containing the sizes of
% each subband.  This is compatible with the MatLab Wavelet toolbox.
%
% LEVS (optional) should be a list of levels to include, or the string
% 'all' (default).  The finest scale is number 1.  The lowpass band
% corresponds to lpyrHt(INDICES)+1.
%
% FILT2 (optional) can be a string naming a standard filter (see
% namedFilter), or a vector which will be used for (separable)
% convolution.  Default = 'binom5'.  EDGES specifies edge-handling,
% and defaults to 'reflect1' (see corrDn).

% Eero Simoncelli, 6/96

function res = reconLpyr(pyr, ind, levs, filt2, edges)

if (nargin < 2)
  error('First two arguments (PYR, INDICES) are required');
end

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

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

if (exist('filt2') ~= 1)
  filt2 = 'binom5';
end

if (exist('edges') ~= 1)
  edges= 'reflect1';
end
%%------------------------------------------------------------

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

if isstr(filt2)
  filt2 = namedFilter(filt2);
end

filt2 = filt2(:);
res_sz = ind(1,:);

if any(levs > 1)

  int_sz = [ind(1,1), ind(2,2)];

  nres = reconLpyr( pyr(prod(res_sz)+1:size(pyr,1)), ...
      ind(2:size(ind,1),:), levs-1, filt2, edges);

  if (res_sz(1) == 1)
    res = upConv(nres, filt2', edges, [1 2], [1 1], res_sz);
  elseif (res_sz(2) == 1)
    res = upConv(nres, filt2, edges, [2 1], [1 1], res_sz);
  else
    hi = upConv(nres, filt2, edges, [2 1], [1 1], int_sz);
    res = upConv(hi, filt2', edges, [1 2], [1 1], res_sz);
  end

else

  res = zeros(res_sz);

end

if any(levs == 1)
  res = res + pyrBand(pyr,ind,1);
end

repo: SelfExSR


text:
function [min] = hr2min(hr)
% Convert time from hours to minutes.
% Chad A. Greene 2012
min = hr*60;

author: Sable

repo: mcbench-benchmarks


text:
function pop = crossoverOp(opt, pop, state)
% Function: pop = crossoverOp(opt, pop, state)
% Description: Crossover operator. All of the individuals would be do crossover, but
%   only "crossoverFraction" of design variables of an individual would changed.
%
%         LSSSSWC, NWPU
%    Revision: 1.1  Data: 2011-07-13
%*************************************************************************

%*************************************************************************
% 1. Check for the parameters
%*************************************************************************
% determine the crossover method
strfun = lower(opt.crossover{1});
numOptions = length(opt.crossover) - 1;
[crossoverOpt{1:numOptions}] = opt.crossover{2:end};

switch( strfun )
    case 'intermediate'
        fun = @crsIntermediate;
    otherwise
        error('NSGA2:CrossoverOpError', 'No support crossover operator!');
end

nVar = opt.numVar;

% "auto" crossover fraction
if( ischar(opt.crossoverFraction) )
    if( strcmpi(opt.crossoverFraction, 'auto') )
        fraction = 2.0 / nVar;
    else
        error('NSGA2:CrossoverOpError', 'The "crossoverFraction" parameter should be scalar or "auto" string.');
    end
else
    fraction = opt.crossoverFraction;
end

for ind = 1:2:length(pop)    % Popsize should be even number
        % Create children
        [child1, child2] = fun( pop(ind), pop(ind+1), fraction, crossoverOpt );

        % Round
        for v = 1:nVar
            if( opt.vartype(v) == 2)
                child1.var(v) = round( child1.var(v) );
                child2.var(v) = round( child2.var(v) );
            end
        end

        % Bounding limit
        child1.var = varlimit(child1.var, opt.lb, opt.ub);
        child2.var = varlimit(child2.var, opt.lb, opt.ub);

        pop(ind)     = child1;
        pop(ind+1)   = child2;

end

function [child1, child2] = crsIntermediate(parent1, parent2, fraction, options)
% Function: [child1, child2] = crsIntermediate(parent1, parent2, fraction, options)
% Description: (For real coding) Intermediate crossover. (Same as Matlab's crossover
%   operator)
%       child = parent1 + rand * Ratio * ( parent2 - parent1)
% Parameters:
%   fraction : crossover fraction of variables of an individual
%   options = ratio
%
%         LSSSSWC, NWPU
%    Revision: 1.1  Data: 2011-07-13
%*************************************************************************

if( length(options)~=1 || ~isnumeric(options{1}))
    error('NSGA2:CrossoverOpError', 'Crossover operator parameter error!');
end

ratio = options{1};

child1 = parent1;
child2 = parent2;

nVar = length(parent1.var);
crsFlag = rand(1, nVar) < fraction;

randNum = rand(1,nVar);     % uniformly distribution

child1.var = parent1.var + crsFlag .* randNum .* ratio .* (parent2.var - parent1.var);
child2.var = parent2.var - crsFlag .* randNum .* ratio .* (parent2.var - parent1.var);

author: Eric-Bradford

repo: TS-EMO


text:
function neqn = p22_neqn ( )

%*****************************************************************************80
%
%% P22_NEQN returns the number of equations for problem P22.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NEQN, the number of equations.
%
  neqn = 2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [gal] = m32gal(m3)
% Convert volume from cubic meters to US liquid gallons.
% Chad Greene 2012
gal = m3*264.17205236;


text:
function [dat, dimord] = read_shm_data(hdr, chanindx, begtrial, endtrial)

% READ_SHM_DATA reads the data in real-time from shared memory
% this is a helper function for READ_DATA

% Copyright (C) 2007-2009, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% this persistent variable is used for caching the 600 packets
% which inproves the throughput when reading overlapping data segments
persistent ctf_shm

% read the data from shared memory, first the meta information only
[msgType msgId sampleNumber numSamples numChannels] = read_ctf_shm;

if isempty(ctf_shm)
  ctf_shm.msgType      = nan(size(msgType));
  ctf_shm.msgId        = nan(size(msgId));
  ctf_shm.sampleNumber = nan(size(sampleNumber));
  ctf_shm.numSamples   = nan(size(numSamples));
  ctf_shm.numChannels  = nan(size(numChannels));
  ctf_shm.data = {};
  ctf_shm.hit   = 0;
  ctf_shm.mis   = 0;
end

% there seems to be a bug in Acq, causing the messageId to wrap around
% hence it cannot be used as index into the packets, so construct a new trial numbering vector
trlNum = nan(size(msgId));
trlNum(msgType==1) = sampleNumber(msgType==1)./numSamples(msgType==1) + 1;

% allocate memory for the data, fill with NaNs
dat = nan(length(chanindx), hdr.nSamples, endtrial-begtrial+1);

% determine which trials/packets to read
sel = find((trlNum>=begtrial) & (trlNum<=endtrial));

% this is for calibrating the integer values to physical values
if isfield(hdr, 'orig') && isfield(hdr.orig, 'gainV')
  gain = diag(hdr.orig.gainV(chanindx));
elseif isfield(hdr, 'orig') && isfield(hdr.orig, 'res4')
  gain = diag([hdr.orig.res4.senres(chanindx).gain]);
end

if length(chanindx)>1
  % this speeds up the multiplication, but would result in a sparse matrix when nchans=1
  gain = sparse(gain);
end

for i=1:length(sel)
  % nchan = numChannels(i);
  % nsmp  = numSamples(i);
  nchan = hdr.nChans;
  nsmp  = hdr.nSamples;
  % buf = read_ctf_shm(sel(i));
  % buf = buf(1:nchan*nsmp);
  if all([msgType(sel(i)) msgId(sel(i)) sampleNumber(sel(i)) numSamples(sel(i)) numChannels(sel(i))] ==  [ctf_shm.msgType(sel(i)) ctf_shm.msgId(sel(i)) ctf_shm.sampleNumber(sel(i)) ctf_shm.numSamples(sel(i)) ctf_shm.numChannels(sel(i))])
    % get the data from the cache
    buf = ctf_shm.data{sel(i)};
    ctf_shm.hit = ctf_shm.hit+1;
  else
    % read the data from shared memory, update the cache
    buf = read_ctf_shm(sel(i),nchan*nsmp);
    ctf_shm.data{sel(i)}         = buf;
    ctf_shm.msgType(sel(i))      = msgType(sel(i));
    ctf_shm.msgId(sel(i))        = msgId(sel(i));
    ctf_shm.sampleNumber(sel(i)) = sampleNumber(sel(i));
    ctf_shm.numSamples(sel(i))   = numSamples(sel(i));
    ctf_shm.numChannels(sel(i))  = numChannels(sel(i));
    ctf_shm.mis                  = ctf_shm.mis+1;
  end
  buf = reshape(buf, nchan, nsmp);
  thistrial = trlNum(sel(i));
  % apply calibration to the selected channels
  dat(:,:,thistrial-begtrial+1) = gain * double(buf(chanindx,:));
end

% if any(isnan(dat(:)))
%   ft_warning('data has been padded with NaNs');
%   fprintf('trials present   = %d - %d\n', min(trlNum), max(trlNum));
%   fprintf('trials requested = %d - %d\n', begtrial, endtrial);
% end

dimord = 'chans_samples_trials';

author: spm


text:
function idx = argmax(x)

    [~,idx] = max(x);

author: gordonberman

repo: MotionMapper


text:
function vAberr=aberCorr(vOrig,obsVel,sunDist)
%%ABERCORR Rotate vectors for stellar aberration, approximately including
%          the effects of the gravitational potential of the Sun if the
%          distances between the observers and the Sun are provided.
%
%INPUTS: vOrig A 3XN set of N vectors pointing from the observers to the
%              objects being observed without corruption due to
%              aberration. The units of the vector do not matter.
%              Normally, one would expect vOrig to be composed of unit
%              vectors as aberration affects the pointing direction of the
%              vectors.
%       obsVel The 3XN matrix of the velocity of each observer with
%              respect to the stellar coordinate system origin in meters
%              per second. For example, if using the barycentric celestial
%              reference system (BCRS), they are the set of velocity
%              vectors of the observers with respect to the barycenter.
%      sunDist An optional 3XN vector of the scalar distances in meters
%              from the Sun to the observer. If this parameter is
%              included, an approximate correction due to the
%              gravitational potential of the Sun is included.
%
%OUTPUTS: vAberr The 3XN set of N vectors rotated to deal with aberration
%                effects.
%
%If all parameters are provided, the function is essentially a wrapper for
%the function iauAb in the International Astronomical Union's Standards of
%Fundamental Astronomy library with an adjustment so that the vector in
%question need not have unit magnitude. If the distance to the Sun is
%omitted, then the standard special relativistic aberration correction
%without any gravitational effects is applied. The standard special
%relativistic aberration correction is described in Chapter 7 of [1]
%
%Note that if the vectors provided are meant to be apparent distances,
%rather than just unit vectors representing directions, the transformation
%does not adjust the magnitudes of the vectors to account for special
%relativistic contraction of space due to the motion of the observer with
%respect to the origin. The magnitudes of the output vectors equal the
%magnitudes of the input vectors.
%
%The algorithm can be compiled for use in Matlab  using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%vAberr=aberCorr(vOrig,obsVel,sunDist);
%or
%vAberr=aberCorr(vOrig,obsVel);
%
%REFERENCES:
%[1] S. E. Urban and K. P.Seidelmann, Eds.,Explanatory Supplement to the
%    Astronomical Almanac, 3rd ed. Mill Valley, CA: University Science
%    Books, 2013.
%
%March 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function circle_integrals_test01 ( )

%*****************************************************************************80
%
%% CIRCLE_INTEGRALS_TEST01 tests CIRCLE01_SAMPLE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 January 2014
%
%  Author:
%
%    John Burkardt
%
  m = 2;
  n = 4192;
  test_num = 20;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  Use CIRCLE01_SAMPLE to compare exact and\n' );
  fprintf ( 1, '  estimated integrals along the circumference \n' );
  fprintf ( 1, '  of the unit circle in 2D.\n' );
%
%  Get sample points.
%
  seed = 123456789;
  [ x, seed ] = circle01_sample ( n, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of sample points used is %d\n', n );
%
%  Randomly choose X, Y exponents.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  If any exponent is odd, the integral is zero.\n' );
  fprintf ( 1, '  We restrict this test to randomly chosen even exponents.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Ex  Ey     MC-Estimate           Exact      Error\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    [ e, seed ] = i4vec_uniform_ab ( m, 0, 5, seed );

    e(1:m) = e(1:m) * 2;

    value = monomial_value ( m, n, e, x );

    result = circle01_length ( ) * sum ( value(1:n) ) / n;
    exact = circle01_monomial_integral ( e );
    error = abs ( result - exact );

    fprintf ( 1, '  %2d  %2d  %14.6g  %14.6g  %10.2e\n', ...
      e(1:m), result, exact, error );

  end

  return
end


text:
function [ n_data, x, fx ] = gamma_log_values ( n_data )

%*****************************************************************************80
%
%% GAMMA_LOG_VALUES returns some values of the Log Gamma function.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      Log[Gamma[x]]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 20;

  fx_vec = [ ...
      0.1524063822430784E+01, ...
      0.7966778177017837E+00, ...
      0.3982338580692348E+00, ...
      0.1520596783998375E+00, ...
      0.0000000000000000E+00, ...
     -0.4987244125983972E-01, ...
     -0.8537409000331584E-01, ...
     -0.1081748095078604E+00, ...
     -0.1196129141723712E+00, ...
     -0.1207822376352452E+00, ...
     -0.1125917656967557E+00, ...
     -0.9580769740706586E-01, ...
     -0.7108387291437216E-01, ...
     -0.3898427592308333E-01, ...
     0.00000000000000000E+00, ...
     0.69314718055994530E+00, ...
     0.17917594692280550E+01, ...
     0.12801827480081469E+02, ...
     0.39339884187199494E+02, ...
     0.71257038967168009E+02 ];

  x_vec = [ ...
      0.20E+00, ...
      0.40E+00, ...
      0.60E+00, ...
      0.80E+00, ...
      1.00E+00, ...
      1.10E+00, ...
      1.20E+00, ...
      1.30E+00, ...
      1.40E+00, ...
      1.50E+00, ...
      1.60E+00, ...
      1.70E+00, ...
      1.80E+00, ...
      1.90E+00, ...
      2.00E+00, ...
      3.00E+00, ...
      4.00E+00, ...
     10.00E+00, ...
     20.00E+00, ...
     30.00E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function [comod] = bz_CrossFreqMod(lfp,phaserange,amprange,varargin)
% [comod] = bz_CrossFreqMod(lfp,phaserange,amprange,flagPlot)
%
%
%This function calculates the modulation index of phase-amplitude between
%phaserange (lower frequencies) to amplitude range (higher frequencies).
%It can really take a long time if you do very small steps of frequency,
%due to wavelet processing each frequency at a time.
%
%%INPUT
%    lfp            a buzcode structure with fields lfp.data,
%                                                   lfp.timestamps
%                                                   lfp.samplingRate
%                   -lfp can also be a [t x 1] timeseries signal. in which
%                   case you need to input 'samplingRate'
%    phaserange     [min:steps:max] array of frequencies to filter for the
%                   phase signal
%    amprange       [min:stepsmax] array of frequencies range for wavelets
%                   for the power signal
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     phaseCh      channel to compute phase. If empty take first channel
%     ampChans     channels to compute amplitude. If empty take first channel
%     makePlot      default true
%    =========================================================================
%
%OUTPUT
%   comod               phase-frequency versus amplitude-frequency
%                       comodulogram matrix in modulation index units
%
%Dependencies
%   bz_Filter
%   bz_WaveSpec
%
%   AntonioFR, 2/2019

%% Parse inputs

p = inputParser;
addParameter(p,'phaseCh',lfp.channels(1),@isnumeric);
addParameter(p,'ampChans',lfp.channels(1),@isnumeric);
addParameter(p,'samplingRate',1250,@isnumeric);
addParameter(p,'makePlot',true,@islogical);

parse(p,varargin{:});
phaseCh = p.Results.phaseCh;
ampChans = p.Results.ampChans;
samplingRate = p.Results.samplingRate;
makePlot = p.Results.makePlot;

%lfp input
if isstruct(lfp)
    data = lfp.data;
    timestamps = lfp.timestamps;
    samplingRate = lfp.samplingRate;
elseif iscell(lfp) %for multiple trials
    celllengths = cellfun(@length,lfp);
    data = vertcat(lfp{:});
elseif isnumeric(lfp)
    data = lfp;
    timestamps = [1:length(lfp)]'./samplingRate;
end

%% Filter LFP for phase
nfreqs = length(amprange);

for bnd = 1:length(phaserange)-1
    filtered_phase(bnd,:) = bz_Filter(lfp,'passband',phaserange(bnd:bnd+1),'filter','fir1','channels',phaseCh);
end

%% Wavelet Transform LFP in intervals
for ch = 1:length(ampChans)
    comod = zeros(length(amprange)-1,length(filtered_phase),length(ampChans));

    lfpCh = lfp; % this should be remove when bz_WaveSpec can take 'channels' input
    lfpCh.data = lfp.data(:,find(lfp.channels == ampChans(ch)));

    for apr = 1:length(amprange)-1
        wavespec_amp = bz_WaveSpec(lfp,'frange',[amprange(apr) amprange(apr+1)],'nfreqs',1);

        wavespec_amp.data = abs(wavespec_amp.data);
        %% Bin phase and power
        numbins = 50;
        phasebins = linspace(-pi,pi,numbins+1);
        phasecenters = phasebins(1:end-1)+(phasebins(2)-phasebins(1));

        for idx = 1:length(filtered_phase)
            [phasedist,~,phaseall] = histcounts(filtered_phase(idx).phase,phasebins);

            phaseAmp = zeros(numbins,1);
            for bb = 1:numbins
                phaseAmp(bb) = mean(wavespec_amp.data(phaseall==bb),1);
            end

            phaseAmp = phaseAmp./sum(phaseAmp,1);
            comod(apr,idx,ch) = sum(phaseAmp.*log(phaseAmp./(ones(numbins,size(phaseAmp,2))/numbins)))/log(numbins);
        end

    end
    ampfreqs = wavespec_amp.freqs;
    clear lfpCh
end

%% Plot
if makePlot

    figure
    for ch = 1:length(ampChans)
        subplot(1,length(ampChans),ch);
        contourf(phaserange(2:end),amprange(2:end),abs(comod(:,:,ch)),20,'LineColor','none');
        colorbar ('SouthOutside'); colormap jet;
        %title(signals)
        if ch > 1
            set(gca,'YTick',[]);
        end
    end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [X,n] = shiftdim(varargin)
% SHIFTDIM (overloaded)

Y = varargin{1};
X = Y;
X.basis = [];
for i = 1:size(Y.basis,2)
    base = reshape(full(Y.basis(:,i)),[X.dim(1) X.dim(2)]);
    [base,n] = shiftdim(base,varargin{2:end});
    X.basis = [X.basis sparse(base(:))];
end
[X.dim(1), X.dim(2)] = size(base);
X.conicinfo = [0 0];
X = clean(X);


text:
function vos_map = obtain_vos_map(img_q,N)
    % build color cube (color histogram)
    [row,col] = size(img_q(:,:,1));
    color_cube = zeros(N,N,N);
    for i = 1 : row
        for j = 1 : col
            cr = img_q(i,j,1);
            cg = img_q(i,j,2);
            cb = img_q(i,j,3);
            color_cube(cr,cg,cb) = color_cube(cr,cg,cb) + 1;
        end
    end

    % label major color
    num_counted_pixels = 0; % record numder of pixels that have been counted
    num_major_color = 0; % record numder of major colors
    prop = 0.95; % care about 95%
    while num_counted_pixels < prop*row*col
        color_most_n = max(max(max(color_cube))); % number of the most color at this iteration
        same_n = length(color_cube(color_cube == color_most_n));
        num_counted_pixels = num_counted_pixels + same_n*color_most_n;
        num_major_color = num_major_color + same_n;
        % label the counted points in color_cube as -1, meaning that these
        % colors are "major"
        color_cube(color_cube == color_most_n) = -1;
    end

    % assign the rest 5% colors to their near major color
    [r_set,g_set,b_set] = ind2sub(size(color_cube),find(color_cube==-1));
    while ~isempty(color_cube(color_cube>0))
        color_most_n = max(max(max(color_cube)));
        [cr,cg,cb] = ind2sub(size(color_cube),find(color_cube==color_most_n));
        color_cube(color_cube==color_most_n) = -2;
        for i = 1 : length(cr)
            dist = zeros(length(r_set),1);
            for j = 1 : length(r_set)
                dist(j) = (cr(i)-r_set(j))^2 + (cg(i)-g_set(j))^2 + (cb(i)-b_set(j))^2;
            end
            label = find(dist == min(dist));
            r_assign = r_set(label(1));
            g_assign = g_set(label(1));
            b_assign = b_set(label(1));
            for m = 1 : row
                for n = 1 : col
                    if img_q(m,n,1) == cr(i) && img_q(m,n,2)==cg(i) && img_q(m,n,3)==cb(i)
                        img_q(m,n,:) = reshape([r_assign,g_assign,b_assign],1,1,3);
                    end
                end
            end
        end
    end
    % compute vos saliency
    [r_set,g_set,b_set] = ind2sub(size(color_cube),find(color_cube==-1));
    sal = zeros(length(r_set),1);
    vos_map = zeros(row,col);
    for i = 1 : length(r_set)
        for j = 1 : length(r_set)
            if j == i
                continue
            end
            sal(i) = sal(i) + sqrt((r_set(i)-r_set(j))^2+(g_set(i)-g_set(j))^2+(b_set(i)-b_set(j))^2);
        end
        for m = 1 : row
            for n = 1 : col
                if img_q(m,n,1)==r_set(i) && img_q(m,n,2)==g_set(i) && img_q(m,n,3)==b_set(i)
                    vos_map(m,n) = sal(i);
                end
            end
        end
    end
    vos_map = mat2gray(vos_map);
end


text:
function [J, grad] = cofiCostFunc(params, Y, R, num_users, num_movies, ...
                                  num_features, lambda)
    %COFICOSTFUNC Collaborative filtering cost function
    %   [J, grad] = COFICOSTFUNC(params, Y, R, num_users, num_movies, ...
    %   num_features, lambda) returns the cost and gradient for the
    %   collaborative filtering problem.
    %

    % Unfold the U and W matrices from params
    X = reshape(params(1:num_movies*num_features), num_movies, num_features);
    Theta = reshape(params(num_movies*num_features+1:end), ...
                    num_users, num_features);

    % You need to return the following values correctly
    J = 0;
    X_Gradient = zeros(size(X));
    Theta_Gradient = zeros(size(Theta));

    % ====================== YOUR CODE HERE ======================
    % Instructions: Compute the cost function and gradient for collaborative
    %               filtering. Concretely, you should first implement the cost
    %               function (without regularization) and make sure it is
    %               matches our costs. After that, you should implement the
    %               gradient and use the checkCostFunction routine to check
    %               that the gradient is correct. Finally, you should implement
    %               regularization.
    %
    % Notes: X - num_movies  x num_features matrix of movie features
    %        Theta - num_users  x num_features matrix of user features
    %        Y - num_movies x num_users matrix of user ratings of movies
    %        R - num_movies x num_users matrix, where R(i, j) = 1 if the
    %            i-th movie was rated by the j-th user
    %
    % You should set the following variables correctly:
    %
    %        X_grad - num_movies x num_features matrix, containing the
    %                 partial derivatives w.r.t. to each element of X
    %        Theta_grad - num_users x num_features matrix, containing the
    %                     partial derivatives w.r.t. to each element of Theta
    %

    J_Temporary = power ((X * Theta' - Y), 2);
    J = sum (sum (J_Temporary (R == 1))) / 2 + ((lambda / 2) .* sum (sum (power (Theta, 2)))) + ((lambda / 2) .* sum (sum (power (X, 2))));

    X_Gradient = ((X * Theta' - Y) .* R) * Theta + lambda .* X;
    Theta_Gradient = ((X * Theta' - Y) .* R)' * X + lambda .* Theta;

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

    grad = [X_Gradient(:); Theta_Gradient(:)];

end

author: UtkarshPathrabe


text:
% This file load eq catalogs in hypoellipes format. Since the formated readin-
% is fast, this is the best way to load data for
% large catalogs.
% Earthquakes within a rectengular box (latmin latmax lonmin lonmax)
% a time window and a magnitude threshold can be selected.
% Note: lat/Lon are assumed to be in minutes in the original catalog
%
% last modified: March 96   stefan wiemer

report_this_filefun(mfilename('fullpath'));

%  close the file in case it is still open
safe_fclose(fid);

% selection criteria
tmin = 1960.0001;
tmax = 2005;
lonmin = -180.0;
lonmax =  180.0;
latmin = -90.0;
latmax = 90.0 ;
Mmin = 0.0
deperr = 1000.5;

a = []; b = [];
n2 = 0;

% open the file and read 2000 lines at a time
[file1,path1] = uigetfile([ '*.sum'],' Earthquake Datafile');
fid = fopen([path1 file1],'r') ;;

while  ferror(fid) == ''
    %  l =fscanf(fid,'%2d%2d%2d%2d%2d%4d%2d%*1c%4d%3d%*1c%4d%5f%2d%*31c%2d%*15c%4d',[14 2000]) ;
    % variabl #   1  2  3  4  5  6  7     8  9      10 11 12    13    14
    % position    2  4  6  8 10 14 16  17 21 24 25  29 34 36 45  49 54  58
    %  l =fscanf(fid,'%2d%2d%2d%2d%2d%4d%2d%*1c%4d%3d%*1c%4d%5f%2d%*9c%4d%*5c%4d',[14 2000]) ;
    %            yr mo da hr mi
    % variabl #   1  2  3  4  5     6  7     8   9     10 11      12     13 14     15
    % position    2  4  6  8 10 11 17  20 21 26 30  31 36 43   81 84  85 87 91  93 97
    l =fscanf(fid,'%4d%2d%2d%2d%2d%*1c%6d%3d%*1c%5d%4d%*1c%5d%7d%*38%3d%*1c%2d%4d%*2c%4d',[15   2000]);

    n2 = n2 +  length(l(1,:));

    % make the catalog in the correct format
    % lon lat yr mo da mag dep hr min dip_direction dip rake misfit
    b = [ -l(9,:)-l(10,:)/6000 ; l(7,:)+l(8,:)/6000 ;...
        l(1,:);l(2,:);l(3,:);l(12,:)/10;l(11,:)/100;l(4,:);...
        l(5,:);l(13,:); l(14,:);l(15,:);l(16,:);l(17,:);l(18,:)];
    b = b';
    % l =  b(:,6) >= Mmin & b(:,1) >= lonmin & b(:,1) <= lonmax & ...
    % b(:,2) >= latmin & b(:,2) <= latmax & b(:,3) <= tmax  & ...
    % b(:,3) >= tmin & b(:,10) <= deperr;
    a = [a ; b];
    disp([ num2str(n2,10) ' earthquakes scanned, ' num2str(length(a)) ' EQ found'])
    %  if max(b(:,3)) >  tmax ; break; end
end
ferror(fid)
fclose(fid);
dep1 = 0.3*(max(a.Depth)-min(a.Depth))+min(a.Depth);
dep2 = 0.6*(max(a.Depth)-min(a.Depth))+min(a.Depth);
dep3 = max(a.Depth);

stri1 = [file1];
tim1 = minti;
tim2 = maxti;
minma2 = minma;
maxma2 = maxma;
minde = min(a.Depth);
maxde = max(a.Depth);
rad = 50.;
ic = 0;
ya0 = 0.;
xa0 = 0.;
iwl3 = 1.;
step = 3;
t1p(1) = 80.;
t2p(1) = 85.;
t3p(1) = 90.;
t4p(1) = 93.;
tresh = 10;

% save the catalog
save newcat.mat a

author: CelsoReyes


text:
function jac = jacmodhs52(p)
  m=5;

  jac(1, 1:m)=[4.0, -1.0, 0.0, 0.0, 0.0];
  jac(2, 1:m)=[0.0, 1.0, 1.0, 0.0, 0.0];
  jac(3, 1:m)=[0.0, 0.0, 0.0, 1.0, 0.0];
  jac(4, 1:m)=[0.0, 0.0, 0.0, 0.0, 1.0];

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [xDeformV,yDeformV,zDeformV] = getDeformationAt(deformS,toPlanC,toScanNum,fromPlanC,fromScanNum,fromXYZv)
% function [xDeformV,yDeformV,zDeformV] = getDeformationAt(deformS,toPlanC,toScanNum,fromPlanC,fromScanNum,fromXYZv)
%
% APA, 08/20/2012

% % Obtain base and moving scan UIDs
% baseScanUID = deformS.baseScanUID;
% movScanUID  = deformS.movScanUID;
%
% Figure out whether an inverse Vector field is required
indexToS = toPlanC{end};
indexFromS = fromPlanC{end};
%
% toScanUID = toPlanC{indexToS.scan}(toScanNum).scanUID;
% fromScanUID = fromPlanC{indexFromS.scan}(fromScanNum).scanUID;
% if baseScanUID == fromScanUID
%     calc_inv_vf_flag = 1;
% else
%     calc_inv_vf_flag = 0;
% end
%
% % Create b-spline coefficients file
% bspFileName = fullfile(getCERRPath,'ImageRegistration','tmpFiles',['bsp_coeffs_',baseScanUID,'_',toScanUID,'.txt']);
% success     = write_bspline_coeff_file(bspFileName,deformS.algorithmParamsS);
%
% % Obtain Vf from b-splice coefficients
% vfFileName = fullfile(getCERRPath,'ImageRegistration','tmpFiles',['vf_',baseScanUID,'_',toScanUID,'.mha']);
% system(['plastimatch xf-convert --input ',escapeSlashes(bspFileName), ' --output ', escapeSlashes(vfFileName), ' --output-type vf'])
%
% % Obtain Inverse Vf if required
% if calc_inv_vf_flag
%     system(['plastimatch vf_invert --input ', escapeSlashes(vfFileName), ' --output ', escapeSlashes(vfFileName), ' --dims="x y z" --origin="x y z" --spacing="x y z"'])
% end
%
% % Read .mha file for Vf
% infoS  = mha_read_header(vfFileName);
% data3M = mha_read_volume(infoS);
% xDeform3M = flipdim(permute(data3M(:,:,:,1),[2,1,3]),3);
% yDeform3M = flipdim(permute(data3M(:,:,:,2),[2,1,3]),3);
% zDeform3M = flipdim(permute(data3M(:,:,:,3),[2,1,3]),3);

vf = get_vector_field(deformS,toPlanC,toScanNum,fromPlanC,fromScanNum);

% Get grid coordinates for "fromScan"
[xV, yV, zV] = getScanXYZVals(fromPlanC{indexFromS.scan}(fromScanNum));

% Interolate fromXYZv points to get Vf based on the "fromScan" grid
xDeformV = finterp3(fromXYZv(:,1), fromXYZv(:,2), fromXYZv(:,3), vf(:,:,:,1), [xV(1)-eps*10^10 xV(2)-xV(1) xV(end)+eps*10^10], [yV(end)-eps*10^10 yV(1)-yV(2) yV(1)+eps*10^10], zV, 0);
yDeformV = finterp3(fromXYZv(:,1), fromXYZv(:,2), fromXYZv(:,3), vf(:,:,:,2), [xV(1)-eps*10^10 xV(2)-xV(1) xV(end)+eps*10^10], [yV(end)-eps*10^10 yV(1)-yV(2) yV(1)+eps*10^10], zV, 0);
zDeformV = finterp3(fromXYZv(:,1), fromXYZv(:,2), fromXYZv(:,3), vf(:,:,:,3), [xV(1)-eps*10^10 xV(2)-xV(1) xV(end)+eps*10^10], [yV(end)-eps*10^10 yV(1)-yV(2) yV(1)+eps*10^10], zV, 0);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [g, gdata, gprior] = rbfgrad(net, x, t)
%RBFGRAD Evaluate gradient of error function for RBF network.
%
%	Description
%	G = RBFGRAD(NET, X, T) takes a network data structure NET together
%	with a matrix X of input vectors and a matrix T of target vectors,
%	and evaluates the gradient G of the error function with respect to
%	the network weights (i.e. including the hidden unit parameters). The
%	error function is sum of squares. Each row of X corresponds to one
%	input vector and each row of T contains the corresponding target
%	vector. If the output function is 'NEUROSCALE' then the gradient is
%	only computed for the output layer weights and biases.
%
%	[G, GDATA, GPRIOR] = RBFGRAD(NET, X, T) also returns separately  the
%	data and prior contributions to the gradient. In the case of multiple
%	groups in the prior, GPRIOR is a matrix with a row for each group and
%	a column for each weight parameter.
%
%	See also
%	RBF, RBFFWD, RBFERR, RBFPAK, RBFUNPAK, RBFBKP
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Check arguments for consistency
switch net.outfn
case 'linear'
   errstring = consist(net, 'rbf', x, t);
case 'neuroscale'
   errstring = consist(net, 'rbf', x);
otherwise
   error(['Unknown output function ', net.outfn]);
end
if ~isempty(errstring);
  error(errstring);
end

ndata = size(x, 1);

[y, z, n2] = rbffwd(net, x);

switch net.outfn
case 'linear'

   % Sum squared error at output units
   delout = y - t;

   gdata = rbfbkp(net, x, z, n2, delout);
   [g, gdata, gprior] = gbayes(net, gdata);

case 'neuroscale'
   % Compute the error gradient with respect to outputs
   y_dist = sqrt(dist2(y, y));
   D = (t - y_dist)./(y_dist+diag(ones(ndata, 1)));
   temp = y';
   gradient = 2.*sum(kron(D, ones(1, net.nout)) .* ...
      (repmat(y, 1, ndata) - repmat((temp(:))', ndata, 1)), 1);
   gradient = (reshape(gradient, net.nout, ndata))';
   % Compute the error gradient
   gdata = rbfbkp(net, x, z, n2, gradient);
   [g, gdata, gprior] = gbayes(net, gdata);
otherwise
   error(['Unknown output function ', net.outfn]);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%MDL_COBRA600 Create model of Adept Cobra 600 manipulator
%
% MDL_COBRA600 is a script that creates the workspace variable c600 which
% describes the kinematic characteristics of the 4-axis Adept Cobra 600
% SCARA manipulator using standard DH conventions.
%
% Also define the workspace vectors:
%   qz         zero joint angle configuration
%
% Notes::
% - SI units are used.
%
% See also SerialRevolute, mdl_puma560akb, mdl_stanford.

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

% MODEL: Adept, Cobra600, 4DOF, standard_DH

% hardstop limits included
links = [
    Revolute('d', 0.387, 'a', 0.325, 'qlim', [-50 50]*pi/180);
    Revolute('a', 0.275, 'alpha', pi, 'qlim', [-88 88]*pi/180);
    Prismatic('qlim', [0 0.210]);
    Revolute()
    ];

c600 = SerialLink(links, 'name', 'Cobra600', 'manufacturer', 'Adept', ...
    'plotopt', {'workspace', [0 0.8 -0.6 0.6 0 0.4]} );

qz = [0 0 0 0];

author: petercorke

repo: robotics-toolbox-matlab


text:
function sphere_cubed_grid_lines_test ( )

%*****************************************************************************80
%
%% SPHERE_CUBED_GRID_LINES_TEST tests SPHERE_CUBED_GRID_LINES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 May 2015
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_CUBED_GRID_LINES_TEST\n' );
  fprintf ( 1, '  SPHERE_CUBED_GRID_LINES defines the lines\n' );
  fprintf ( 1, '  on a cubed sphere grid.\n' );
  fprintf ( 1, '  Each cube face is divided into %dx%d subfaces\n', n, n  );

  point_num = sphere_cubed_grid_point_count ( n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of points is %d\n', point_num );

  xyz = sphere_cubed_grid_points ( n, point_num );

  line_num = sphere_cubed_grid_line_count ( n );

  fprintf ( 1, '  The number of grid lines is %d\n', line_num );

  line_data = sphere_cubed_grid_lines ( n, line_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Line     Start              End\n' );
  fprintf ( 1, '  Index    X    Y   Z         X  Y   Z\n' );
  fprintf ( 1, '\n' );

  for i = 1 : min ( 10, line_num )
    fprintf ( 1, '\n' )
    fprintf ( 1, '  %4d  %10f  %10f  %10f    %10f  %10f  %10f\n', ...
      i, line_data(i,1:3,1), line_data(i,1:3,2) );
  end

  return
end


text:
function [R,T,Xc,best_solution,opt]=efficient_pnp_gauss(x3d_h,x2d_h,A)

% EFFICIENT_PNP_GAUSS Main Function to solve the PnP problem
%       as described in:
%
%       Francesc Moreno-Noguer, Vincent Lepetit, Pascal Fua.
%       Accurate Non-Iterative O(n) Solution to the PnP Problem.
%       In Proceedings of ICCV, 2007.
%
%       Note: In this version of the software we perform a final
%       optimization using Gauss-Newton,which is not described in the
%       paper.
%
%       x3d_h: homogeneous coordinates of the points in world reference
%       x2d_h: homogeneous position of the points in the image plane
%       A: intrincic camera parameters
%       R: Rotation of the camera system wrt world reference
%       T: Translation of the camera system wrt world reference
%       Xc: Position of the points in the camera reference
%       best solution: dimension of the kernel for the best solution
%                     (before applying Gauss Newton).
%       opt: some parameters of the optimization process
%
% Copyright (C) <2007>  <Francesc Moreno-Noguer, Vincent Lepetit, Pascal Fua>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the version 3 of the GNU General Public License
% as published by the Free Software Foundation.
%
% 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/>.
%
% Francesc Moreno-Noguer, CVLab-EPFL, October 2007.
% fmorenoguer@gmail.com, http://cvlab.epfl.ch/~fmoreno/

Xw=x3d_h(:,1:3);
U=x2d_h(:,1:2);

THRESHOLD_REPROJECTION_ERROR=20;%error in degrees of the basis formed by the control points.
%If we have a larger error, we will compute the solution using a larger
%number of vectors in the kernel

%define control points in a world coordinate system (centered on the 3d
%points centroid)
Cw=define_control_points();

%compute alphas (linear combination of the control points to represent the 3d
%points)
Alph=compute_alphas(Xw,Cw);

%Compute M
M=compute_M_ver2(U,Alph,A);

%Compute kernel M
Km=kernel_noise(M,4); %in matlab we have directly the funcion km=null(M);

%1.-Solve assuming dim(ker(M))=1. X=[Km_end];------------------------------
dim_kerM=1;
X1=Km(:,end);
[Cc,Xc,sc]=compute_norm_sign_scaling_factor(X1,Cw,Alph,Xw);

[R,T]=getrotT(Xw,Xc);  %solve exterior orientation
err(1)=reprojection_error_usingRT(Xw,U,R,T,A);

sol(1).Xc=Xc;
sol(1).Cc=Cc;
sol(1).R=R;
sol(1).T=T;
sol(1).error=err(1);
sol(1).betas=[1];
sol(1).sc=sc;
sol(1).Kernel=X1;

%2.-Solve assuming dim(ker(M))=2------------------------------------------
Km1=Km(:,end-1);
Km2=Km(:,end);

%control points distance constraint
D=compute_constraint_distance_2param_6eq_3unk(Km1,Km2);
dsq=define_distances_btw_control_points();
betas_=inv(D'*D)*D'*dsq;
beta1=sqrt(abs(betas_(1)));
beta2=sqrt(abs(betas_(3)))*sign(betas_(2))*sign(betas_(1));
X2=beta1*Km1+beta2*Km2;

[Cc,Xc,sc]=compute_norm_sign_scaling_factor(X2,Cw,Alph,Xw);

[R,T]=getrotT(Xw,Xc);  %solve exterior orientation
err(2)=reprojection_error_usingRT(Xw,U,R,T,A);

sol(2).Xc=Xc;
sol(2).Cc=Cc;
sol(2).R=R;
sol(2).T=T;
sol(2).error=err(2);
sol(2).betas=[beta1,beta2];
sol(2).sc=sc;
sol(2).Kernel=[Km1,Km2];

%3.-Solve assuming dim(ker(M))=3------------------------------------------
if min(err)>THRESHOLD_REPROJECTION_ERROR %just compute if we do not have good solution in the previus cases

    Km1=Km(:,end-2);
    Km2=Km(:,end-1);
    Km3=Km(:,end);

    %control points distance constraint
    D=compute_constraint_distance_3param_6eq_6unk(Km1,Km2,Km3);
    dsq=define_distances_btw_control_points();
    betas_=inv(D)*dsq;
    beta1=sqrt(abs(betas_(1)));
    beta2=sqrt(abs(betas_(4)))*sign(betas_(2))*sign(betas_(1));
    beta3=sqrt(abs(betas_(6)))*sign(betas_(3))*sign(betas_(1));

    X3=beta1*Km1+beta2*Km2+beta3*Km3;

    [Cc,Xc,sc]=compute_norm_sign_scaling_factor(X3,Cw,Alph,Xw);

    [R,T]=getrotT(Xw,Xc);  %solve exterior orientation
    err(3)=reprojection_error_usingRT(Xw,U,R,T,A);

    sol(3).Xc=Xc;
    sol(3).Cc=Cc;
    sol(3).R=R;
    sol(3).T=T;
    sol(3).error=err(3);
    sol(3).betas=[beta1,beta2,beta3];
    sol(3).sc=sc;
    sol(3).Kernel=[Km1,Km2,Km3];

end

%4.-Solve assuming dim(ker(M))=4------------------------------------------
if min(err)>THRESHOLD_REPROJECTION_ERROR %just compute if we do not have good solution in the previus cases
    Km1=Km(:,end-3);
    Km2=Km(:,end-2);
    Km3=Km(:,end-1);
    Km4=Km(:,end);

    D=compute_constraint_distance_orthog_4param_9eq_10unk(Km1,Km2,Km3,Km4);
    dsq=define_distances_btw_control_points();
    lastcolumn=[-dsq',0,0,0]';
    D_=[D,lastcolumn];
    Kd=null(D_);

    P=compute_permutation_constraint4(Kd);
    lambdas_=kernel_noise(P,1);
    lambda(1)=sqrt(abs(lambdas_(1)));
    lambda(2)=sqrt(abs(lambdas_(6)))*sign(lambdas_(2))*sign(lambdas_(1));
    lambda(3)=sqrt(abs(lambdas_(10)))*sign(lambdas_(3))*sign(lambdas_(1));
    lambda(4)=sqrt(abs(lambdas_(13)))*sign(lambdas_(4))*sign(lambdas_(1));
    lambda(5)=sqrt(abs(lambdas_(15)))*sign(lambdas_(5))*sign(lambdas_(1));

    betass_=lambda(1)*Kd(:,1)+lambda(2)*Kd(:,2)+lambda(3)*Kd(:,3)+lambda(4)*Kd(:,4)+lambda(5)*Kd(:,5);
    beta1=sqrt(abs(betass_(1)));
    beta2=sqrt(abs(betass_(5)))*sign(betass_(2));
    beta3=sqrt(abs(betass_(8)))*sign(betass_(3));
    beta4=sqrt(abs(betass_(10)))*sign(betass_(4));
    X4=beta1*Km1+beta2*Km2+beta3*Km3+beta4*Km4;

    [Cc,Xc,sc]=compute_norm_sign_scaling_factor(X4,Cw,Alph,Xw);

    [R,T]=getrotT(Xw,Xc);  %solve exterior orientation
    err(4)=reprojection_error_usingRT(Xw,U,R,T,A);

    sol(4).Xc=Xc;
    sol(4).Cc=Cc;
    sol(4).R=R;
    sol(4).T=T;
    sol(4).error=err(4);
    sol(4).betas=[beta1,beta2,beta3,beta4];
    sol(4).sc=sc;
    sol(4).Kernel=[Km1,Km2,Km3,Km4];
end

%5.-Gauss Newton Optimization------------------------------------------------------
[min_err,best_solution]=min(err);
Xc=sol(best_solution).Xc;
R=sol(best_solution).R;
T=sol(best_solution).T;
Betas=sol(best_solution).betas;
sc=sol(best_solution).sc;
Kernel=sol(best_solution).Kernel;

if best_solution==1
    Betas=[0,0,0,Betas];
elseif best_solution==2
    Betas=[0,0,Betas];
elseif best_solution==3
    Betas=[0,Betas];
end

Km1=Km(:,end-3);
Km2=Km(:,end-2);
Km3=Km(:,end-1);
Km4=Km(:,end);
Kernel=[Km1,Km2,Km3,Km4];

%refine the solution iterating over the betas
Beta0=Betas/sc;
[Xc_opt,R_opt,T_opt,err_opt,iter]=optimize_betas_gauss_newton(Kernel,Cw,Beta0,Alph,Xw,U,A);

%Just update R,T,Xc if Gauss Newton improves results (which is almost
%always)
if err_opt<min_err
    R=R_opt;
    T=T_opt;
    Xc=Xc_opt;
end

opt.Beta0=Beta0;
opt.Kernel=Kernel;
opt.iter=iter;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [R, T]=getrotT(wpts,cpts)

% This routine solves the exterior orientation problem for a point cloud
%  given in both camera and world coordinates.

% wpts = 3D points in arbitrary reference frame
% cpts = 3D points in camera reference frame

n=size(wpts,1);
M=zeros(3);

ccent=mean(cpts);
wcent=mean(wpts);

for i=1:3
  cpts(:,i)=cpts(:,i)-ccent(i)*ones(n,1);
  wpts(:,i)=wpts(:,i)-wcent(i)*ones(n,1);
end
for i=1:n
   M=M+cpts(i,:)'*wpts(i,:);
end
[U S V]=svd(M);
R=U*V';
if det(R)<0
  R=-R;
end
T=ccent'-R*wcent';
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [err,Urep]=reprojection_error_usingRT(Xw,U,R,T,A)

%clear all; close all; load reprojection_error_usingRT;
n=size(Xw,1);

P=A*[R,T];
Xw_h=[Xw,ones(n,1)];
Urep_=(P*Xw_h')';

%project reference points into the image plane
Urep=zeros(n,2);
Urep(:,1)=Urep_(:,1)./Urep_(:,3);
Urep(:,2)=Urep_(:,2)./Urep_(:,3);

%reprojection error
err_=sqrt((U(:,1)-Urep(:,1)).^2+(U(:,2)-Urep(:,2)).^2);
err=sum(err_)/n;

author: cvlab-epfl

repo: EPnP

save_path: github-repos/MATLAB/cvlab-epfl-EPnP


text:
% Kbeta
%
% Computes weighted sum of kernel matrices
%
% Usage:
%	K = Kbeta(Ks,w);
%	K = Kbeta(Ks,w,symmetric);
%
% Input:
%	Ks - Gram matrices [n x n x M]
%	w - weights [M x 1]
%	symmetric (=false) - optional argument. If 1 the Ks
%		matrices are assumed to be symmetric which results in a
%		small speedup.
%
% Output:
%	K - weighted kernel matrix equiv with
%		K=0;for m=1:M, K=K+ w(m)*Ks(:,:,m);end
%
% Example:
%
%	K = Kbeta(randn(100,100,10),rand(10,1));
%
% Peter Gehler 07/2008 pgehler@tuebingen.mpg.de

author: BatzoglouLabSU

repo: SIMLR


text:
close all;
clear all;
clc;
rng('default');

png_export = true;
pdf_export = false;

load ('bin/figure_1_spherical_dict_model_1bp_success_with_k.mat');

mf = spx.graphics.Figures();

mf.new_figure('Recovery probability with K for BP-MMV');
hold all;
legends = cell(1, num_ss);
for ns=1:num_ss
    S = Ss(ns);
    plot(Ks, bp_success_with_k(ns, :));
    legends{ns} = sprintf('S=%d', S);
end
grid on;
xlabel('Sparsity Level');
ylabel('Empirical Recovery Rate');
legend(legends);

author: indigits

repo: sparse-plex


text:
function [phiNorm] = normalize_angle(phi)
%Normalize phi to be between -pi and pi

while(phi>pi)
	phi = phi - 2*pi;
endwhile

while(phi<-pi)
	phi = phi + 2*pi;
endwhile
phiNorm = phi;

end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
function B=evalBSplinePolys(t,x,k,idx)
%%EVALBSPLINEPOLYS Given a set of at least 2*(k-1) knots in non-decreasing
%                  order (some can be repeated), evaluate all associated
%                  scalar order-k b-spline polynomials that might be
%                  nonzero. These are B_{idx-k+1,k}(x) to B_{idx,k}(x). The
%                  first subscript refers to the knots used in determining
%                  the polynomial and the second subscript refers to the
%                  order of the polynomial. The center index of the knots
%                  is idx and must be k-1<=idx<length(t)-(k-1). These
%                  polynomial evaluations can be used in B-spline
%                  interpolation. This function only gives accurate results
%                  if t(idx)<=x<=t(idx+1). Outside of that range, consider
%                  the function evalBSplinePoly, which evaluates the
%                  polynomials one at a time.
%
%INPUTS: t A 1XnumKnots or numKnotsX1 vector containing the knots. The
%          knots must be sorted in ascending order. It is possible for the
%          knots to be repeated. However, it is required that
%          t(idx)<t(idx+1). If more than the minimum number of knots needed
%          for a given order k are passed, then the extra knots will be
%          ignored. Which knots are "extra" depends on the idx input. These
%          must be real.
%        x The numPointsX1 or 1XnumPoints set of points at which the
%          polynomials should be evaluated. It is assumed that
%          t(idx)<=x<=t(idx+1) for all points. Outside of this range, the
%          results are invalid. These must be real.
%        k Optionally, the order of the B-spline polynomials. If this
%          parameter is omitted, then the default of fix(length(t)/2)+1 are
%          used. This is the highest order that can be used with the given
%          number of knots. This is a positive real value.
%      idx This optional input specifies where in terms of the knots the
%          B-spline polynomials are centered. The default if omitted or an
%          empty matrix is passed is k-1. If t is length 2*(k-1), then idx
%          must be k-1 or be omitted. In general idx is limited to
%          k-1<=idx<length(t)-(k-1).
%
%OUTPUTS: B A kXnumPoints matrix contianing the B-spline polynomials values
%           B_{idx-k+1,k}(x) to B_{idx,k}(x) evaluated at each value in x.
%           These values are real.
%
%B-spline polynomials are discussed in detail in Chapter IX of [1]. The are
%most simply expressed by the recursion starting with B_{idx,0}(x)=1 for
%t(idx)<=x<=t(idx+1) (and is zero otherwise) and then
%B_{idx,k}(x)=((x-t(idx))/(t(idx+k)-t(idx)))*B_{idx,k-1}(x)+((t(idx+k+1)-x)/(t(idx+k+1)-t(idx+1)))*B_{idx+1,k-1}(x)
%An algorithm to efficiently implement such a recursion is described in
%Chapter X of [1] and is implemented here. This function has been
%implemented such that when faced with knots such that t(idx)=t(idx+1), all
%of the B values are zero, because the non-finite term is replaced with
%zero. Except in such an instance of repeated knots, sum(B)==1.
%
%REFERENCES:
%[1] C. de Boor, A Practical Guide to Splines. New York: Springer-Verlag,
%    1978.
%
%April 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

kKnots=fix(length(t)/2)+1;

if(nargin<3||isempty(k))
    k=kKnots;
end

if(kKnots<k)
   error('There are not enough points to obtain B-spline polynomials of the desired order.')
end

if(nargin<4||isempty(idx))
   idx=k-1;
end

if(idx<k-1)
    error('The value of iIdx is too small.')
elseif(idx>length(t)-(k-1))
    error('The value of iIdx is too large.')
end

%This is doing if(any(x<t(idx)|x>t(idx+1)) However, finite precision errors
%can make things slightly off and create unnecessary warnings, so the
%comparisons with a tolerance help to avoid that.
epsMult=4;
if(any(-(x-t(idx))>epsMult*eps(x)|(x-t(idx+1))>epsMult*eps(x)))
    warning('Points are outside of the valid region implied by iIdx. Function results can be inaccurate.')
end

x=x(:)';
numPoints=length(x);

B=zeros(k,numPoints);

deltaR=zeros(k-1,numPoints);
deltaL=zeros(k-1,numPoints);

B(1,:)=1;
for j=1:(k-1)
    deltaR(j,:)=t(idx+j)-x;
    deltaL(j,:)=x-t(idx+1-j);

    recurVal=zeros(1,numPoints);
    for r=1:j
        term=B(r,:)./(deltaR(r,:)+deltaL(j+1-r,:));
        B(r,:)=recurVal+deltaR(r,:).*term;

        recurVal=deltaL(j+1-r,:).*term;
    end
    B(j+1,:)=recurVal;
end

end

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


text:
function denoised = NCSR_WRAP(image_path, sigma)

    Noisy_Image = imread(image_path);
    par     =   Parameters_setting( sigma );
    par.nim =   double( Noisy_Image );

    denoised = uint8(255 * mat2gray(NCSR_Denoising( par )));

end

author: lbasek

repo: image-denoising-benchmark


text:
%% Machine Learning Online Class
%  Exercise 8 | Anomaly Detection and Collaborative Filtering
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     estimateGaussian.m
%     selectThreshold.m
%     cofiCostFunc.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% ================== Part 1: Load Example Dataset  ===================
%  We start this exercise by using a small dataset that is easy to
%  visualize.
%
%  Our example case consists of 2 network server statistics across
%  several machines: the latency and throughput of each machine.
%  This exercise will help us find possibly faulty (or very fast) machines.
%

fprintf('Visualizing example dataset for outlier detection.\n\n');

%  The following command loads the dataset. You should now have the
%  variables X, Xval, yval in your environment
load('ex8data1.mat');

%  Visualize the example dataset
plot(X(:, 1), X(:, 2), 'bx');
axis([0 30 0 30]);
xlabel('Latency (ms)');
ylabel('Throughput (mb/s)');

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

%% ================== Part 2: Estimate the dataset statistics ===================
%  For this exercise, we assume a Gaussian distribution for the dataset.
%
%  We first estimate the parameters of our assumed Gaussian distribution,
%  then compute the probabilities for each of the points and then visualize
%  both the overall distribution and where each of the points falls in
%  terms of that distribution.
%
fprintf('Visualizing Gaussian fit.\n\n');

%  Estimate my and sigma2
[mu, sigma2] = estimateGaussian(X);

%  Returns the density of the multivariate normal at each data point (row)
%  of X
p = multivariateGaussian(X, mu, sigma2);

%  Visualize the fit
visualizeFit(X,  mu, sigma2);
xlabel('Latency (ms)');
ylabel('Throughput (mb/s)');

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

%% ================== Part 3: Find Outliers ===================
%  Now you will find a good epsilon threshold using a cross-validation set
%  probabilities given the estimated Gaussian distribution
%

pval = multivariateGaussian(Xval, mu, sigma2);

[epsilon, F1] = selectThreshold(yval, pval);
fprintf('Best epsilon found using cross-validation: %e\n', epsilon);
fprintf('Best F1 on Cross Validation Set:  %f\n', F1);
fprintf('   (you should see a value epsilon of about 8.99e-05)\n\n');

%  Find the outliers in the training set and plot the
outliers = find(p < epsilon);

%  Draw a red circle around those outliers
hold on
plot(X(outliers, 1), X(outliers, 2), 'ro', 'LineWidth', 2, 'MarkerSize', 10);
hold off

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

%% ================== Part 4: Multidimensional Outliers ===================
%  We will now use the code from the previous part and apply it to a
%  harder problem in which more features describe each datapoint and only
%  some features indicate whether a point is an outlier.
%

%  Loads the second dataset. You should now have the
%  variables X, Xval, yval in your environment
load('ex8data2.mat');

%  Apply the same steps to the larger dataset
[mu, sigma2] = estimateGaussian(X);

%  Training set
p = multivariateGaussian(X, mu, sigma2);

%  Cross-validation set
pval = multivariateGaussian(Xval, mu, sigma2);

%  Find the best threshold
[epsilon, F1] = selectThreshold(yval, pval);

fprintf('-> Best epsilon found using cross-validation: %e\n', epsilon);
fprintf('-> Best F1 on Cross Validation Set:  %f\n', F1);
fprintf('# Outliers found: %d\n', sum(p < epsilon));
fprintf('   (you should see a value epsilon of about 1.38e-18)\n\n');
pause

author: worldveil

repo: coursera-ml


text:
function d = sine_taper_scaled(n, k)

% Compute Riedel & Sidorenko sine tapers.
% sine_taper_scaled(n, k) produces the first 2*k tapers of length n,
% returned as the columns of d. The norm of the tapers will not be 1. The
% norm is a function of the number of the taper in the sequence. This is to
% mimick behavior of the scaling of the resulting powerspectra prior to
% april 29, 2011. Before april 29, 2011, equivalent scaling was applied to
% the powerspectra of the tapered data segments, prior to averaging.

% Copyright (C) 2006, Tom Holroyd
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if nargin < 2
  ft_error('usage: sine_taper_scaled(n, k)');
end

k = round(k * 2);
if k <= 0 || k > n
  ft_error('sine_taper_scaled: k is %g, must be in (1:n)/2', k)
end

x = (1:k) .* (pi / (n + 1));
d = sqrt(2 / (n + 1)) .* sin((1:n)' * x);

scale = zeros(k);
for i = 1:k
  scale(i,i) = sqrt(1 - ((i-1)./(k-1)).^2);
end
d = d * scale;

return

author: fieldtrip

repo: fieldtrip


text:
% 3 DOF mobile robot example.
% by Almaskhan Baimyshev
% 08/28/2013

map = struct('name', 'bench_june1.mat', 'start_point', [-14.5 -7.5], 'goal_point', [10 -0.65]);
max_iter = 20e3;
is_benchmark = false;
rand_seed = 40;
variant = 'FNSimple3D';
result = rrt_star(map, max_iter, is_benchmark, rand_seed, variant);

author: olzhas


text:
function pass = test_galleryball( )

names = {'deathstar','gaussian','helmholtz','moire','peaks',...
             'roundpeg','solharm','stripes','wave'};

N = length(names);
% Test construction of each gallery function.
for k = 1:N
    pass(k) = doesNotCrash(names{k});
end

if (nargout > 0)
    pass = all(pass(:));
end
end

function pass = doesNotCrash(name)
try
    fn = cheb.galleryball(name);  % Test returning the function
    pass = true;
catch ME %#ok<NASGU>
    pass = false;
end
end

author: chebfun

repo: chebfun


text:
function [RHSdiv, RHSdivx, RHSdivy] = divergenceTerm2D(F)
% This function calculates the divergence of a field using its face
% average flux vector facevariable, which is a face vector
%
% SYNOPSIS:
%   [RHSdiv, RHSdivx, RHSdivy] = divergenceTerm2D(F)
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% extract data from the mesh structure
Nx = F.domain.dims(1);
Ny = F.domain.dims(2);
G=reshape(1:(Nx+2)*(Ny+2), Nx+2, Ny+2);
DX = repmat(F.domain.cellsize.x(2:end-1), 1, Ny);
DY = repmat(F.domain.cellsize.y(2:end-1)', Nx, 1);

% define the vector of cell index
row_index = reshape(G(2:Nx+1,2:Ny+1),Nx*Ny,1); % main diagonal (only internal cells)

% reassign the east, west, north, and south flux vectors for the
% code readability
Fe = F.xvalue(2:Nx+1,:);		Fw = F.xvalue(1:Nx,:);
Fn = F.yvalue(:,2:Ny+1);       Fs = F.yvalue(:,1:Ny);

% compute the divergence
div_x = (Fe - Fw)./DX;
div_y = (Fn - Fs)./DY;

% define the RHS Vector
RHSdiv = zeros((Nx+2)*(Ny+2),1);
RHSdivx = zeros((Nx+2)*(Ny+2),1);
RHSdivy = zeros((Nx+2)*(Ny+2),1);

% assign the values of the RHS vector
RHSdiv(row_index) = reshape(div_x+div_y,Nx*Ny,1);
RHSdivx(row_index) = reshape(div_x,Nx*Ny,1);
RHSdivy(row_index) = reshape(div_y,Nx*Ny,1);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% change log:
%   040204 smw stub file
%
% file: get_ci_nap_1s.m
%
function varargout = get_ci_nap_1s(varargin)
  error('get_ci_nap_1s is not part of current release.');
return;

author: Sable

repo: mcbench-benchmarks


text:
function [compartments, uniqueCompartments, abbr, uniqueAbbr] = getCompartment(mets)
% Gets the compartment for each metabolite, and the unique compartments,
% from metabolite abbreviation(s), each of which must have compartment
% symbol concatentated on the right hand side (i.e. `metAbbr[*]`).
%
% USAGE:
%
%    [compartments, uniqueCompartments] = getCompartment(mets)
%
% INPUT:
%    mets:                  char array with a single metabolite abbreviation
%                           or
%                           `m x 1` cell array of metabolite abbreviations
%
% OUTPUTS:
%    compartments:          char array with a single compartment identifier
%                           or
%                           `m x 1` cell array of compartment identifiers
%
%    uniqueCompartments:    char array with a single compartment identifier
%                           or
%                           cell array of unique compartment identifiers
%    abbr:                  char array with a single metabolite
%                           abbreviation, without compartment
%                           or
%                           `m x 1` cell array of metabolite abbreviations,
%                           without compartments
%
% .. Author:
%       - Ronan M.T. Fleming
%       - Hulda SH, Nov. 2012   Switched from for loop to regular expression

bool=0;
if ischar(mets)
    tmp{1}=mets;
    mets=tmp;
    bool=1;
end

pat = '(?<abbr>[^\[]+)\[(?<compartment>[^\]]+)\]';
metStruct = regexp(mets,pat,'names'); % m x 1 cell array with fields abbr and compartment in each cell
metStruct = [metStruct{:}]'; % Convert from cell array to double
compartments = {metStruct.compartment}; % Concatenate compartment fields
compartments = reshape(compartments,length(compartments),1);
abbr = {metStruct.abbr}; % Concatenate compartment fields
abbr = reshape(abbr,length(abbr),1);
uniqueCompartments = unique(compartments);
uniqueAbbr = unique(abbr);

if bool==1
    compartments = compartments{1};
    uniqueCompartments = uniqueCompartments{1};
    abbr = abbr{1};
    uniqueAbbr = uniqueAbbr{1};
end

author: opencobra

repo: cobratoolbox


text:
clear, clc, close all

% load a .wav file
[x, fs] = audioread('track.wav');     % get the samples of the .wav file
x = x(:, 1);                        % get the first channel
xmax = max(abs(x));                 % find the maximum abs value
x = x/xmax;                         % scalling the signal

% define analysis parameters
xlen = length(x);                   % length of the signal
wlen = 1024;                        % window length (recomended to be power of 2)
h = wlen/4;                         % hop size (recomended to be power of 2)
nfft = 4096;                        % number of fft points (recomended to be power of 2)

% define the coherent amplification of the window
K = sum(hamming(wlen, 'periodic'))/wlen;

% perform STFT
[s, f, t] = stft(x, wlen, h, nfft, fs);

% take the amplitude of fft(x) and scale it, so not to be a
% function of the length of the window and its coherent amplification
s = abs(s)/wlen/K;

% correction of the DC & Nyquist component
if rem(nfft, 2)                     % odd nfft excludes Nyquist point
    st(2:end, :) = s(2:end, :).*2;
else                                % even nfft includes Nyquist point
    s(2:end-1, :) = s(2:end-1, :).*2;
end

% convert amplitude spectrum to dB (min = -120 dB)
s = 20*log10(s + 1e-6);

% plot the spectrogram
figure(1)
imagesc(t, f, s)
set(gca,'YDir','normal')
set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
xlabel('Time, s')
ylabel('Frequency, Hz')
title('Amplitude spectrogram of the signal')

handl = colorbar;
set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
ylabel(handl, 'Magnitude, dB')

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function gaussian = gaussianFunction
s = 0.1;
gaussian = @(txi,psi) 1/sqrt(2*pi*s^2)*exp(-(txi - psi).^2/(2*s^2));
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% S-SOFM Toolbox
% Version 1.0(beta) 23-Mar-2006
%
% Tested under Matlab 7.0.0.19920 (R14)
%
% This toolbox contains a set of functions and a GUI which can be used to
% create glyphs from spherical Self-Organizing Feature Maps (SOFMs). It is
% freely distributable for educational and research purposes.
%
%
% Example:
%
% load henon-1024-4.mat                   % loads the data to be visualized
% load c4-24.mat                          %loads the S-SOFM structure
% [w,g,r]=trainGlyph(P,X,C,0.5,20,'plot');%trains the S-SOFM
% glyph(X,r);                             % plots the S-SOFM
%
%
%  Graphical User Interface
%     ssofm           - Displays the mail GUI window.
%     ssofmAbout      - Displays the copyright infoprmation.
%
%  List of functions.
%    trig             - Creates/Plots a tessellated shpere.
%    glyph            - Creates/Plots a glyph.
%    trainSphSOFM     - Create/Train a Spherical SOFM.
%    adaptSphSOFM     - Adapt a Spherical SOFM.
%    updateSphSOFM    - Update the weights of a Spherical SOFM.
%    Lcurve           - Estimates the L-curve of Spherical SOFMs.
%    sphereneigh      - Tracks neighboring points on a tessellated shpere.
%    stdrc            - Auxiliary function for range/color calculation.
%    colorscatter     - Plots a 2D or 3D colored scatter diagram.
%
%  Structure files. (mat)
%    c0-1             - Nodes and neighborhood related to trig(0)
%    c1-3             - Nodes and neighborhood related to trig(1)
%    c2-6             - Nodes and neighborhood related to trig(2)
%    c3-12            - Nodes and neighborhood related to trig(3)
%    c4-24            - Nodes and neighborhood related to trig(4)
%
%  Data files. (mat)
%    henon-1024-4     - Henon map on 4 dimensions, 1024 vectors
%    henon-2048-4     - Henon map on 4 dimensions, 2048 vectors
%    henon-ikeda-1024 - Coupled Henon-Ikeda map, 1024 vectros
%    henon-ikeda-4096 - Coupled Henon-Ikeda map, 4096 vectros
%    ikeda-1024-4     - Ikeda map on 4 dimensions, 1024 vectors
%    ikeda-4096-4     - Ikeda map on 4 dimensions, 4096 vectors
%    logistic-1024-4  - Logistic map on 4 dimensions, 1024 vectors
%    logistic-4096-4  - Logistic map on 4 dimensions, 4096 vectors
%
%  Images. (jpg)
%    Icosahedron0     - The original Icosahedron
%    Icosahedron1     - The original Icosahedron subdivided 1 time
%    Icosahedron2     - The original Icosahedron subdivided 2 times
%    Icosahedron3     - The original Icosahedron subdivided 3 times
%    Icosahedron4     - The original Icosahedron subdivided 4 times
%
%
%
% Authors:
%
% Archana P. Sangole, PhD., P.E. (TX chapter)
% School of Physical & Occupational Therapy
% McGill University
% 3654 Promenade Sir-William-Osler
% Montreal, PQ, H3G 1Y5
% e-mail: archana.sangole@mail.mcgill.ca
%
% CRIR, Rehabilitation Institute of Montreal
% 6300 Ave Darlington
% Montreal, PQ, H3S 2J5
% Tel: 514.340.2111 x2188
% Fax: 514.340.2154
%
%
%
% Alexandros Leontitsis, PhD
% Department of Education
% University of Ioannina
% 45110 - Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% Alexandros Leontitsis is grateful to the Greek State Scholarships
% Foundation (IKY) for supporting his work.