text:
function [p,stats] = colamd2 (S, knobs)
%COLAMD2 Column approximate minimum degree permutation.
%    P = COLAMD2(S) returns the column approximate minimum degree permutation
%    vector for the sparse matrix S.  For a non-symmetric matrix S, S(:,P)
%    tends to have sparser LU factors than S.  The Cholesky factorization of
%    S(:,P)'*S(:,P) also tends to be sparser than that of S'*S.  The ordering
%    is followed by a column elimination tree post-ordering.
%
%    Note that this function is the source code for the built-in MATLAB colamd
%    function.  It has been renamed here to colamd2 to avoid a filename clash.
%    colamd and colamd2 are identical.
%
%    See also COLAMD, AMD, SYMAMD, SYMAMD2.
%
%    Example:
%            P = colamd2 (S)
%            [P, stats] = colamd2 (S, knobs)
%
%    knobs is an optional one- to three-element input vector.  If S is m-by-n,
%    then rows with more than max(16,knobs(1)*sqrt(n)) entries are ignored.
%    Columns with more than max(16,knobs(2)*sqrt(min(m,n))) entries are
%    removed prior to ordering, and ordered last in the output permutation P.
%    Only completely dense rows or columns are removed if knobs(1) and knobs(2)
%    are < 0, respectively.  If knobs(3) is nonzero, stats and knobs are
%    printed.  The default is knobs = [10 10 0].  Note that knobs differs from
%    earlier versions of colamd.

%    Copyright 1998-2007, Timothy A. Davis, and Stefan Larimore
%    Developed in collaboration with J. Gilbert and E. Ng.
%
%    Acknowledgements: This work was supported by the National Science
%       Foundation, under grants DMS-9504974 and DMS-9803599.

%-------------------------------------------------------------------------------
% Perform the colamd ordering:
%-------------------------------------------------------------------------------

if (nargout <= 1 & nargin == 1)						    %#ok
    p = colamd2mex (S) ;
elseif (nargout <= 1 & nargin == 2)					    %#ok
    p = colamd2mex (S, knobs) ;
elseif (nargout == 2 & nargin == 1)					    %#ok
    [p, stats] = colamd2mex (S) ;
elseif (nargout == 2 & nargin == 2)					    %#ok
    [p, stats] = colamd2mex (S, knobs) ;
else
    error ('colamd:  incorrect number of input and/or output arguments') ;
end

%-------------------------------------------------------------------------------
% column elimination tree post-ordering:
%-------------------------------------------------------------------------------

[ignore, q] = etree (S (:,p), 'col') ;
p = p (q) ;

%    stats is an optional 20-element output vector that provides data about the
%    ordering and the validity of the input matrix S.  Ordering statistics are
%    in stats (1:3).  stats (1) and stats (2) are the number of dense or empty
%    rows and columns ignored by COLAMD and stats (3) is the number of
%    garbage collections performed on the internal data structure used by
%    COLAMD (roughly of size 2.2*nnz(S) + 4*m + 7*n integers).
%
%    MATLAB built-in functions are intended to generate valid sparse matrices,
%    with no duplicate entries, with ascending row indices of the nonzeros
%    in each column, with a non-negative number of entries in each column (!)
%    and so on.  If a matrix is invalid, then COLAMD may or may not be able
%    to continue.  If there are duplicate entries (a row index appears two or
%    more times in the same column) or if the row indices in a column are out
%    of order, then COLAMD can correct these errors by ignoring the duplicate
%    entries and sorting each column of its internal copy of the matrix S (the
%    input matrix S is not repaired, however).  If a matrix is invalid in other
%    ways then COLAMD cannot continue, an error message is printed, and no
%    output arguments (P or stats) are returned.  COLAMD is thus a simple way
%    to check a sparse matrix to see if it's valid.
%
%    stats (4:7) provide information if COLAMD was able to continue.  The
%    matrix is OK if stats (4) is zero, or 1 if invalid.  stats (5) is the
%    rightmost column index that is unsorted or contains duplicate entries,
%    or zero if no such column exists.  stats (6) is the last seen duplicate
%    or out-of-order row index in the column index given by stats (5), or zero
%    if no such row index exists.  stats (7) is the number of duplicate or
%    out-of-order row indices.
%
%    stats (8:20) is always zero in the current version of COLAMD (reserved
%    for future use).

author: fieldtrip

repo: fieldtrip


text:
function cst = matRad_dummyCst(ct)
% matRad function to create a dummy cst struct for a ct
%
% call
%   cst = matRad_dummyCst(ct)
%
% input
%   ct: matRad ct struct
%
% output
%   cst: matRad cst struct
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

warning('Did not find RTSS. Creating dummy segmentation for matRad.');

% allocate
cst = cell(1,6);

% fill
cst{1,1}          = 0; % first organ has number 0
cst{1,2}          = 'dummyContour';
cst{1,3}          = 'OAR';
cst{1,4}{1}       = find(ct.cubeHU{1}>0.1);
cst{1,5}.Priority = 1;

% set default parameter for biological planning
cst{1,5}.alphaX   = 0.1;
cst{1,5}.betaX    = 0.05;
cst{1,5}.Visible  = 1;
cst{1,5}.visibleColor = [0 0 0];

% define no objcetives
cst{1,6}          = [];

repo: matRad


text:
function h = max(f, g, dim)
%MAX   Maximum value of a CHEBFUN3 in one of the directions.
%   MAX(F) returns a CHEBFUN2 representing the maximum of the CHEBFUN3
%   object F along the x direction, i.e, MAX(F) = @(y,z) max(F(:, y, z)).
%
%   MAX(F, [], dim) returns a CHEBFUN2 representing the maximum of F along
%   the dimension DIM, where DIM = 1 means along the x-direction, DIM = 2
%   for the y-direction, and DIM = 3 means along the z-direction.
%
%   WARNING: This function is not always accurate to the expected
%   precision.
%
%   Use MAX3 for the global maximum.
%
% See also CHEBFUN3/MAX2 and CHEBFUN3/MAX3.

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

% Empty check:
if ( isempty(f) )
    error('CHEBFUN:CHEBFUN3:max:input', 'CHEBFUN3 is empty');
end

% Default to max of one chebfun3:
if ( nargin < 2 )
    g = [];
end

% Default to maximum along the x direction:
if ( nargin < 3 )
    dim = 1;
end

% Do not allow max(F, G):
if ( nargin > 1 && ~isempty(g) )
    error('CHEBFUN:CHEBFUN3:max:twoCHEBFUN3Inputs', ...
        'Unable to maximize two CHEBFUN3 objects.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We have no idea how to achieve this in an efficient way. This
% is an attempt to return a result, but typically it won't be accurate to
% more than 4-5 digits.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dom = f.domain;
n = 129;
if ( dim == 1 )
    vals = sample(f, n, n, n);
    h = chebfun2(squeeze(max(vals, [], 1)), dom(3:6));
    h = simplify(h);
elseif ( dim == 2 )
    vals = sample(f, n, n, n);
    h = chebfun2(squeeze(max(vals, [], 2)), [dom(1:2), dom(5:6)]);
    h = simplify(h);
elseif ( dim == 3 )
    vals = sample(f, n, n, n);
    h = chebfun2(max(vals, [], 3), dom(1:4));
    h = simplify(h);
elseif ( dim == 0 )
    error('CHEBFUN:CHEBFUN3:max:dim', ...
        'Dimension argument must be a positive integer scalar within indexing range.')
else
   % return the CHEBFUN3 object. This is analogous to that MAX() command in
   % MATLAB.
   h = f;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% gcontour
% Below is a demonstration of the features of the |gcontour| function

%%
clear; close all; clc;

%% Syntax
% |[varargout]=gcontour(varargin);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function data = PackMinibatchTensorVariableLength(feat, nSeqInMB, nSeq, isVAD)
nBatch = ceil(nSeq/nSeqInMB);
nStream = size(feat,1);

if size(feat,2)<nSeq    % some cell have more than 1 sequence. we need to make each cell contains just 1 seq first.
    feat2 = cell(size(feat,1), nSeq);
    % to be implemented later

end

randSegIdx = randperm(nSeq);
feat(:,randSegIdx) = feat;
for i=1:nStream
    if isVAD(i); continue; end
    for j=1:nBatch
        curr_feat = feat(i, nSeqInMB*(j-1)+1 : min(nSeq,nSeqInMB*j));
        data{i,j} = cell2mat_tensor3D(curr_feat, -1e10);
    end
end

end

author: singaxiong

repo: SignalGraph


text:
% Chapter 14 - Poincare Maps and Nonautonomous Systems in the Plane.
% Programs_14c - Phase portraits for nonautonomous systems.
% Copyright Birkhauser 2013. Stephen Lynch.

% Phase portrait for the Duffing system.

deq=@(t,x) [x(2);x(1)-0.3*x(2)-(x(1))^3+0.5*cos(1.25*t)];
options=odeset('RelTol',1e-4,'AbsTol',1e-4);
[t,xx]=ode45(deq,[0 200],[1,0],options);

plot(xx(:,1),xx(:,2))

fsize=15;
axis([-2 2 -2 2])
xlabel('x','FontSize',fsize)
ylabel('y','FontSize',fsize)

% End of Programs_14c.

author: Sable

repo: mcbench-benchmarks


text:
function [EEG,ref,latency,snr] = generate_synth(EEG,options)
%   GENERATE_SYNTH
%       [EEG,REF,LATENCY,SNR] = GENERATE_SYNTH(EEG,OPTIONS)
%
%   Created by Alexandre Gramfort on 2009-05-10.
%   Copyright (c)  Alexandre Gramfort. All rights reserved.

% $Id: generate_synth.m 2 2009-06-16 19:24:10Z gramfort $
% $LastChangedBy: gramfort $
% $LastChangedDate: 2009-06-16 15:24:10 -0400 (Mar, 16 jui 2009) $
% $Revision: 2 $

if nargin<2
    options.null = 0;
end

if ~isfield(options, 'mean_latency')
    options.mean_latency = 0.3;
end
mean_latency = options.mean_latency;

if ~isfield(options, 'mean_sigma')
    options.mean_sigma = 0.05;
end
mean_sigma = options.mean_sigma;

if ~isfield(options, 'mean_frequency')
    options.mean_frequency = 5;
end
mean_frequency = options.mean_frequency;

if ~isfield(options, 'mean_strech')
    options.mean_strech = 0;
end
mean_strech = options.mean_strech;

if ~isfield(options, 'std_latency')
    options.std_latency = 0;
end
std_latency = options.std_latency;

if ~isfield(options, 'std_sigma')
    options.std_sigma = 0;
end
std_sigma = options.std_sigma;

if ~isfield(options, 'std_frequency')
    options.std_frequency = 0;
end
std_frequency = options.std_frequency;

if ~isfield(options, 'std_strech')
    options.std_strech = 0;
end
std_strech = options.std_strech;

if ~isfield(options, 'snr')
    options.snr = 2;
end
snr = options.snr;

if ~isfield(options, 'use_ar')
    options.use_ar = false;
end
use_ar = options.use_ar;

if isfield(options, 'null')
    options = rmfield(options,'null');
end

nepochs = EEG.trials;

R = TruncatedGaussian(1,[-2 2],[nepochs,1]);

latency   = mean_latency   + std_latency * R;
sigma     = mean_sigma     + std_sigma * R;
frequency = mean_frequency + std_frequency * R;
strech    = mean_strech    + std_strech * R;

% fname_suffix = [ '_', ...
%                  'mfreq_',num2str(mean_frequency),'_', ...
%                  'stdfreq_',num2str(std_frequency),'_', ...
%                  'mlat_',num2str(mean_latency),'_', ...
%                  'stdlat_',num2str(std_latency),'_', ...
%                  'msigma_',num2str(mean_sigma),'_', ...
%                  'stdsigma_',num2str(std_sigma),'_', ...
% ];

% noise AR filter config
b = fir1(1024, .5);
[d,p0] = lpc(b,7);
noise = zeros([1 EEG.pnts EEG.trials]);

dphase = pi/3;
gabor = @(time,freq,sigma,dphase) cos(2*pi*freq*time - dphase) .* exp(-time.^2./2./sigma.^2);
ref = gabor(EEG.times./ (1 + mean_strech)-mean_latency,mean_frequency,mean_sigma,dphase);
ref = ref ./ norm(ref);

for iepoch = 1:nepochs
    time_epoch = EEG.times ./ (1 + strech(iepoch)) - latency(iepoch);
    EEG.data(1,:,iepoch) = gabor(time_epoch,frequency(iepoch),sigma(iepoch),dphase);
    EEG.data(1,:,iepoch) = EEG.data(1,:,iepoch) ./ norm(EEG.data(1,:,iepoch));

    if use_ar
        noise(1,:,iepoch) =  filter(1,d,sqrt(p0)*randn(EEG.pnts,1))'; % AR noise
    else
        noise(1,:,iepoch) =  randn(1,size(noise,2)); % white noise
    end
    noise(1,:,iepoch) = noise(1,:,iepoch) ./ norm(noise(1,:,iepoch));
end

noise = noise ./ norm(noise(:));
noise = norm(EEG.data(:)) * noise ./ snr;

EEG.data = EEG.data + noise;

EEG.times = EEG.times;
EEG.setname = 'synth_realignment';
EEG.xmin = 0;
EEG.xmax = 1;
EEG.icaact = [];
EEG.icawinv = [];
EEG.icasphere = [];
EEG.icaweights = [];

end %  function

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [x, funVal, ValueL]=tree_mcLeastR(A, y, z, opts)
%
%%
% Function mcLeastR:
%      Least Squares Loss for Multi-class (task) Learning
%             with the tree structured group Lasso Regularization
%
%% Problem
%
%  min  1/2 || A x - y||^2 + z * sum_i sum_j w_j ||x^i_{G_j}||
%
%  x^i denotes the i-th row of x
%
%  G_j's are nodes with tree structure
%
%  We assume that the tasks are of a tree structure
%
%  The tree structured group information is contained in
%  opts.ind, which is a 3 x nodes matrix, where nodes denotes the number of
%  nodes of the tree.
%
%  opts.ind(1,:) contains the starting index
%  opts.ind(2,:) contains the ending index
%  opts.ind(3,:) contains the corresponding weight (w_j)
%
%  Note:
%  1) If each element of x^j is a leaf node of the tree and the weight for
%  this leaf node are the same, we provide an alternative "efficient" input
%  for this kind of node, by creating a "super node" with
%  opts.ind(1,1)=-1; opts.ind(2,1)=-1; and opts.ind(3,1)=the common weight.
%
%  2) If the features are well ordered in that, the features of the left
%  tree is always less than those of the right tree, opts.ind(1,:) and
%  opts.ind(2,:) contain the "real" starting and ending indices. That is to
%  say, x^j( opts.ind(1,j):opts.ind(2,j) ) denotes x^j_{G_j}. In this case,
%  the entries in opts.ind(1:2,:) are within 1 and k (the number of tasks).
%
%
%  If the features are not well ordered, please use the input opts.G for
%  specifying the index so that
%   x^j( opts.G ( opts.ind(1,j):opts.ind(2,j) ) ) denotes x^j_{G_j}.
%  In this case, the entries of opts.G are within 1 and n, and the entries of
%  opts.ind(1:2,:) are within 1 and length(opts.G).
%
% The following example shows how G and ind works:
%
% G={ {1, 2}, {4, 5}, {3, 6}, {7, 8},
%     {1, 2, 3, 6}, {4, 5, 7, 8},
%     {1, 2, 3, 4, 5, 6, 7, 8} }.
%
% ind={ [1, 2, 100]', [3, 4, 100]', [5, 6, 100]', [7, 8, 100]',
%       [9, 12, 100]', [13, 16, 100]', [17, 24, 100]' },
%
% where each node has a weight of 100.
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size mxk)
%  z -        Regularization parameter (z >=0)
%  opts-      optional inputs (default value: opts=[])
%
%% Output parameters:
%  x-         Solution (of size n x k)
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on October 3, 2010.
%
%% Related papers
%
% [1] Jun Liu and Jieping Ye, Moreau-Yosida Regularization for
%     Grouped Tree Structure Learning, NIPS 2010
%
%% Related functions:
%
%  sll_opts
%
%%

%% Verify and initialize the parameters
%%
if (nargin <4)
    error('\n Inputs: A, y, z, and opts (.ind) should be specified!\n');
end

[m,n]=size(A);

k=size(y,2);

if (size(y,1) ~=m)
    error('\n Check the length of y!\n');
end

if (z<0)
    error('\n z should be nonnegative!\n');
end

opts=sll_opts(opts); % run sll_opts to set default values (flags)

% restart the program for better efficiency
%  this is a newly added function
if (~isfield(opts,'rStartNum'))
    opts.rStartNum=opts.maxIter;
else
    if (opts.rStartNum<=0)
        opts.rStartNum=opts.maxIter;
    end
end

%% Detailed initialization
%%
%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Group & Others

% Initialize ind
if (~isfield(opts,'ind'))
    error('\n In tree_mcLeastR, the field .ind should be specified');
else
    ind=opts.ind;

    if (size(ind,1)~=3)
        error('\n Check opts.ind');
    end
end

GFlag=0;
% if GFlag=1, we shall apply general_altra
if (isfield(opts,'G'))
    GFlag=1;

    G=opts.G;
    if (max(G) >k || max(G) <1)
        error('\n The input G is incorrect. It should be within %d and %d',1,k);
    end
end

%% Starting point initialization

% compute AT y
if (opts.nFlag==0)
    ATy =A'* y;
elseif (opts.nFlag==1)
    ATy= A'* y -  mu' * sum(y, 1);    ATy=ATy./repmat(nu, 1, k);
else
    invNu=y./repmat(nu, 1, k);        ATy=A'*invNu-mu' * sum(invNu, 1);
end

% process the regularization parameter
if (opts.rFlag==0)
    lambda=z;
else % z here is the scaling factor lying in [0,1]
%     if (z<0 || z>1)
%         error('\n opts.rFlag=1, and z should be in [0,1]');
%     end

    computedFlag=0;
    if (isfield(opts,'lambdaMax'))
        if (opts.lambdaMax~=-1)
            lambda=z*opts.lambdaMax;
            computedFlag=1;
        end
    end

    if (~computedFlag)
        if (GFlag==0)
            lambda_max=findLambdaMax_mt(ATy, n, k, ind, size(ind,2));
        else
            lambda_max=general_findLambdaMax_mt(ATy, n, k, G, ind, size(ind,2));
        end

        % As .rFlag=1, we set lambda as a ratio of lambda_max
        lambda=z*lambda_max;
    end
end

% The following is for computing lambdaMax
% we use opts.lambdaMax=-1 to show that we need the computation.
%
% One can use this for setting up opts.lambdaMax
if (isfield(opts,'lambdaMax'))

    if (opts.lambdaMax==-1)
        if (GFlag==0)
            lambda_max=findLambdaMax_mt(ATy, n, k, ind, size(ind,2));
        else
            lambda_max=general_findLambdaMax_mt(ATy, n, k, G, ind, size(ind,2));
        end

        x=lambda_max;
        funVal=lambda_max;
        ValueL=lambda_max;

        return;
    end
end

% initialize a starting point
if opts.init==2
    x=zeros(n,k);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if ( size(x,1)~=n || size(x,2)~=k )
            error('\n Check the input .x0');
        end
    else
        x=ATy;  % if .x0 is not specified, we use ratio*ATy,
        % where ratio is a positive value
    end
end

% compute Ax= A * x
if (opts.nFlag==0)
    Ax=A* x;
elseif (opts.nFlag==1)
    invNu=x./repmat(nu, 1, k);  mu_invNu=mu * invNu;
    Ax=A*invNu -repmat(mu_invNu, m, 1);
else
    Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./repmat(nu, 1, k);
end

if (opts.init==0)
    % ------  This function is not available
    %
    % If .init=0, we set x=ratio*x by "initFactor"
    % Please refer to the function initFactor for detail
    %

    % Here, we only support starting from zero, due to the complex tree
    % structure

    x=zeros(n,k);
end

%% The main program

bFlag=0; % this flag tests whether the gradient step only changes a little

% The Armijo Goldstein line search scheme + accelearted gradient descent
if (opts.mFlag==0 && opts.lFlag==0)

    L=1;
    % We assume that the maximum eigenvalue of A'A is over 1

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,k);

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute the gradient (g) at s
        As=Ax + beta* (Ax-Axp);

        % compute AT As
        if (opts.nFlag==0)
            ATAs=A'*As;
        elseif (opts.nFlag==1)
            ATAs= A'* As -  mu' * sum(As, 1);  ATAs=ATAs./repmat(nu, 1, k);
        else
            invNu=As./repmat(nu, 1, k);        ATAs=A'*invNu-mu' * sum(invNu, 1);
        end

        % obtain the gradient g
        g=ATAs-ATy;

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the L1/Lq-norm regularized projection
            v=s-g/L;

            % tree overlapping group Lasso projection
            ind_work(1:2,:)=ind(1:2,:);
            ind_work(3,:)=ind(3,:) * (lambda / L);

            if (GFlag==0)
                x=altra_mt(v, n, k, ind_work, size(ind_work,2));
            else
                x=general_altra_mt(v, n, k, G, ind_work, size(ind_work,2));
            end

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute Ax= A * x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./repmat(nu, 1, k);  mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./repmat(nu, 1, k);
            end

            Av=Ax -As;
            r_sum=norm(v,'fro')^2; l_sum=norm(Av,'fro')^2;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is ||Av||_2^2 <= L * ||v||_2^2
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end

        ValueL(iterStep)=L;

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        xxp=x-xp;   Axy=Ax-y;
        funVal(iterStep)=norm(Axy,'fro')^2/2;

        for i=1:n
            xRow=x(i,:);

            if (GFlag==0)
                tree_norm=treeNorm(xRow, k, ind, size(ind,2));
            else
                tree_norm=general_treeNorm(xRow, k, G, ind, size(ind,2));
            end

            funVal(iterStep)=funVal(iterStep)+ lambda*tree_norm;
        end

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=norm(xp,'fro');    norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end

         % restart the program every opts.rStartNum
        if (~mod(iterStep, opts.rStartNum))
            alphap=0; alpha=1;
            xp=x; Axp=Ax; xxp=zeros(n,k); L =L/2;
        end
    end
else
    error('\n The function does not support opts.mFlag neq 0 & opts.lFlag neq 0!');
end

author: thomaskuestner


text:
function CPD = reset_ess(CPD)
% RESET_ESS Reset the Expected Sufficient Statistics of a hhmm F node.
% CPD = reset_ess(CPD)

CPD.sub_CPD_term = reset_ess(CPD.sub_CPD_term);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function s = movingstd(x,k,windowmode)
% movingstd: efficient windowed standard deviation of a time series
% usage: s = movingstd(x,k,windowmode)
%
% Movingstd uses filter to compute the standard deviation, using
% the trick of std = sqrt((sum(x.^2) - n*xbar.^2)/(n-1)).
% Beware that this formula can suffer from numerical problems for
% data which is large in magnitude.
%
% At the ends of the series, when filter would generate spurious
% results otherwise, the standard deviations are corrected by
% the use of shorter window lengths.
%
% arguments: (input)
%  x   - vector containing time series data
%
%  k   - size of the moving window to use (see windowmode)
%        All windowmodes adjust the window width near the ends of
%        the series as necessary.
%
%        k must be an integer, at least 1 for a 'central' window,
%        and at least 2 for 'forward' or 'backward'
%
%  windowmode - (OPTIONAL) flag, denotes the type of moving window used
%        DEFAULT: 'central'
%
%        windowmode = 'central' --> use a sliding window centered on
%            each point in the series. The window will have total width
%            of 2*k+1 points, thus k points on each side.
%
%        windowmode = 'backward' --> use a sliding window that uses the
%            current point and looks back over a total of k points.
%
%        windowmode = 'forward' --> use a sliding window that uses the
%            current point and looks forward over a total of k points.
%
%        Any simple contraction of the above options is valid, even
%        as short as a single character 'c', 'b', or 'f'. Case is
%        ignored.
%
% arguments: (output)
%  s   - vector containing the windowed standard deviation.
%        length(s) == length(x)

% check for a windowmode
if (nargin<3) || isempty(windowmode)
  % supply the default:
  windowmode = 'central';
elseif ~ischar(windowmode)
  error 'If supplied, windowmode must be a character flag.'
end
% check for a valid shortening.
valid = {'central' 'forward' 'backward'};
windowmode = lower(windowmode);
ind = strmatch(windowmode,valid);
if isempty(ind)
  error 'Windowmode must be a character flag: ''c'', ''b'', or ''f''.'
else
  windowmode = valid{ind};
end

% length of the time series
n = length(x);

% check for valid k
if (nargin<2) || isempty(k) || (rem(k,1)~=0)
  error 'k was not provided or not an integer.'
end
switch windowmode
  case 'central'
    if k<1
      error 'k must be at least 1 for windowmode = ''central''.'
    end
    if n<(2*k+1)
      error 'k is too large for this short of a series and this windowmode.'
    end
  otherwise
    if k<2
      error 'k must be at least 2 for windowmode = ''forward'' or ''backward''.'
    end
    if (n<k)
      error 'k is too large for this short of a series.'
    end
end

% Improve the numerical analysis by subtracting off the series mean
% this has no effect on the standard deviation.
x = x - mean(x);

% we will need the squared elements
x2 = x.^2;

% split into the three windowmode cases for simplicity
A = 1;
switch windowmode
  case 'central'
    B = ones(1,2*k+1);
    s = sqrt((filter(B,A,x2) - (filter(B,A,x).^2)*(1/(2*k+1)))/(2*k));
    s(k:(n-k)) = s((2*k):end);
  case 'forward'
    B = ones(1,k);
    s = sqrt((filter(B,A,x2) - (filter(B,A,x).^2)*(1/k))/(k-1));
    s(1:(n-k+1)) = s(k:end);
  case 'backward'
    B = ones(1,k);
    s = sqrt((filter(B,A,x2) - (filter(B,A,x).^2)*(1/k))/(k-1));
end

% special case the ends as appropriate
switch windowmode
  case 'central'
    % repairs are needed at both ends
    for i = 1:k
      s(i) = std(x(1:(k+i)));
      s(n-k+i) = std(x((n-2*k+i):n));
    end
  case 'forward'
    % the last k elements must be repaired
    for i = (k-1):-1:1
      s(n-i+1) = std(x((n-i+1):n));
    end
  case 'backward'
    % the first k elements must be repaired
    for i = 1:(k-1)
      s(i) = std(x(1:i));
    end
end


text:
% Function to do creaky voice detection using the algorithm described in
% Ishi et al 2008
%
% Description
% Function to do creaky voice detection using the algorithm described in
% Ishi et al 2008
%
%
% Inputs
%  x        : [samples] [Nx1]  Speech signal
%  fs       : [Hz]      [1x1]  Sampling frequency
%  plots    : [integer] [1x1]  Option for plots (1==do plot)
%
% Outputs
%
% Example
%  Please see the HOWTO_glottalsource.m example file.
%
% References
%  [1] Ishi, C., Sakakibara, K-I, Ishiguro, H., (2008) `A method for
%       automatic detection of vocal fry', IEEE TASLP, 16(1), 47-56.
%  [2] Drugman, T., Kane, J., Gobl, C., `Automatic Analysis of Creaky
%       Excitation Patterns', Submitted to Computer Speech and
%       Language.
%  [3] Kane, J., Drugman, T., Gobl, C., (2013) `Improved automatic
%       detection of creak', Computer Speech and Language 27(4), pp.
%       1028-1047.
%  [4] Drugman, T., Kane, J., Gobl, C., (2012) `Resonator-based creaky
%       voice detection', Interspeech 2012, Portland, Oregon, USA.
%
% Copyright (c) 2013 University of Mons, FNRS & 2013 Trinity College Dublin
%
% License
%  This code is a part of the GLOAT toolbox with the following
%  licence:
%  This program is free software: you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation, either version 3 of the License, or
%  (at your option) any later version.
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Authors
%  Thomas Drugman <thomas.drugman@umons.ac.be> & John Kane <kanejo@tcd.ie>

function [creak,t_IFP,IFP,PwP,IPS,t_pow,veryShort_powCont_dB,peak_idx,peak_pow,peak_idx_samp] = ...
    ishi_creak_detection(x,fs,plots)

%% Thresholds from Ishi et al (2008)
localpeak_thresh = 2; % 2 dB threshold,
PwPthresh = 7; % Power peaks above 7 dB for creak
IFPthresh = 0.5; % Infraframe perdiocity below 0.5 for creak
IPSthresh = 0.5; % Interpulse similarity above 0.5 for creak
Tmax = 100*(fs/1000); % Maximum interpulse interval 100 ms (i.e. F0 = 10 Hz)

%% Bandpass filter signal 100 ~ 1500 Hz
Rp = 3; Rs = 40;
Wp = [100 1500]/(fs/2);
Ws = [50 3000]/(fs/2);
[n,Wn] = buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(n,Wn);
x_filt=filtfilt(b,a,x); % Zero-phase filtering

%% (1) Get PwP values from "very short-term" power contour
veryShort_len = 4*(fs/1000); % 4ms frame length for "very short-term" analysis
veryShort_shift = 2*(fs/1000); % 2ms shift for "very short-term" analysis
veryShort_powCont = zeros(1,ceil((length(x)-veryShort_len)/veryShort_shift));
t_pow = zeros(1,ceil((length(x)-veryShort_len)/veryShort_shift));
start=1;
finish=start+veryShort_len-1;

% Get very short term power contour
n=1;
x_filt2 = x_filt.^2;
while finish <= length(x_filt)
    veryShort_powCont(n) = mean(x_filt2(start:finish));
    t_pow(n)=mean([start finish]);
    start = start + veryShort_shift;
    finish=start+veryShort_len-1;
    n=n+1;
end
clear x_filt2;

veryShort_powCont_dB = 20*log10(veryShort_powCont);

% Find peaks in power contour
[peak_idx,peak_pow]=v_findpeaks(veryShort_powCont_dB);

for n=1:length(peak_pow)

    if peak_idx(n)-3 > 0
        start=peak_idx(n)-3;
    else start=1;
    end
    if peak_idx(n)+3<=length(veryShort_powCont_dB)
        finish=peak_idx(n)+3;
    else finish=length(veryShort_powCont_dB);
    end

    if peak_pow(n)-localpeak_thresh < veryShort_powCont_dB(start) || ...
            peak_pow(n)-localpeak_thresh < veryShort_powCont_dB(finish)
        peak_pow(n)=NaN;
        peak_idx(n)=NaN;
    end
end

peak_pow(isnan(peak_pow))=[];
peak_idx(isnan(peak_idx))=[];

% Suggestion from Ishi
maxPow=max(peak_pow);
max2remove_dB=30;
peak_idx(peak_pow<maxPow-max2remove_dB)=[];
peak_pow(peak_pow<maxPow-max2remove_dB)=[];

% Get power rising and power falling values
PwP.rise =zeros(1,length(peak_idx));
PwP.fall =zeros(1,length(peak_idx));
PwP.idx=peak_idx;

for n=1:length(peak_idx)
    if peak_idx(n)-5 > 0
        start = peak_idx(n)-5;
    else start = 1;
    end
    if peak_idx(n)+5 <= length(veryShort_powCont_dB)
        finish = peak_idx(n)+5;
    else finish = length(veryShort_powCont_dB);
    end

    % power rising
    PwP.rise(n) = max(peak_pow(n)- ...
        veryShort_powCont_dB(start:peak_idx(n)-1));
    % power falling
    PwP.fall(n) = max(peak_pow(n)- ...
        veryShort_powCont_dB(peak_idx(n)+1:finish));
end

%% (2) Calculate Intraframe Periodicity Measure (IFP)
peak_idx_samp=round(t_pow(peak_idx));
[IFP,t_IFP,IFP_peaks] = get_ifp(x_filt,fs,IFPthresh,t_pow,peak_idx);

%% (3) Calculate interpulse similarity measure (IPS)
IPS = zeros(1,length(peak_idx));
N = 15*(fs/1000); % 15 ms frame length for similarity analysis
maxLags = 15*(fs/1000); % 15 ms frame length for similarity analysis

for n=2:length(IPS);
    if peak_idx_samp(n)-peak_idx_samp(n-1) > Tmax
        IPS(n) = 0;
    else

        if peak_idx_samp(n)+round(N/2)>length(x)
            break
        end
        if peak_idx_samp(n-1)-round(N/2) < 1 || peak_idx_samp(n)+round(N/2) > length(x_filt)
            IPS(n) = 0;
        else
            peak_seg_cur = x_filt(peak_idx_samp(n)-round(N/2):peak_idx_samp(n)+round(N/2));
            peak_seg_prev = x_filt(peak_idx_samp(n-1)-round(N/2):peak_idx_samp(n-1)+round(N/2));
            Rxx = xcorr(peak_seg_cur,peak_seg_prev,maxLags,'coeff'); % Normalised cross correlation
            IPS(n)=max(Rxx);
        end
    end
end

%% Decide on candidate peaks based on set threholds
creak_peak_idx=zeros(1,length(peak_idx_samp));
creak_peak_idx(PwP.rise>=PwPthresh & PwP.fall>=PwPthresh & ...
    IFP_peaks<=IFPthresh & IPS>=IPSthresh)=1; % binary vector
creak=zeros(1,length(x));
creak_dur_ms=[];
n=1;
cnt=1;
while n < length(creak_peak_idx)
    if creak_peak_idx(n)==1 && creak_peak_idx(n+1)==1
        start = peak_idx_samp(n);
        n=n+1;
        while creak_peak_idx(n)==1 && peak_idx_samp(n) ...
                - peak_idx_samp(n-1) < Tmax && n < length(creak_peak_idx)
            n=n+1;
        end
        finish= peak_idx_samp(n-1);

        creak(start:finish)=1;
        cnt=cnt+1;

    end
    n=n+1;
end

%creak=creak.*VUV; % Remove creak detected in unvoiced areas
IFP(IFP<0)=0;
IFP(IFP>1)=1;
%% Plots
if nargin > 2
    if plots==1
        t=(1:length(x))/fs;
       % figure,
       x=x/max(x);
        subplot(411), plot(t,x/max(x)), hold on, plot(t,(creak)*.8,'r'), hold off
         axis tight, ylabel('Amplitude','FontSize',14)
         title('Speech waveform with creak regions','FontSize',14,'FontWeight','bold'), ylim([min(x)*1.1 max(x)*1.1])
        subplot(412), plot(t(round(t_pow)),veryShort_powCont_dB), hold on,
         plot(t(peak_idx_samp),peak_pow,'rx'),
         plot(t(peak_idx_samp(PwP.rise>PwPthresh)),peak_pow(PwP.rise>PwPthresh),'gx'),
         hold off, axis tight,
         ylabel('Power (dB)','FontSize',14)
         title('Very short power contour (dB) with Local peaks','FontSize',14,'FontWeight','bold')
        subplot(413), plot(t(t_IFP),IFP),hold on,% plot([0 length(x)], ...
           % [IFPthresh IFPthresh],'k'), plot(t_IFP,IFP_orig,'r')
            hold off, axis tight, title('Intraframe periodicity','FontSize',14,'FontWeight','bold')
            ylabel('IFP','FontSize',14)
        subplot(414), stem(t(peak_idx_samp),IPS), hold on,  plot(t([1 length(x)]), ...
            [IPSthresh IPSthresh],'k'), hold off, axis tight,
        ylabel('IPS','FontSize',14),xlabel('Time (seconds)','FontSize',14)
        title('Interpulse similarity','FontSize',14,'FontWeight','bold')
    end
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function yaw_start(nlp, bounds)

    domain = nlp.Plant;

     % yaw initialized zero
    x = domain.States.x;
    yaw = x('BaseRotZ');
    yaw_fun = SymFunction(['yaw_',domain.Name], yaw, {x});
    addNodeConstraint(nlp, yaw_fun, {'x'}, 'first',  ...
        bounds.constrBounds.yaw_initial.lb, ...
        bounds.constrBounds.yaw_initial.ub,'Linear');

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function value = i4vec_any_lt ( n, a, b )

%*****************************************************************************80
%
%% I4VEC_ANY_LT: ( any ( A < B ) ) for I4VEC's.
%
%  Discussion:
%
%    An I4VEC is a vector of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 April 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries.
%
%    Input, integer A(N), the first vector.
%
%    Input, integer B(N), the second vector.
%
%    Output, logical I4VEC_ANY_LT is 1 if any entry
%    of A is less than the corresponding entry of B.
%
  value = any ( a(1:n) < b(1:n) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ cphd_meta ] = cphd_update_meta( cphd_meta, version_str )
%CPHD_UPDATE_META Update a CPHD metadata structure to latest
%
% Basically a cut-and-paste from sicd_update_meta.
%
% Author: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Update metadata structure to version 1.0 if necessary
vers_parts = str2double(regexp(version_str,'\.','split'));
if all(isfinite(vers_parts)) % Valid version format
    if compare_version(vers_parts, 1) % Version 1.0
        % do nothing
    elseif compare_version(vers_parts, [0 3]) % Version 0.3
        cphd_meta = cphd_update_meta_0_3(cphd_meta);
    else % Either older version or misformed version string
        error('CPHD_UPDATE_META:CPHD_VERSION','Unrecognized CPHD version.');
    end
end

end

% Is vers_ref is new than or equal to vers_comapre
function bool_out = compare_version(vers_ref, vers_compare)
    if numel(vers_ref)>numel(vers_compare)
        vers_compare(numel(vers_ref))=0;
    elseif numel(vers_compare)>numel(vers_ref)
        vers_ref(numel(vers_compare))=0;
    end
    tenpowers = 10.^(0:-1:-(numel(vers_ref)-1)).';
    bool_out = (sign(vers_compare - vers_ref) * tenpowers) <= 0;
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function img = bb_burn(img, bb_in, width)

if ~exist('width','var')
    width = 0;
end

if isempty(bb_in)
    return;
end

for i = 1:size(bb_in,2)

    bb = bb_in(:,i);
    bb(1) = max([1 bb(1)]);
    bb(2) = max([1 bb(2)]);
    bb(3) = min([size(img,2), bb(3)]);
    bb(4) = min([size(img,1), bb(4)]);

    img(bb(2):bb(2)+width,bb(1):bb(3)) = 255;
    img(bb(4)-width:bb(4),bb(1):bb(3)) = 255;
    img(bb(2):bb(4),bb(1):bb(1)+width) = 255;
    img(bb(2):bb(4),bb(3)-width:bb(3)) = 255;

end

repo: OpenTLD


text:
%PRDOWNLOAD Download well defined data and toolboxes
%
%   STATUS = PRDOWNLOAD(URL,DIRNAME)
%
% INPUT
%   URL          String containing URL of file to be downloaded
%   DIRNAME      Final directory for download, created if necessary.
%                Default: directory of calling function
%
% DESCRIPTION
% The URL will be downloaded in directory DIRNAME (to be created if
% needed). The resulting file will be uncompressed in case of a zip-, gz-
% or tar-file.
%
% The main purpose of this routine is to download missing datafiles or
% datasets from PRDATAFILES and PRDATASETS.
%
% SEE ALSO (PRTools Guide)
% PRDATASETS, PRDATAFILES

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com

function [status,dirname] = prdownload(url,dirname,siz)

% create target directory if needed
if nargin < 3, siz = []; end
if nargin > 1
  if exist(dirname,'dir') ~= 7
    success = mkdir(dirname);
    if ~success
    	maindir = input(['It is not possible to create a directory here.' newline ...
        'Please give another location: '],'s');
      dirname = fullfile(maindir,dirname);
      mkdir(dirname);
    end
  end
else
  dirname = callerdir;
end

% download
[dummy,ff,xx]= fileparts(url);
filename = [ff xx];
if isempty(siz)
  disp(['Downloading ' filename ' ....'])
else
  disp(['Downloading ' filename ' (' num2str(siz) ' MB) ....'])
end
rfilename = fullfile(dirname,filename);
if ~usejava('jvm') & isunix
  disp('isunix')
  disp(['wget -q -O ' rfilename [' '] url])
	stat = unix(['wget -q -O ' rfilename [' '] url]);
	status = (stat == 0);
else
  if verLessThan('matlab','8')
    [f,status] = urlwrite(url,rfilename);
  else
    [f,status] = urlwrite(url,rfilename,'TimeOut',20);
  end
end
if status == 0
	error('Server unreachable, timed out or file not found')
end

decompress(rfilename);

function dirname = callerdir

ss = dbstack;
if length(ss) < 3
  % no caller, commandline call
  dirname = pwd;
else
  dirname = fileparts(ss(3).name);
end

function decompress(file)

[dirname,filename,ext] = fileparts(file);
if any(strcmp(ext,{'.zip','.gz','.tgz','.tar'}))
  disp('Decompression ....')
  if strcmp(ext,'.zip')
    if ~usejava('jvm') & isunix
      [stat,s] = unix(['unzip ' file ' -d ' dirname]);
    else
      unzip(file,dirname);
    end
  elseif strcmp(ext,'.gz')
    gunzip(file,dirname);
  elseif strcmp(ext,'.tar')| strcmp(ext,'.tgz')
    untar(file,dirname);
  end
  delete(file);
  decompress(fullfile(dirname,filename));
end

author: marianux

repo: ecg-kit

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


text:
function correlationVariables = cookbook(corr)
   %% Correlation object cookbook
   % This is a short demonstration of some of the major tools in the waveform
   % correlation toolbox. The overarching goal of the toolbox is to provide a
   % set of functions for performing routine manipulations on a large set of
   % waveforms (tens to thousands) which may bare some similarity to each
   % other - cross correlation, stacking, grouping, plotting, interferometry
   % etc. Detailed usages can be found for each function with the HELP
   % command. At a minimum, read HELP CORRELATION for the basic architecture
   % of a CORRELATION object. Detailed installation and version notes can be
   % found in GISMO/@correlation/README.txt
   %
   % Data can be imported to a CORRELATION object from several sources.
   % Primary sources are intended to be Antelope and Winston databases. Since
   % trace handling is carried out with the WAVEFORM toolbox though, any data
   % that can be converted to a WAVEFORM object can be converted to a
   % CORRELATION object, including SAC files. All uses of CORRELATION require
   % that the WAVEFORM toolbox be installed.
   %
   % [NOTE that the set(gcf,'Position',...) lines used repeatedly below are
   % just for display purposes in this document.]

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   %% Use traces from the demo dataset
   % For the purposes of this cookbook example we use the built-in correlation
   % demo dataset of 100 traces drawn from a broadband sensor AU13 at
   % Augustine volcano.

   c = NewCorrelation.demo();

   %% Load traces from an outside datasource
   % This approach requires that the Antelope toolbox for Matlab be installed.
   % The waveform object must also be installed. Once a list of "trigger
   % times" has been defined (in Matlab date format), reading traces can
   % be accomplished with a single command. Loading a few thousand short
   % traces into a single correlation object is no problem. See HELP
   % DATASOURCE and HELP SCNLOBJECT for a full description of these
   % prerequisite tools.
   %
   %>> time_str = [
   %     '01/11/2006 03:05:14.363'
   %     '01/11/2006 03:06:24.336'
   %     ...
   %     '01/12/2006 06:25:59.710'
   %     '01/12/2006 06:27:23.797'
   %     ];
   %>> time_num = datenum(time_str);
   %>> ds = datasource('antelope','archive');
   %>> scnl = scnlobject('sta','chan');
   %>> c = NewCorrelation(ds,scnl,time_num,-5,10)

   %% How to use a correlation object
   % List what functions can be used on a CORRELATION object. Display the 7
   % fields contained with each CORRELATION object.

   methods(c)
   disp(c)

   %% Plot raw data
   % Use 'wig' option for a traditional wiggle plot. 'sha' displays shaded
   % traces often useful when viewing more than ~50 traces at a
   % time.

   plot(c,'wig');
   set(gcf,'Position',[50 50 700 400]);

   %% Crop, taper, filter traces
   % This particular use of crop is unnecessary but it demonstrates cropping
   % the data (in this case on the low end) and zero padding (in this case on
   % the high end.) BUTTER applies a zero-phase Butterworth filter to the
   % data. It is important to at least apply a high pass filter to data before
   % cross-correlating as long period rolls can greatly bias the correlation.

   c = c.crop(-4,12);
   c = taper(c);
   c = butter(c,[.5 10]);

   %% Cross correlate
   % Cross correlate waveforms, identifying the maximum correlation value and
   % the relative time offset between each pair of traces to acheive the
   % maximum correlation. For many hundreds of events, this can be a time
   % consuming process. XCORR provides status information along the way.

   c = xcorr(c,[0 3.5]);

   %% Similarity matrix
   % Sort traces by time and plot the correlation matrix - the set of maximum
   % cross-correlation values between each pair of traces. To access the
   % correlation matrix directly, first extract it from the correlation object
   % using the get command.

   c = sort(c);
   plot(c,'corr');
   set(gcf,'Position',[50 50 500 400]);
   corr_matrix = c.corrmatrix;
   corr_matrix(1:5,1:5)

   %% Lag matrix
   % Similar to above, this function plots the matrix of lag times (in
   % seconds) that yeild the maximum correlation between each pair of traces.

   plot(c,'lag');
   set(gcf,'Position',[50 50 500 400]);
   lag_matrix = c.lags;
   lag_matrix(1:5,1:5)

   %% Realign traces
   % Apply the time corrections in the lag matrix to the trigger times. This
   % will result in highly correlated traces being well aligned. Plot the
   % results as a shaded waveform plot.

   c = adjusttrig(c,'MIN',1);
   plot(c,'sha');
   set(gcf,'Position',[50 50 700 400]);

   %% hierarchical cluster tree
   % Plot a hierarchical cluster tree relationship (denrogram) between traces.
   % Plot reorders traces such that they correspond with the ordering of
   % traces on the cluster tree.

   c = linkage(c);
   plot(c,'den');
   close(gcf) %NOTE: why is this closed?
   set(gcf,'Position',[50 50 600 400]);

   %% Event clusters
   % Trim the cluster tree relationship into discrete clusters of events. Subset
   % just the events of the largest cluster.

   c = cluster(c,.8);
   index = find(c,'CLUST',1);
   c1 = subset(c,index);
   plot(c1,'wig');
   set(gcf,'Position',[50 50 700 400]);

   %% Coverting waveform and correlation objects
   % Inside the CORRELATION object, waveforms are stored as Seismic Traces.
   % It is possible to extract and then reinsert them. This is use for
   % carrying out trace manipulations that do not exist in the CORRELATION
   % toolbox. Be aware if reinserting waveforms that it is possible to
   % manipulate waveforms in ways that may be incompatible with the
   % correlation object, removing waveforms, for example, or altering
   % timestamps.

   T  = c1.traces;
   for n = 1:numel(T)
      T(n) = T(n).^2 .* sign(T(n).data);
   end
   c2 = c1;
   c2.traces = T;

   %% Sign traces
   % SIGN effectively removes the ampltiude information from traces. In some
   % cases the can improve correlations in the presence of large amplitude
   % transients.

   c2 = sign(c1);
   plot(c2,'wig',.3);
   set(gcf,'Position',[50 50 700 400]);

   %% Stack traces
   % Crop, normalize and stack events. The stack is appended as an additional
   % traces at the bottom. The first NORM call is to ensure that each trace is
   % given equal weight in the stack. The second NORM call is to normalize the
   % new trace stack appended to the end. STACK does not, by default, divide
   % the summation by the number of traces.
   c1 = crop(c1,-4,9);
   c1 = norm(c1);
   c1 = stack(c1);
   c1 = norm(c1);

   %% Find residual waveform
   % Subtract the stacked waveform from all other traces. The final trace is
   % zeroed out because it was subtracted from itself (the stack).

   c2 = norm(c1);
   c2 = minus(c2);
   plot(c2,'raw',.5)
   set(gcf,'Position',[50 50 700 400]);

   %% Interferogram
   % Calculate the correlation and lag values on narrow time windows within
   % each trace relative to a master event. By default, the final trace is
   % used as the master event, which in this case is the trace stack. Both the
   % maximum correlations and the lag times can be plotted behind a wiggle
   % trace plot, or they can be exported as matrices. This can be a
   % computationally intensive routine. Feedback is given to show progress
   % through the calculation. See HELP CORRELATION/PLOT for an explanation of
   % the lag interferogram color scale.

   c1 = xcorr(c1,[1 3]);       % align traces on the initial wavelet
   c1 = adjusttrig(c1);        %     "     "
   c1 = interferogram(c1);
   plot(c1,'int',1,'corr');    % plot the correlation value interferogram
   set(gcf,'Position',[50 50 700 400]);
   plot(c1,'int',1,'lag',.01);    % plot the lag value interferogram
   set(gcf,'Position',[50 50 700 400]);

   %% Occurence plot
   % Display stacked traces for the 10 largest clusters and plot them against
   % the time history of the events.

   c = crop(c,-3,10);
   plot(c,'occurence',1,1:10);
   set(gcf,'Position',[50 50 700 400]);

   %% Overlay traces
   % Plot all traces from the largest cluster together with the stack of
   % the traces.

   index = find(c,'CLUST',1);
   plot(c,'overlay',1,index);
   set(gcf,'Position',[50 50 700 200]);

   correlationVariables.c1 = c1;
   correlationVariables.c2 = c2;
   if exist('corr','var'),
   correlationVariables.corr = corr;
   end;
   correlationVariables.corr_matrix = corr_matrix;
   correlationVariables.index = index;
   correlationVariables.n = n;
   correlationVariables.T = T;

end

author: geoscience-community-codes

repo: GISMO


text:
function [L,C] = kmeans(X,k)
%KMEANS Cluster multivariate data using the k-means++ algorithm.
%   [L,C] = kmeans(X,k) produces a 1-by-size(X,2) vector L with one class
%   label per column in X and a size(X,1)-by-k matrix C containing the
%   centers corresponding to each class.

%   Version: 2013-02-08
%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%
%   References:
%   [1] J. B. MacQueen, "Some Methods for Classification and Analysis of
%       MultiVariate Observations", in Proc. of the fifth Berkeley
%       Symposium on Mathematical Statistics and Probability, L. M. L. Cam
%       and J. Neyman, eds., vol. 1, UC Press, 1967, pp. 281-297.
%   [2] D. Arthur and S. Vassilvitskii, "k-means++: The Advantages of
%       Careful Seeding", Technical Report 2006-13, Stanford InfoLab, 2006.

L = [];
L1 = 0;

while length(unique(L)) ~= k

    % The k-means++ initialization.
    C = X(:,1+round(rand*(size(X,2)-1)));
    L = ones(1,size(X,2));
    for i = 2:k
        D = X-C(:,L);
        D = cumsum(sqrt(dot(D,D,1)));
        if D(end) == 0, C(:,i:k) = X(:,ones(1,k-i+1)); return; end
        C(:,i) = X(:,find(rand < D/D(end),1));
        [~,L] = max(bsxfun(@minus,2*real(C'*X),dot(C,C,1).'));
    end

    % The k-means algorithm.
    while any(L ~= L1)
        L1 = L;
        for i = 1:k, l = L==i; C(:,i) = sum(X(:,l),2)/sum(l); end
        [~,L] = max(bsxfun(@minus,2*real(C'*X),dot(C,C,1).'),[],1);
    end

end


text:
function prtPlotUtilDataSetExploreGuiWithNavigation(ds,AdditionalOptions)
% Internal function,
% xxx Need Help xxx - see prtDataSetClass.explore

% AdditionalOptions must be structure with the field additionalOnClickFunction
% additionalOnClickFunction is a function handles that accepts 3 input
% arguments:  DataSet, clickedObsInd, currentlyDisplayedFeatures
%
% Example:
% ds = prtDataGenIris;
% O.additionalOnClickFunction = @(~,clickInd,~)disp(clickInd);
% prtPlotUtilDataSetExploreGuiWithNavigation(ds,O)
% % when you click the index clicked will be displayed in the command
% % window.

if nargin < 2 || isempty(AdditionalOptions)
    AdditionalOptions.additionalOnClickFunction = @(DataSet,clickedObsInd,currentlyDisplayedFeatures)([]); %Nothign extra
end
assert(isstruct(AdditionalOptions) && isscalar(AdditionalOptions) && isfield(AdditionalOptions,'additionalOnClickFunction'),'prt:prtDataSetClassExplore','AdditionalOptions must be a scalar structure with the field "additionalOnClickFunction".');
assert(isa(AdditionalOptions.additionalOnClickFunction,'function_handle') && nargin(AdditionalOptions.additionalOnClickFunction)==3,'prt:prtDataSetClassExplore','AdditionalOptions.additionalOnClickFunction must be a function handle that accepts two input arguments: DataSet,clickedObsInd,currentlyDisplayedFeatures.');

if ds.nFeatures > 1
    plotInds = [1 2 0];
elseif ds.nFeatures > 0
    plotInds = [1 0 0];
else
    error('prt:prtDataSetClassExplore','Dataset has zero features and cannot be explored');
end

if strcmpi(get(gcf,'tag'),'prtDataSetExplorerControl')
    % Current figure is nav controls
    % so make a new one
    figure
end

% Current figure is new (made by gcf or figure calls above) or existed
plotAxes = gca; % Will make an axes if one isn't already there.
if strcmpi(get(plotAxes,'tag'),'prtDataSetExploreAxes')
    % We are calling explore into an existing explore axes
    % This is bad. We need to destroy the existing nav control for this
    % explorer window (if it still exists).
    oldNavHandle = get(plotAxes,'UserData');
    if ishandle(oldNavHandle)
        close(oldNavHandle)
    end
end

% Make control panel figure
% Make axes current
plotAxesFig = gcf;
set(plotAxesFig,'visible','off');

legend('off');

navFigSize = [300 250];
navFigPad = [18 54];

plotAxesFigPos = get(plotAxesFig,'position');
% navFigPosTop = plotAxesFigPos(2)+plotAxesFigPos(4);
% navFigPosLeft = plotAxesFigPos(1)+plotAxesFigPos(3);

oldUnits = get(plotAxes,'units');
set(plotAxes,'Units','pixels');
plotAxesOuterPos = get(plotAxes,'outerposition');
set(plotAxes,'units',oldUnits);

navFigPosTop = plotAxesFigPos(2)+plotAxesOuterPos(2)-1+plotAxesOuterPos(4) + navFigPad(2);
navFigPosLeft = plotAxesFigPos(1)+plotAxesOuterPos(1)-1+plotAxesOuterPos(3) + navFigPad(1);

% Make sure we aren't off the screen
screenSize = get(0,'screensize');
navFigPosTop = min([navFigPosTop, screenSize(4)-50]);
navFigPosLeft = min([navFigPosLeft, screenSize(3)-navFigSize(2)-60]);

navFigPos = cat(2,navFigPosLeft, navFigPosTop-navFigSize(2), navFigSize(1), navFigSize(2));

navFigH = figure('Name','PRT Dataset Explorer Controls',...
    'Menu','none',...
    'Toolbar','none',...
    'Units','pixels',...
    'Position',navFigPos,...
    'NextPlot','new',...
    'tag','prtDataSetExplorerControl',...
    'DockControls','off',...
    'visible','off');

try
    set(navFigH,'Number','Off')
catch
    try
        set(navFigH,'NumberTitle','off');
    end
end

% Just in case anyone tries to plot in this window we will plot that inside
% an invisible axes
invisibleAxes = axes('parent',navFigH,...
    'units','pixels',...
    'position',[0 0  1 1],...
    'visible','off',...
    'handlevisibility','on'); %#ok<NASGU>

tabGroupH = createTabGroup(ds);

PlotHandles = []; % Add to all workspaces
infoClickedInd = []; % Add to all workspaces
legendHandle = []; % Add to all workspaces
legendStruct = []; % Add to all workspaces
remakePlot();

% We don't display the plots until now because it is kind of slow to make
% all of the gui elements and MATLAB will sometimes try to draw the half
% completed windows
set(plotAxesFig,'visible','on');
set(navFigH,'visible','on');

    function plotAxesDeleteFunction(myHandle, evenData) %#ok<INUSD>
        try %#ok<TRYNC>
            close(navFigH)
        end
    end
    function uitableEditFun(source, eventData) %#ok<INUSD>
        updateVisibleClasses();
    end

    function featureSelectPopupCallback(myHandle, eventData, varargin)  %#ok<INUSL>
        cVal = get(myHandle,'value');
        axisInd = varargin{1};

        if axisInd > 1
            cVal = cVal - 1;
        end
        plotInds(axisInd) = cVal;

        remakePlot();
    end

    function updateVisibleClasses()
        cData = get(tabGroupH.uitable,'data');
        displayLogical = cat(1,cData{:,2});

        if isempty(displayLogical)
            displayLogical = true;
            % Unlabeled dataset
        end

        for iClass = 1:length(PlotHandles)
            if displayLogical(iClass)
                set(PlotHandles(iClass),'visible','on');
            else
                set(PlotHandles(iClass),'visible','off');
            end
        end
        drawnow;
    end

    function remakePlot()

        actualPlotDims = plotInds(plotInds>=1);
        axes(plotAxes); %#ok<MAXES>

        % Before we plot we get the legend properties (it may not exist)
        legendHandle = legend(plotAxes);
        legendStruct = get(legendHandle);

        PlotHandles = plot(ds,actualPlotDims);
        % Set the plotAxes delete function so that we delete the controls
        % when necessary.
        % We have to reset this each time
        setAxesProperties();

        updateVisibleClasses();

        % Make sure the navFig is on top
        figure(navFigH);
    end

    function setAxesProperties()
        set(plotAxes,'deleteFcn',@plotAxesDeleteFunction)
        set(plotAxes,'tag','prtDataSetExploreAxes')
        set(plotAxes,'UserData',navFigH);
        set(plotAxes,'ButtonDownFcn',@plotAxesOnClick);
        set(PlotHandles,'HitTest','off');

        % I know the function is called setAxesProperties but this is also
        % a great place to put the legend back where we want, amoung other
        % things.
        legendHandle = legend(plotAxes);
        if ~isempty(legendStruct)
            % legendStruct has all of the properties of the old legend
            % we will reset some of the new legend's properties to match
            % the old one. This ends the annoyances of moving legends and
            % having them move back to their default location.
            set(legendHandle,'Position',legendStruct.Position,...
                             'Location',legendStruct.Location,...
                             'Orientation',legendStruct.Orientation,...
                             'String',legendStruct.String);
        end
    end

    function plotAxesOnClick(myHandle,eventData)  %#ok<INUSD>
        actualPlotDims = plotInds(plotInds>=1);
        if ~ishandle(navFigH)
            return
        end

        cData = get(tabGroupH.uitable,'data');
        displayLogical = cat(1,cData{:,2});

        if all(displayLogical)
            data = ds.getFeatures(actualPlotDims);

            if size(data,2) == 1
                data = cat(2,data,ds.getTargets());
            end

        else
            actualInds = ismember(ds.getTargetsClassInd, find(displayLogical));
            data = ds.getObservations(actualInds,actualPlotDims);

            if size(data,2) == 1
                data = cat(2,data,ds.getTargets(actualInds));
            end

        end

        [rP,rD] = rotateDataAndClick(data);

        dist = prtDistanceEuclidean(rP,rD);
        [minDist,clickedObsInd] = min(dist);  %#ok<ASGLU>

        if ~all(displayLogical)
            actualInds = find(actualInds);
            clickedObsInd = actualInds(clickedObsInd);
        end

        displayInfo(clickedObsInd);

        AdditionalOptions.additionalOnClickFunction(ds,clickedObsInd,plotInds);

    end

    function displayInfo(clickedObsInd)
        obsName = ds.getObservationNames(clickedObsInd);

        if ds.isLabeled
            cClassInd = ds.getTargetsClassInd(clickedObsInd);
        else
            cClassInd = [];
        end

        if isempty(cClassInd)
            cString = sprintf('Closest Observation:\n\t Name: %s',obsName{1});
        else
            className = ds.getClassNamesByClassInd(cClassInd);

            if ~iscell(className)
                className = {className};
            end
            cString = sprintf('Closest Observation:\n\t Name: %s\n\tClass: %s',obsName{1},className{1});
        end

        if length(ds.observationInfo) > clickedObsInd
            infoClickedInd = clickedObsInd;
            set(tabGroupH.infoButton,'enable','on');
            if ~isdeployed % can't use evalc if deployed
                s = ds.observationInfo(infoClickedInd);
                infoStr = evalc('display(s)');
                infoStr = infoStr(6:end); % 6 here is for the 's =   '

                cString = cat(2, cString, sprintf('\n\nObservation Info:'),infoStr);
            end
        else
            set(tabGroupH.infoButton,'enable','off');
        end

        set(tabGroupH.infoText,'string',cString);

    end

    function infoToWorkspaceCallback(myHandle,eventData)  %#ok<INUSD>

        c = ds.observationInfo(infoClickedInd);

        assignin('base','prtPlotUtilDataSetExploreGuiWithNavigationTempVar',c);
        openvar('prtPlotUtilDataSetExploreGuiWithNavigationTempVar');
    end

    function [rotatedData, rotatedClick] = rotateDataAndClick(data)
        % Used internally; from Click3dPoint from matlab central;
        % See prtExternal.ClickA3DPoint.()
        %
        % The copyright info from that file is below.
        %
        % Copyright (c) 2009, Babak Taati 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.

        point = get(plotAxes, 'CurrentPoint'); % mouse click position
        camPos = get(plotAxes, 'CameraPosition'); % camera position
        camTgt = get(plotAxes, 'CameraTarget'); % where the camera is pointing to

        camDir = camPos - camTgt; % camera direction
        camUpVect = get(plotAxes, 'CameraUpVector'); % camera 'up' vector

        % build an orthonormal frame based on the viewing direction and the
        % up vector (the "view frame")
        zAxis = camDir/norm(camDir);
        upAxis = camUpVect/norm(camUpVect);
        xAxis = cross(upAxis, zAxis);
        yAxis = cross(zAxis, xAxis);

        rot = [xAxis; yAxis; zAxis]; % view rotation

        % the clicked point represented in the view frame
        rotatedClick = rot * point' ;
        rotatedClick = rotatedClick(1:2,1)';

        if size(data,2) < 3
            if size(data,2) < 2
                rot = rot(1);
            else
                rot = rot(1:2,1:2);
            end
        end

        % the point cloud represented in the view frame
        rotatedData = (rot * data')';

        if size(rotatedData,2) > 2
            rotatedData = rotatedData(:,1:2);
        end

        if size(data,2) < 2
            rotatedClick = rotatedClick(1);
        end
    end

    function H = createTabGroup(ds)
        % Make control panel uitabs and uipanels
        H.navTabGroup = prtUtilUitabgroup('parent',navFigH);
        H.navTab(1) = prtUtilUitab(H.navTabGroup, 'title', 'Features');
        H.navTab(2) = prtUtilUitab(H.navTabGroup, 'title', 'Classes');
        H.navTab(3) = prtUtilUitab(H.navTabGroup, 'title', 'Info');

        H.navPanel(1) = uipanel(H.navTab(1));
        H.navPanel(2) = uipanel(H.navTab(2));
        H.navPanel(3) = uipanel(H.navTab(3));

        % Make Feature Selection panel
        featureNames = ds.getFeatureNames;
        featureNames = featureNames(:);
        H.xAxesLabel = uicontrol(H.navTab(1),'style','text',...
            'units','normalized',...
            'position',[0.025 0.875 0.95 0.1],...
            'string','X Axis',...
            'FontWeight','bold',...
            'horizontalAlignment','left',...
            'FontUnits','normalized',...
            'FontSize',0.75);

        H.xPopUp = uicontrol(H.navTab(1),'style','popup',...
            'units','normalized',...
            'position',[0.025 0.675 0.95 0.2],...
            'string',featureNames,...
            'FontUnits','normalized',...
            'FontSize',0.5,...
            'callback',{@featureSelectPopupCallback 1});

        H.yAxesLabel = uicontrol(H.navTab(1),'style','text',...
            'units','normalized',...
            'position',[0.025 0.55 0.95 0.1],...
            'string','Y Axis',...
            'FontWeight','bold',...
            'horizontalAlignment','left',...
            'FontUnits','normalized',...
            'FontSize',0.75);

        H.yPopUp = uicontrol(H.navTab(1),'style','popup',...
            'units','normalized',...
            'position',[0.025 0.325 0.95 0.225],...
            'string',cat(1,{'None'}, featureNames),...
            'FontUnits','normalized',...
            'FontSize',0.5,...
            'callback',{@featureSelectPopupCallback 2});

        H.zAxesLabel = uicontrol(H.navTab(1),'style','text',...
            'units','normalized',...
            'position',[0.025 0.225 0.95 0.1],...
            'string','Z Axis',...
            'FontWeight','bold',...
            'horizontalAlignment','left',...
            'FontUnits','normalized',...
            'FontSize',0.75);

        H.zPopUp = uicontrol(H.navTab(1),'style','popup',...
            'units','normalized',...
            'position',[0.025 0.025 0.95 0.2],...
            'string',cat(1,{'None'}, featureNames),...
            'FontUnits','normalized',...
            'FontSize',0.5,...
            'callback',{@featureSelectPopupCallback 3});

        set(H.xPopUp,'value',plotInds(1)) % Always have at least 1 features
        if plotInds(2) > 0
            set(H.yPopUp,'value',plotInds(2)+1)
        else
            set(H.yPopUp,'value',1) % None;
        end
        if plotInds(3) > 0
            set(H.zPopUp,'value',plotInds(3)+1)
        else
            set(H.zPopUp,'value',1) % None;
        end

        if ds.hasUnlabeled
            tableData = cat(2,cat(1,ds.getClassNames,{prtPlotUtilUnlabeledLegendString}),num2cell(true(ds.nClasses+1,1)));
        else
            tableData = cat(2,ds.getClassNames,num2cell(true(ds.nClasses,1)));
        end

        % Make classes uitable with tick boxes
        H.uitable = uitable('parent',H.navPanel(2),...
            'units','Normalized',...
            'position',[0.025 0.025 0.95 0.95],...
            'columnFormat',{'char','logical'},...
            'ColumnEditable',[false true],...
            'RowName',[],...
            'ColumnName',{'Class' 'Plot'},...
            'FontUnits','normalized',...
            'FontSize',0.1,...
            'SelectionHighlight','off',...
            'CellEditCallback',@uitableEditFun,...
            'data',tableData);

        % Column widths must be specified in pixels
        set(H.uitable,'units','pixels')
        uitablePos = floor(get(H.uitable,'position'));
        leftColumnWidth = floor(uitablePos(3)*0.7);
        set(H.uitable,'ColumnWidth', {leftColumnWidth uitablePos(3)-leftColumnWidth-22});
        set(H.uitable,'units','normalized')

        H.infoText = uicontrol('style','edit',...
                                'Max',200,'enable','inactive',...
                               'parent',H.navPanel(3),...
                               'units','normalized',...
                               'position',[0.025 0.125 0.95 0.85],...
                               'FontSize',10,...
                               'FontName',get(0,'FixedWidthFontName'),...
                               'HorizontalAlignment','Left',...
                               'string','Click in the axes to inspect observations.');

        H.infoButton = uicontrol('parent',H.navPanel(3),...
                                 'string','View Obs. Info',...
                                 'units','normalized',...
                                 'position',[0.2 0.025 0.6 0.095],...
                                 'callback',@infoToWorkspaceCallback,...
                                 'enable','off');

    end

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [F] = sTopMain(fileToRead)
%STOPMAIN Plots an animated spinning top with Cardan mounting from raw
%animation data.
%
%   example call: sTopMain('animation_A.dat')
%
%   See PDF documentation for details.
%
%   Alexander Erlich (alexander.erlich@gmail.com)

rawData = importdata(fileToRead);
time = rawData(:,1);
phi = rawData(:,2);
theta = -1.*rawData(:,3);
psi = rawData(:,4);

sTopEuler;  % creating surf plots for the spinning Top in its initial position
            % (phi, theta, psi) = (0, 0, 0)

for i=1:length(phi)-1
    if (ishandle(H1)~= 0)

        % computing the new position
        [Xtemp1,Ytemp1,Ztemp1]=...
            multiplyEuMat(EuMat(phi(i),theta(i),psi(i)),X1,Y1,Z1);
        [Xtemp2,Ytemp2,Ztemp2]=...
            multiplyEuMat(EuMat(phi(i),theta(i),0),X2,Y2,Z2);
        [Xtemp3,Ytemp3,Ztemp3]=...
            multiplyEuMat(EuMat(phi(i),0,0),X3,Y3,Z3);
        [Xtempbar1,Ytempbar1,Ztempbar1]=...
            multiplyEuMat(EuMat(phi(i),theta(i),psi(i)),Xbar1,Ybar1,Zbar1);
        [Xtempbar2,Ytempbar2,Ztempbar2]=...
            multiplyEuMat(EuMat(phi(i),theta(i),psi(i)),Xbar2,Ybar2,Zbar2);

        % updating the new position
        set(H1,'XData',Xtemp1,'YData',Ytemp1,'ZData',Ztemp1);
        set(H2,'XData',Xtemp2,'YData',Ytemp2,'ZData',Ztemp2);
        set(H3,'XData',Xtemp3,'YData',Ytemp3,'ZData',Ztemp3);
        set(bar1,'XData',Xtempbar1,'YData',Ytempbar1,'ZData',Ztempbar1);
        set(bar2,'XData',Xtempbar2,'YData',Ytempbar2,'ZData',Ztempbar2);

        drawnow % drawing the new position
        pause(time(i+1)-time(i)); % pausing until next frame is shown
    end
end


text:
function [x, fx, exitFlag] = bisection(f,lb,ub,target,options)
% BISECTION Fast and robust root-finding method that handles n-dim arrays.
%
%   [x,fVal,ExitFlag] = BISECTION(f,LB,UB,target,options) finds x +/- TolX
%   (LB < x < UB) such that f(x) = target +/- TolFun.
%
%   x = BISECTION(f,LB,UB) finds the root(s) of function f on the interval
%   [LB, UB], i.e. finds x such that f(x) = 0 where LB <= x <= UB. f will
%   never be evaluated outside of the interval specified by LB and UB. f
%   should have only one root and f(UB) and f(LB) must bound it. Elements
%   of x are NaN for instances where a solution could not be found.
%
%   x = BISECTION(f,LB,UB,target) finds x such that f(x) = target.
%
%   x = BISECTION(f,LB,UB,target,TolX) will terminate the search when the
%   search interval is smaller than TolX (TolX must be positive).
%
%   x = BISECTION(f,LB,UB,target,options) solves with the default
%   parameters replaced by values in the structure OPTIONS, an argument
%   created with the OPTIMSET function. Used options are TolX and TolFun.
%   Note that OPTIMSET will not allow arrays for tolerances, so set the
%   fields of the options structure manually for non-scalar TolX or TolFun.
%
%   [x,fVal] = BISECTION(f,...) returns the value of f evaluated at x.
%
%   [x,fVal,ExitFlag] = BISECTION(...) returns an ExitFlag that describes
%   the exit condition of BISECTION. Possible values of elements of
%   ExitFlag and the corresponding exit conditions are
%
%       1   Search interval smaller than TolX.
%       2   Function value within TolFun of target.
%       3   Search interval smaller than TolX AND function value within
%           TolFun of target.
%      -1   No solution found.
%
%   Any or all of f(scalar), f(array), LB, UB, target, TolX, or TolFun may
%   be scalar or n-dim arrays. All non-scalar arrays must be the same size.
%   All outputs will be this size.
%
%   Default values are target = 0, TolX = 1e-6, and TolFun = 0.
%
%   There is no iteration limit. This is because BISECTION (with a TolX
%   that won't introduce numerical issues) is guaranteed to converge if f
%   is a continuous function on the interval [UB, LB] and f(x)-target
%   changes sign on the interval.
%
%   The bisection method is very robust root-finding method. The absolute
%   error is halved at each step so the method converges linearly. However,
%   Brent's method (such as implemented in FZERO) can converge
%   superlinearly and is as robust. FZERO also has more features and input
%   checking, so use BISECTION in cases where either the optimization
%   toolbox is unavailable or if FZERO would have to be implemented in a
%   loop to solve multiple cases, in which case BISECTION will be much
%   faster because of vectorization.
%
%   Define LB, UB, target, TolX, and TolFun for each specific application
%   using great care for the following reasons:
%     - There is no iteration limit, so given an unsolvable task, such as
%       TolX = TolFun = 0, BISECTION remains in an unending loop.
%     - Spacing between very large floating point numbers is likely to be
%       greater than TolX.
%     - There is no initial check to make sure that f(x) - target changes
%       sign between LB and UB.
%     - Very large or very small numbers can introduce numerical issues.
%
%   Example 1: Find cube root of array 'target' without using NTHROOT and
%   compare speed to using FZERO.
%       options = optimset('TolX', 1e-9);
%       target = [(-100:.1:100)' (-1000:1:1000)'];
%
%       tic;
%       xfz = zeros(size(target));
%       for ii = 1:numel(target)
%           xfz(ii) = fzero(@(x) x.^3-target(ii), [-20 20], options);
%       end
%       fzero_time = toc
%
%       tic;
%       xbis = bisection(@(x) x.^3, -20, 20, target, options);
%       bisection_time = toc
%
%       fprintf('FZERO took %0.0f times longer than BISECTION.\n',...
%                   fzero_time/bisection_time)
%
%   Example 2: Find roots by varying the function coefficients.
%       [A, B] = meshgrid(linspace(1,2,6), linspace(4,12,10));
%       f = @(x) A.*x.^0.2 + B.*x.^0.87 - 15;
%       xstar = bisection(f,0,5);
%
%   See also FZERO, FMINBND, OPTIMSET, FUNCTION_HANDLE.
%
%   [x,fVal,ExitFlag] = BISECTION(f,LB,UB,target,options)

%   Copyright 2010-2013 Sky Sartorius
%   Author  - Sky Sartorius
%   Contact - www.mathworks.com/matlabcentral/fileexchange/authors/101715

% --- Process inputs. ---
% Set default values
tolX    = 1e-6;
tolFun  = 0;
if nargin == 5
    if isstruct(options)
        if isfield(options,'TolX') && ~isempty(options.TolX)
            tolX = options.TolX;
        end
        if isfield(options,'TolFun') && ~isempty(options.TolFun)
            tolFun = options.TolFun;
        end
    else
        tolX = options;
    end
end
if nargin<4 || isempty(target); target=0; end

ub_in = ub; lb_in = lb;
f = @(x) f(x) - target;

% --- Flip UB and LB if necessary. ---
isFlipped = lb>ub;
if any(isFlipped(:))
    ub(isFlipped) = lb_in(isFlipped);
    lb(isFlipped) = ub_in(isFlipped);
    ub_in = ub; lb_in = lb;
end

% --- Make sure everything is the same size for a non-scalar problem. ---
if isscalar(lb) && isscalar(ub)
    % Test if f returns multiple outputs for scalar input.
    if ~isscalar(target)
        ub = ub + zeros(size(target));
    else
        jnk = f(ub);
        if ~isscalar(jnk)
            ub = ub + zeros(size(jnk));
        end
    end
end

% Check if lb and/or ub need to be made into arrays.
if isscalar(lb) && ~isscalar(ub)
    lb = lb + zeros(size(ub));
elseif ~isscalar(lb) && isscalar(ub)
    ub = ub + zeros(size(lb));
end

testconvergence
% --- Iterate ---
while any(stillNotDone(:))
    bigger  = fx.*f(ub) > 0;
    ub(bigger)= x(bigger);
    lb(~bigger)= x(~bigger);

    testconvergence;
end

    function testconvergence
        x=(ub+lb)/2;
        fx=f(x);
        outsideTolFun =  abs(fx)  > tolFun;
        outsideTolX =   (ub - lb) > tolX;
        stillNotDone = outsideTolX & outsideTolFun;
    end

% --- Check that f(x+tolX) and f(x-tolX) have opposite sign. ---
fu = f(min(x+tolX,ub_in));
fl = f(max(x-tolX,lb_in));
unboundedRoot = (fu.*fl) > 0;

% Throw out unbounded results if not meeting TolFun convergence criteria.
x(unboundedRoot & outsideTolFun) = NaN;

% --- Catch NaN elements of UB, LB, target, or other funky stuff. ---
x(isnan(fx)) = NaN;

% --- Characterize results. ---
fx = fx + target;
if nargout > 2
    exitFlag                                    = +~outsideTolX;
    exitFlag(~outsideTolFun)                    =  2;
    exitFlag(~outsideTolFun & ~outsideTolX)     =  3;
    exitFlag(isnan(x))                          = -1;
end

end

% V2: July     2010
% V3: December 2012
% don't remember when
%   typo line 39; added fn handle to see also; made array in example 2
%   smaller; changed wording in example 1
% 2013-08-23
%  -changed scalar*ones(...) calls to scalar+zeros(...) calls based on
%   http://undocumentedmatlab.com/blog/allocation-performance-take-2/
%  -rearranged help block and formatted a tiny bit


text:
function [Rob,Sen,Lmk,Obs,Frm,Fac] = solveGraphSchur(Rob,Sen,Lmk,Obs,Frm,Fac,options)

% SOLVEGRAPHSCHUR Solves the SLAM graph using Schur decomposition.
%   [Rob,Sen,Lmk,Obs,Frm,Fac] = solveGraphSchur(Rob,Sen,Lmk,Obs,Frm,Fac)
%   solves the graph-SLAM problem using Schur decomposition of the
%   Hessian matrix.
%
%   IMPORTANT NOTE: This method is illustrative and constitutes the
%   motivation for this toolbox. One can achieve better performances, both
%   in computing time and possibly in robustness and accuracy, by using
%   Matlab's built-in nonlinear optimization tools, such as LSQNONLIN.
%
%   See courseSLAM.pdf in the documentation for details about the Schur
%   decomposition for solving the graph-SLAM problem.
%
%   See also ERRORSTATEJACOBIANS, UPDATESTATES, COMPUTEERROR,
%   COMPUTERESIDUAL, COLAMD, '\', MLDIVIDE, LSQNONLIN.

% Copyright 2015-    Joan Sola @ IRI-UPC-CSIC.

global Map

% Control of iterations and exit conditions
n_iter      = options.niterations; % exit criterion of number of iterations
target_dres = options.target_dres; % exit criterion for error variation
target_res  = options.target_res;  % exit criterion for current residual
res_old     = 1e10;                   % last iteration's error

% Map states range
Map.mr = find(Map.used);

for it = 1:n_iter

    % Compute Jacobians for projection onto the manifold
    [Frm,Lmk] = errorStateJacobians(Frm,Lmk);

    % Build Hessian H and rhs vector b, in global Map
    Fac       = buildProblem(Rob,Sen,Lmk,Obs,Frm,Fac);

    % Get partition ranges
    fr = ranges(Frm);
    lr = ranges(Lmk);

    % Get Schur complement -- The schur complement Spp is the sqrt factor
    [Map.sSff, Map.iHll] = schurc(Map.H(Map.mr,Map.mr),lr,fr,1);

    % Solve frames subproblem
    bf        = Map.b(fr) - Map.H(fr,lr) * Map.iHll * Map.b(lr);
    y         = - Map.sSff'\bf;
    Map.x(fr) = Map.sSff\y;

    % Solve landmarks subproblem
    bl        = Map.b(lr) + Map.H(lr,fr) * Map.x(fr);
    Map.x(lr) = - Map.iHll * bl;

    % Update nominal states
    [Rob,Lmk,Frm]  = updateStates(Rob,Lmk,Frm);

    % Check resulting errors
    [res, err_max] = computeResidual(Rob,Sen,Lmk,Obs,Frm,Fac);
    dres           = res - res_old;
    res_old        = res;

    if ( ( -dres <= target_dres ) || (err_max <= target_res) ) %&& ( abs(derr) < target_derr) )
        break;
    end

end

end

function Fac = buildProblem(Rob,Sen,Lmk,Obs,Frm,Fac)

% BUILDPROBLEM Build least squares problem's matrix H and vector b
%   Fac = BUILDPROBLEM(Rob,Sen,Lmk,Obs,Frm,Fac) Builds the least squares
%   problem's matrix H and vector b for a solution using sparse Schur
%   factorization of H.

global Map

% Reset Hessian and rhs vector
Map.H(Map.mr,Map.mr) = 0;
Map.b(Map.mr)    = 0;

% Iterate all factors
for fac = find([Fac.used])

    % Extract some pointers
    rob    = Fac(fac).rob;
    sen    = Fac(fac).sen;
    lmk    = Fac(fac).lmk;
    frames = Fac(fac).frames;

    % Compute factor error, info mat, and Jacobians
    [Fac(fac), e, W, ~, J1, J2, r1, r2] = computeError(...
        Rob(rob),       ...
        Sen(sen),       ...
        Lmk(lmk),       ...
        Obs(sen,lmk),   ...
        Frm(frames),    ...
        Fac(fac));

    % Compute sparse Hessian blocks
    H_11 = J1' * W * J1;
    H_12 = J1' * W * J2;
    H_22 = J2' * W * J2;

    % Compute rhs vector blocks
    b1   = J1' * W * e;
    b2   = J2' * W * e;

    % Update H and b
    Map.H(r1,r1) = Map.H(r1,r1) + H_11;
    Map.H(r1,r2) = Map.H(r1,r2) + H_12;
    Map.H(r2,r1) = Map.H(r2,r1) + H_12';
    Map.H(r2,r2) = Map.H(r2,r2) + H_22;

    Map.b(r1,1)  = Map.b(r1,1)  + b1;
    Map.b(r2,1)  = Map.b(r2,1)  + b2;

end

end

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function plot(rsam_vector, varargin)
    % RSAM/PLOT plot rsam data
    % handle = plot(rsam_vector, varargin)
    % Properties include:
    %   yaxistype, h, addgrid, addlegend, fillbelow, plotspikes, plottransients, plottremor
    % to change where the legend plots set the global variable legend_ypos
    % a positive value will be within the axes, a negative value will be below
    % default is -0.2. For within the axes, log(20) is a reasonable value.
    % yaxistype is like 'logarithmic' or 'linear'
    % h is an axes handle (or an array of axes handles)
    % use h = generatePanelHandles(numgraphs)

    % Glenn Thompson 1998-2009
    %
    % % GTHO 2009/10/26 Changed marker size from 5.0 to 1.0
    % % GTHO 2009/10/26 Changed legend position to -0.2

    % parse variable input arguments
    p = inputParser;
    p.addParameter('addgrid',false);
    p.addParameter('addlegend', true);
    p.addParameter('fillbelow', false);
    p.addParameter('yaxistype','linear'); % or log
    p.addParameter('symbol','-'); % or log
    p.addParameter('h', []);
    p.parse(varargin{:});
    addgrid = p.Results.addgrid;
    addlegend = p.Results.addlegend;
    fillbelow = p.Results.fillbelow;
    yaxistype = p.Results.yaxistype;
    symbol = p.Results.symbol;
    h = p.Results.h;
    legend_ypos = -0.2;

    % colours to plot each station
    lineColour={[0 0 0]; [0 0 1]; [1 0 0]; [0 1 0]; [.4 .4 0]; [0 .4 0 ]; [.4 0 0]; [0 0 .4]; [0.5 0.5 0.5]; [0.25 .25 .25]};

    % units - so that we put different ones on different figures
    units = {rsam_vector.units};
    unique_units = unique(units);
    %previousfignum = get(gcf,'Number');
    previousfignum = get_highest_figure_number();

    % Plot the data graphs
    for c = 1:length(rsam_vector)
        self = rsam_vector(c);
%         thisfignum = find(ismember(unique_units, self.units)) + previousfignum;
%         figure(thisfignum);
        try
            figure(previousfignum + c);
        catch
            figure;
        end

        hold on;
        t = self.dnum;
        y = self.data;

        debug.print_debug(10,sprintf('Data length: %d',length(y)));

        if length(y)>0
            %figure
            %if strcmp(rsam_vector(c).units, 'Hz')

            if strcmp(yaxistype,'linear')

                % plot on a linear axis, with station name as a y label
                % datetick too, add measure as title, fiddle with the YTick's and add max(y) in top left corner
                if ~p.Results.fillbelow
                    %handlePlot = plot(t, y, symbol, 'Color', lineColour{c});
                    handlePlot = plot(t, y, symbol, 'Color', lineColour{1});
                else
                    %handlePlot = fill([min(t) t max(t)], [min([y 0]) y min([y 0])], lineColour{c});
                    handlePlot = fill([min(t) t max(t)], [min([y 0]) y min([y 0])], lineColour{1});
                end

                % if c ~= numel(rsam_vector)
                %     set(gca,'XTickLabel','');
                % end

                % yt=get(gca,'YTick');
                % ytinterval = (yt(2)-yt(1))/2;
                % yt = yt(1) + ytinterval: ytinterval: yt(end);
                % ytl = yt';
                % ylim = get(gca, 'YLim');
                % set(gca, 'YLim', [0 ylim(2)],'YTick',yt);
                % %ylabelstr = sprintf('%s.%s %s (%s)', self.sta, self.chan, self.measure, self.units);
    %             ylabelstr = sprintf('%s', self.sta);
    %             ylabel(ylabelstr);
                ylabel(sprintf('%s',self.units))
    %             datetick('x','keeplimits');
                a = axis;
                datetick('x')
                set(gca,'XLim',[a(1) a(2)]);
                xlabel('Date/Time');
                if addlegend
                    legend();
                end

            else

                % make a logarithmic plot, with a marker size and add the station name below the x-axis like a legend
                y = log10(y);  % use log plots

                %                 handlePlot = plot(t, y, symbol, 'Color', lineColour{c},...
                %                    'MarkerSize', 1.0);
                handlePlot = plot(t, y, symbol, 'Color', lineColour{1},...
                   'MarkerSize', 1.0);
                if strfind(self.measure, 'dr')
                    %ylabel(sprintf('%s (cm^2)',self(c).measure));
                    %ylabel(sprintf('D_R (cm^2)',self(c).measure));
                    Yticks = [0.01 0.02 0.05 0.1 0.2 0.5 1 2 5 10 20 50 ];
                    Ytickmarks = log10(Yticks);
                    for count = 1:length(Yticks)
                        Yticklabels{count}=num2str(Yticks(count),3);
                    end
                    set(gca, 'YLim', [min(Ytickmarks) max(Ytickmarks)],...
                       'YTick',Ytickmarks,'YTickLabel',Yticklabels);
                end
                axis tight
                a = axis;
                datetick('x')
                set(gca,'XLim',[a(1) a(2)]);
    %
                xlabel(sprintf('Date/Time starting at %s',datestr(self.snum)))
                ylabel(sprintf('log(%s)',self.units))

            end

            if p.Results.addgrid
                grid on;
            end
            if p.Results.addlegend && ~isempty(y)
                xlim = get(gca, 'XLim');
                legend_ypos = 0.9;
                legend_xpos = c/10;
            end

            datetick('x','keeplimits')
            if ~strcmp(self.ChannelTag.string(), '...')
                tstr = sprintf('%s %s\n%s %.0f s',self.ChannelTag.string(), datestr(t(1)), self.measure, round(self.sampling_interval) );
            else
                tstr = sprintf('%s\n%s %.0f s',self.files.file, datestr(a(1)), round(self.sampling_interval) );
            end
            title(tstr)
        end
    end
end


text:
function [deg,min,sec]=deg2DegMinSec(deg)
%%DEG2DEGMINSEC Convert angular degrees to a representation in integer
%               degrees, arcminutes, and decimal arcseconds. This can be
%               useful for displaying angular quantities in a user-frieldy
%               manner.
%
%INPUTS: deg A scalar or matrix of angles degrees.
%
%OUTPUTS: deg The integer part of the input deg.
%         min The integer number of arcminutes in deg.
%         sec A decimal number of arcseconds.
%
%If a negative number of degrees is entered, then degrees, arcminutes and
%arcseconds will all come back negative.
%
%There are 60 arcminutes per degree and 60 arcseconds per arcminute.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    fracPart=deg-fix(deg);%Get rid of integer degrees.
    deg=fix(deg);
    min=fix(fracPart*60);
    fracPart=fracPart-(min/60);%Get rid of integer minutes.
    sec=60*60*fracPart;
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  img = makeMontageIfiles(baseFileName, imDim, imList, numDigits, byteFlag)
% img = makeMontageIfiles(baseFileName, imDim, imList, [numDigits], [byteFlag])
%
% Reads in reconstructed Pfile images and compiles
% them into an image montage.
%
% numDigits is the number of digits in the numerical
% filename extension. Default is 3 (eg. I.001).
%
% EXAMPLE:
%
% if your I-files are in a directory '/red/u4/mri/someIfiles' and
% have the standard 'I.001, I.002, I.003, ...' names, and they are
% 256x256 images, and there are 100 of them, then try this:
%
%     img = makeMontageIfiles('/red/u4/mri/someIfiles/I', [256 256], [1:100])
%
% To view the resulting image, use something like:
%
%    imagesc(img); colormap(gray); axis image;
%
%
% see also makeCubeIfiles

if(~exist('numDigits') | isempty(numDigits))
    numDigits = 3;
end
if(~exist('byteFlag')) byteFlag = []; end
im3d = makeCubeIfiles(baseFileName, imDim, imList, numDigits, byteFlag);
img = makeMontage(im3d);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft