text:
%% Copyright (C) 2014-2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym isnan (@var{x})
%% Return true if a symbolic expression is Not-a-Number.
%%
%% Example:
%% @example
%% @group
%% A = [sym(1) sym(0)/0 sym(1)/0; sym(nan) 1 2]
%%   @result{} A = (sym 2×3 matrix)
%%       ⎡ 1   nan  zoo⎤
%%       ⎢             ⎥
%%       ⎣nan   1    2 ⎦
%% isnan(A)
%%   @result{} ans =
%%       0   1   0
%%       1   0   0
%% @end group
%% @end example
%%
%% Note that the return is of type logical.
%%
%% @seealso{@@sym/isinf, @@sym/double}
%% @end defmethod

function r = isnan(x)

  if (nargin ~= 1)
    print_usage ();
  end

  r = uniop_bool_helper(x, 'lambda a: a is sp.nan');

end

%!shared x,zoo,oo,snan
%! oo = sym(inf);
%! zoo = sym('zoo');
%! x = sym('x');
%! snan = sym(nan);

%!test
%! % various ops that give nan
%! assert (isnan(0*oo))
%! assert (isnan(0*zoo))
%! assert (isnan(snan))
%! assert (isnan(snan-snan))
%! assert (isnan(oo+snan))
%! assert (isnan(oo-oo))
%! assert (isnan(oo-zoo))
%! assert (isnan(oo+zoo))
%! assert (~isnan(oo))
%! assert (~isnan(zoo))
%! assert (~isnan(oo+oo))

%!test
%! % more ops give nan
%! assert(isnan(x+snan))
%! assert(isnan(x*snan))
%! assert(isnan(0*snan))
%! assert(isnan(x+nan))
%! assert(isnan(x*nan))
%! assert(isnan(sym(0)*nan))

%!test
%! % array
%! assert (isequal(  isnan([oo zoo]),    [0 0]  ))
%! assert (isequal(  isnan([10 snan]),   [0 1]  ))
%! assert (isequal(  isnan([snan snan]), [1 1]  ))
%! assert (isequal(  isnan([snan x]),    [1 0]  ))

%!test
%! % sub in to algebraic expression gives nan
%! y = x - oo;
%! y = subs(y, x, oo);
%! assert(isnan(y))

%!test
%! % Must not contain string 'symbol'; these all should make an
%! % actual nan.  Actually a ctor test, not isnan.
%! y = sym(nan);
%! assert (isempty (strfind (sympy (y), 'Symbol')))
%! y = sym('nan');
%! assert (isempty (strfind (sympy (y), 'Symbol')))
%! y = sym('NaN');
%! assert (isempty( strfind (sympy (y), 'Symbol')))

repo: octsympy


text:
function mags = gtmmag(net, latent_data)
%GTMMAG	Magnification factors for a GTM
%
%	Description
%	 MAGS = GTMMAG(NET, LATENTDATA) takes a GTM structure NET, and
%	computes the magnification factors for each point the latent space
%	contained in LATENTDATA.
%
%	See also
%	GTM, GTMPOST, GTMLMEAN
%

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

errstring = consist(net, 'gtm');
if ~isempty(errstring)
  error(errstring);
end

Jacs = rbfjacob(net.rbfnet, latent_data);
nlatent = size(latent_data, 1);
mags = zeros(nlatent, 1);
temp = zeros(net.rbfnet.nin, net.rbfnet.nout);
for m = 1:nlatent
  temp = squeeze(Jacs(m, :, :));  % Turn into a 2d matrix
  mags(m) = sqrt(det(temp*temp'));
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function y = cvx_subs2str( x, mask, fieldalt )
if nargin < 2,
    mask = [ 1, 1, 1 ];
end
if nargin < 3,
    fieldalt = 0;
end
if ~isstruct( x ),
    error( 'First arugment must be a structure.' );
end
y = '';
needfield = fieldalt ~= 0;
for k = 1 : length( x ),
    try
        tp = x(k).type;
    catch %#ok
        error( 'Invalid subscript structure: field "type" is missing.' );
    end
    if ~ischar( tp ) || size( tp, 1 ) ~= 1,
        error( 'Invalid subscript entry #%d: field "type" must be a string', k );
    end
    try
        sb = x(k).subs;
    catch %#ok
        error( 'Invalid subscript structure: field "subs" is missing.' );
    end
    switch tp,
        case '.',
            if ~ischar( sb ) || size( sb, 1 ) ~= 1,
                error( 'Invalid subscript entry #%d: field name must be a string.', k );
            elseif ~isvarname( sb ),
                error( 'Invalid subscript entry #%d: invalid field name: %s', k, sb );
            end
            y = [ y, '.', sb ]; %#ok
            needfield = 0;
        case { '()', '{}' },
            if ~mask( 2 )&& tp(1) == '(',
                error( 'Invalid subscript entry #%d: array subscripts not allowed here.', k );
            elseif ~mask( 3 ) && tp(1) == '{',
                error( 'Invalid subscript entry #%d: cell subscripts not allowed here.', k );
            elseif needfield,
                error( 'Invalid subscript entry #%d: structure field expected here.', k );
            elseif ~iscell( sb ),
                error( 'Invalid subscript entry #%d: field "subs" must be a cell array.', k );
            elseif fieldalt,
                needfield = 1;
            end
            for j = 1 : length( sb ),
                if isnumeric( sb{j} ),
                    sb{j} = sprintf('%g',sb{j});
                elseif ~ischar( sb{j} ) || size( sb{j}, 1 ) ~= 1,
                    error( 'Invalid subscript entry #%d: invalid cell/array subscript', j );
                end
            end
            sb = sprintf( '%s,', sb{:} );
            y = [ y, tp(1), sb(1:end-1), tp(2) ]; %#ok
        otherwise,
            error( 'Invalid subscript entry #%d: invalid subscript tp: %s', k, tp );
    end
end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [km3] = l2km3(l)
% Convert volume from liters to cubic kilometers.
% Chad Greene 2012
km3 = l*1e-12;

author: Sable

repo: mcbench-benchmarks


text:
function  p = adjustMaxTime(p,maxtime,elapsed)

% Allow the solver at-least 1 second
remaining = maxtime-elapsed;
budgetForSolver = max(1,ceil(remaining));

% TODO: Support more solvers
p.options.cplex.timelimit = budgetForSolver;
p.options.mosek.MSK_DPAR_MIO_MAX_TIME = budgetForSolver;
p.options.gurobi.TimeLimit = budgetForSolver;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function frame = getframe_nosteal_focus(fig, resolution, dpi)
%GETFRAME_NOSTEAL_FOCUS Same as getframe, but does not focus figure
%
%  GETFRAME_NOSTEAL_FOCUS captures the current figure at 1200x900
%  resolution and 96 dpi
%
%  GETFRAME_NOSTEAL_FOCUS(FIG) captures the figure FIG at 1200x900
%  resolution and 96 dpi
%
%  GETFRAME_NOSTEAL_FOCUS(FIG, RESOLUTION) captures the figure FIG
%  at the resolution RESOLUTION, e.g. RESOLUTION = [1200 900], and
%  96 dpi
%
%  GETFRAME_NOSTEAL_FOCUS(FIG, RESOLUTION, DPI) captures the figure FIG
%  at the resolution RESOLUTION, and DPI dots per inch.
%
%See also: getframe
%
%Written by Erik Sällström

    % Set current figure to the figure to capture if not specified
    if (~exist('fig', 'var'))
        fig = gcf;
    end

    % Set default resolution if not set
    if (~exist('resolution', 'var'))
        resolution = [1200, 900];
    end

    % Set default dpi if not set
    if (~exist('dpi', 'var'))
       dpi = 96;
    end

    % Get a filename for a temporary file
    tempfilename = [tempname '.tif'];

    % Set the figure size on paper!
    set(fig, 'PaperPosition', [1 1 resolution/dpi]);

    % Print figure to a file
    print('-dtiffnocompression', ['-r' num2str(dpi)], ['-f' ...
                        num2str(fig)], tempfilename);

    % Load picture
    frame.cdata = imread(tempfilename);
    frame.colormap = [];

    % Delete temporary file
    delete(tempfilename);

end


text:
function y = gskewness(x,n)
%GSKEWNESS Skewness of a grouped sample.
%  In some scientific works, once the data have been gathered from a
%  population of interest, it is often difficult to get a sense of what
%  the data indicate when they are presented in an unorganized fashion.
%  Assembling the raw data into a meaningful form, such as a frequency
%  distribution, makes the data easier to understand and interpret. It is
%  in the context of frequency distributions that the importance of
%  conveying in a succinct way numerical information contained in the data
%  is encountered.
%  So, grouped data is data that has been organized into groups known as
%  classes.  The raw dataset can be organized by constructing a table
%  showing the frequency distribution of the variable (whose values are
%  given in the raw dataset). Such a frequency table is often referred to
%  as grouped data.
%  Here, we developed a m-code to calculate the skewness of a grouped data.
%  One can input the returns or modified vectors n and xout containing the
%  frequency counts and the bin locations of the hist m-function, in a
%  column form matrix.
%  GSKEWNESS(X,0) adjusts the skewness for bias (correction by small sample
%  size). GSKEWNESS(X,1) is the same as GSKEWNESS(X), and does not adjust
%  for bias.
%  If skewness = 0, the data are perfectly symmetrical. But a skewness of
%  exactly zero is quite unlikely for real-world data. Here, we use the
%  Bulmers' rule of thumb criterium (1979), about how we can to
%  interpret the skewness number:
%  - Less than ?1 or greater than +1, the distribution is highly skewed
%  - Between ?1 and ?0.5 or between +0.5 and +1, the distribution is
%    moderately skewed
%  - Between ?0.5 and +0.5, the distribution is approximately symmetric
%
%  Skewness calculation uses the formula,
%
%           g1 =  m3/m2^1.5, do not adjusted for bias
%
%           G1 = SQRT(N*(N - 1))/(N - 2) * g1, adjusted for bias
%
%  where:
%  m2 = second moment of the sample about its mean
%  m3 = third moment of the sample about its mean
%  N  = sample size
%
%  Syntax: function y = gskewness(x,n)
%
%  Inputs:
%       x - data matrix (Size of matrix must be n-by-2; absolut frequency=
%           column 1, class mark=column 2)
%       n - adjusted for bias = 0 (default), do not adjust for bias = 1
%  Outputs:
%       y  - skewness of the values in x
%
%  Example: Suppose we have the next frequency table:
%
%                     ----------------
%                       MC         F
%                     ----------------
%                       61         5
%                       64        18
%                       67        42
%                       70        27
%                       73         8
%                     ----------------
%
%  Taken from: http://www.tc3.edu/instruct/sbrown/stat/shape.htm
%
%  Where we are interested to get the skewness value adjusted for bias.
%
%  Data input:
%  x=[61 5;64 18;67 42;70 27;73 8];
%
%  Calling on Matlab the function:
%          y = gskewness(x,0)
%
%  Answer is:
%
%  y = -0.1098
%
%  Created by A. Trujillo-Ortiz, R. Hernandez-Walls and
%             C.M. Espinosa-Lagunes
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             atrujo@uabc.edu.mx
%
%  Copyright (C)  September 24, 2012.
%
%  To cite this file, this would be an appropriate format:
%  Trujillo-Ortiz, A., R. Hernandez-Walls and C.M. Espinosa-Lagunes. (2012).
%     gskewness:Skewness of a grouped sample. [WWW document].
%     URL http://www.mathworks.com/matlabcentral/fileexchange/
%     38319-gskewness
%
%  References:
%  Bulmer, M. G. (1979), Principles of Statistics. NY:Dover Books on
%             Mathematics.
%  Jayaraman, K. (1999), A Statistical Manual for Foresty Research. Foresty
%             Research Support Programme for Asia and the Pacific. FAO-
%             Corporate Document Repository. Forestry Statistics and Data
%             Collection.
%             URL http://www.fao.org/DOCREP/003/X6831E/X6831E00.HTM
%             PDF ftp://ftp.fao.org/docrep/fao/003/X6831E/X6831E00.pdf
%

if  nargin < 2,
    n = 1; %default
end

c = size(x,2);

if c ~= 2
    error('stats:gskewness:BadData','X must have two colums.');
end

mc = x(:,1); %class mark
f = x(:,2); %absolut frequency
s = sum(f.*mc);
m1 = s/sum(f);
m2 = sum(f.*(mc - m1).^2)/sum(f);
m3 = sum(f.*(mc - m1).^3)/sum(f);
g1 = m3/m2^1.5;

if n == 0;
    y = sqrt(sum(f)*(sum(f) - 1))/(sum(f) - 2) * g1; %skewness adjusted for
                                                     %bias
else n = 1; %default
    y = g1; %skewness not adjusted for bias
end

return,


text:
% SEEMOVIE - see an EEG movie produced by EEGMOVIE
%
% Usage: >> seemovie(Movie,ntimes,Colormap)
%
% Inputs:
%         Movie    = Movie matrix returned by EEGMOVIE
%         ntimes   = Number of times to display {0 -> -10}
%                    If ntimes < 0, movie will play forward|backward
%         Colormap = Color map returned by EEGMOVIE {0 -> default}
%
% Author: Scott Makeig & Colin Humphries, CNL / Salk Institute, 6/3/97
%
% See also: EEGMOVIE

% Copyright (C) 6/3/97 Scott Makeig & Colin Humphries, CNL / Salk Institute, La Jolla CA
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. 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 HOLDER 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.

% 10-31-97 changed test for user-defined colormap -ch & sm
% 1-8-98 added '\n' at end, improved help msg -sm
% 01-25-02 reformated help, added link -ad

function seemovie(Movie,ntimes,Colormap)

fps = 10;   % projetion speed (requested frames per second)

if nargin<1
   help seemovie
   return
end
if nargin<3
    Colormap = 0;
end
if nargin<2
	ntimes = -10;    % default to playing forward|backward endlessly
end
if ntimes == 0
	ntimes = -10;
end

clf
axes('Position',[0 0 1 1]);
if size(Colormap,2) == 3 % if colormap user-defined
	colormap(Colormap)
else
	colormap([jet(64);0 0 0]);    % set up the default topoplot color map
end

if ntimes > 0,
 fprintf('Movie will play slowly once, then %d times faster.\n',ntimes);
else
 fprintf('Movie will play slowly once, then %d times faster forwards and backwards.\n',-ntimes);
end
if abs(ntimes) > 7
	fprintf('   Close figure to abort:  ');
end

%%%%%%%%%%%%%%%%%%%%%%%% Show the movie %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

movie(Movie,ntimes,fps);

%%%%%%%%%%%%%%%%%%%%%%%%    The End     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if abs(ntimes) > 7
	fprintf('\n');
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [varargout] = stdDevDose(planC, varargin)
%Calculate the std dev of the dose (non-variance corrected).
%Stand alone metric :       stdDevDose(planC, structNum, doseNum);
%Request m object   :       stdDevDose(planC, 'getnewmetric');
%Evaluate m object  :       stdDevDose(planC, m, 'evaluate');
%LM:  6 Oct 06, JOD: new metric.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};
optS = planC{indexS.CERROptions};

if(nargin > 1)
    call = varargin{end};
else
    warning('Incorrect Usage, try: stdDevDose(planC, structNum, doseNum);');
    return;
end

switch upper(call)
    case 'GETNEWMETRIC'
        m.name = 'stdDevDose';
        m.valueV = [];
        m.description = 'Returns the standard deviation (non-variance corrected) dose in specified structure';
        m.functionName = @stdDevDose;
        m.note = '';
        m.params(1).name = 'Structure';
        m.params(1).type = 'DropDown';
        m.params(1).list = {planC{indexS.structures}.structureName};
        m.params(1).value = 1;
        m.params(2).name = 'IsTarget?';
        m.params(2).type = 'DropDown';
        m.params(2).list ={'Target', 'Not Target'};
        m.params(2).value = 2;
        m.units = [];
        m.range = [0 inf];
        m.doseSets = [1];
        varargout = {planC,m};
        return;

    case 'EVALUATE'
        m = varargin{1};
        structName = planC{indexS.structures}(m.params(1).value).structureName;
        m.note = structName;
        m.valueV = [];

        maxD = -inf;
        for i=1:length(m.doseSets)
            [planC, doseBinsV, volsHistV] = getDVHMatrix(planC, m.params(1).value, m.doseSets(i));
            ans = calc_stdDevDose(doseBinsV, volsHistV);
            m.valueV = [m.valueV ans];
            maxD = max(maxD, max(planC{indexS.dose}(m.doseSets(i)).doseArray(:)));
        end
        if m.params(2).value == 1
            %Is target, high dose is best.
            m.range = [0 max(m.valueV)];
        elseif m.params(2).value == 2
            m.range = [maxD min(m.valueV)];
        end
        m.units = getDoseUnitsStr(i,planC);
        varargout = {planC,m};
        return;

    otherwise  %No specific call has been made, assume user just wants raw metric.
        if(nargin == 3)
            structNum = varargin{1};
            doseNum = varargin{2};
            [planC, doseBinsV, volsHistV] = getDVHMatrix(planC, structNum, doseNum);
            varargout = {calc_stdDevDose(doseBinsV, volsHistV)};
        else
            error('Not enough parameters, try: stdDevDose(planC, structNum, doseNum);')
            return;
        end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [K, C] = gp_trcov(gp, x1, predcf)
%GP_TRCOV  Evaluate training covariance matrix (gp_cov + noise covariance).
%
%  Description
%    K = GP_TRCOV(GP, TX, PREDCF) takes in Gaussian process GP and
%    matrix TX that contains training input vectors to GP. Returns
%    (noiseless) covariance matrix K for latent values, which is
%    formed as a sum of the covariance matrices from covariance
%    functions in gp.cf array. Every element ij of K contains
%    covariance between inputs i and j in TX. PREDCF is an array
%    specifying the indexes of covariance functions, which are used
%    for forming the matrix. If not given, the matrix is formed
%    with all functions.
%
%    [K, C] = GP_TRCOV(GP, TX, PREDCF) returns also the (noisy)
%    covariance matrix C for observations y, which is sum of K and
%    diagonal term, for example, from Gaussian noise.
%
%  See also
%    GP_SET, GPCF_*
%
% Copyright (c) 2006-2010, 2016 Jarno Vanhatalo
% Copyright (c) 2010 Tuomas Nikoskinen

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

% no covariance functions?
if length(gp.cf)==0 || (nargin>2 && ~isempty(predcf) && predcf(1)==0)
    K=[];
    C=[];
    if nargout>1 && isfield(gp.lik.fh,'trcov')
        C=sparse(0);
        % Add Gaussian noise to the covariance
        C = C + gp.lik.fh.trcov(gp.lik, x1);
        if ~isempty(gp.jitterSigma2)
            C=C+gp.jitterSigma2;
        end
    end
    return
end

[n,m]=size(x1);
ncf = length(gp.cf);
% Evaluate the covariance without noise
K = sparse(n,n);
if isfield(gp,'deriv') && gp.deriv  % derivative observations in use
    ind_Ddim = x1(:,gp.deriv);
    ind_Ddim_derivs = ind_Ddim(ind_Ddim>0);
    uDdim = unique(ind_Ddim_derivs);
    x1 = x1(:,setdiff(1:m,gp.deriv));   % Take only the non-index columns
    if any(strcmp(gp.type,{'FIC' 'PIC' 'PIC_BLOCK' 'CS+FIC' 'VAR' 'DTC' 'SOR'}))
        error('derivative observations have not been implemented for sparse GPs')
    end
end
% check whether predcf is used
if nargin < 3 || isempty(predcf)
    predcf = 1:ncf;
end
% loop through covariance functions
for i=1:length(predcf)
    gpcf = gp.cf{predcf(i)};
    if isfield(gp.lik, 'int_magnitude') && gp.lik.int_magnitude
        if ~isfield(gp,'comp_cf') || (isfield(gp,'comp_cf') && sum(gp.comp_cf{1}==predcf(i)))
            gpcf.magnSigma2=1;
        end
    end

    % derivative observations in use
    if isfield(gp,'deriv') && gp.deriv
        if (~isfield(gpcf, 'selectedVariables') || any(ismember(gpcf.selectedVariables,uDdim)))
            % !!! Note. The check whether to calculate derivative
            % matrices for a covariance function could/should be made
            % nicer
            if size(x1,2) <2    % One dimensional input
                Kff = gpcf.fh.trcov(gpcf, x1(ind_Ddim==0,:));
                Kdf = gpcf.fh.ginput4(gpcf, x1(ind_Ddim==1,:), x1(ind_Ddim==0,:));
                D = gpcf.fh.ginput2(gpcf, x1(ind_Ddim==1,:), x1(ind_Ddim==1,:));

                Kdf=Kdf{1};
                Kfd = Kdf';
                Kdd=D{1};
                % Add all the matrices into a one K matrix
                % K = K + [Kff Kfd; Kdf Kdd];
                Ktemp = [Kff Kfd; Kdf Kdd];
            else
                Ktemp = [];
                % the block of covariance matrix
                Ktemp(ind_Ddim==0,ind_Ddim==0) = gpcf.fh.trcov(gpcf, x1(ind_Ddim==0,:));
                for u1 = 1:length(uDdim)
                    % the blocks on the left side, below Kff
                    if sum(ind_Ddim==0)>0
                        Kdf = gpcf.fh.ginput4(gpcf, x1(ind_Ddim==uDdim(u1),:), x1(ind_Ddim==0,:), uDdim(u1));
                        Ktemp(ind_Ddim==uDdim(u1),ind_Ddim==0) = Kdf{1};
                        Ktemp(ind_Ddim==0,ind_Ddim==uDdim(u1)) = Kdf{1}';
                    end
                    D = gpcf.fh.ginput2(gpcf, x1(ind_Ddim==uDdim(u1),:), x1(ind_Ddim==uDdim(u1),:), uDdim(u1));
                    Ktemp(ind_Ddim==uDdim(u1),ind_Ddim==uDdim(u1)) = D{1};

                    uDdim2 = uDdim(u1+1:end);
                    for u2=1:length(uDdim2)
                        Kdf2 = gpcf.fh.ginput3(gpcf, x1(ind_Ddim==uDdim(u1),:) ,x1(ind_Ddim==uDdim2(u2),:), uDdim(u1), uDdim2(u2));
                        Ktemp(ind_Ddim==uDdim(u1),ind_Ddim==uDdim2(u2)) = Kdf2{1};
                        Ktemp(ind_Ddim==uDdim2(u2),ind_Ddim==uDdim(u1)) = Kdf2{1}';
                    end

                end
            end
        else
            Ktemp = zeros(n,n);
            Ktemp(ind_Ddim==0,ind_Ddim==0) = gpcf.fh.trcov(gpcf, x1(ind_Ddim==0,:));
        end
        K= K+ Ktemp;
    else
        % Regular GP without derivative observations
        K = K + gpcf.fh.trcov(gpcf, x1);
    end
end

n = size(K,1);
n1 = n+1;
if ~isempty(gp.jitterSigma2)
    if issparse(K)
        K = K + sparse(1:n,1:n,gp.jitterSigma2,n,n);
    else
        K(1:n1:end)=K(1:n1:end) + gp.jitterSigma2;
    end
end
if nargout>1
    C=K;
    if isfield(gp.lik.fh,'trcov')
        C = C + gp.lik.fh.trcov(gp.lik, x1);
    end
end

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
% [INPUT]
% f = A valid figure instance.

function maximize_figure(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('f',@(x)validateattributes(x,{'matlab.ui.Figure'},{'scalar'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    f = ipr.f;

    nargoutchk(0,0);

    maximize_figure_internal(f);

end

function maximize_figure_internal(f)

    pause(0.1);

    if (verLessThan('MATLAB','9.4'))
        frame = get(f,'JavaFrame'); %#ok<JAVFM>
        set(frame,'Maximized',true);
    else
        set(f,'WindowState','maximized');
    end

    drawnow();

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
%kkrect2pol 'Convert (real, imaginary) to (r, theta) Coordinates'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros krect2pol.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kkrect2pol(i, {'i','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% krect2pol - Convert (real, imaginary) to (r, theta) Coordinates
%
%  DESCRIPTION
% The "Rect. to Polar" operator converts data from the rectangular to the
% polar coordinate system for each complex data point in the input data
% object, \fBInput".  The conversion from rectangular to polar will result
% in angles measured in radians.  The data type of
% the output data object, \fBOutput", will be the same as the data type
% of \fBInput".
%
% Executing "Rect. to Polar" runs the program \fIkcmplx\fP
% with the -r2p flag.
%
%  "Map Data" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_1input
%
%  "Validity Mask" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "Location and Time Data" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::kcmplx, DATAMANIP::kpol2rect
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkrect2pol(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkrect2pol(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output'};
maxval={0,0};
minval={0,0};
istoggle=[0,0];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kcmplx"  -r2p'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function distplot(y,cutoff,class,nid)

%DISTPLOT plots the vector y versus the index.
% At the height of the cutoff value a red vertical line is plotted.
%
% Required input arguments:
%       y  : the vector to be plotted
%  cutoff  : a cutoff value
%
% Optional input arguments:
%   class  : the class of the y-vector
%     nid  : number of points to be identified (Default value: 3)
%
% I/O: distplot(y,cutoff,class,nid)
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Last update: 24/11/2003

set(gcf,'Name', 'Index plot of the distances', 'NumberTitle', 'off');
n=length(y);
if nargin==2
    laby='Distance';
    nid=3;
elseif nargin==3
    nid=3;
end
ymax=max([max(y),cutoff,2.5])*1.05;
plot(1:n,y,'o')
box on
xlabel('Index')
if strcmp(class,'MCDCOV')
    laby='Robust distance';
elseif strcmp(class,'COV')
    laby='Mahalanobis distance';
end
ylabel(laby)
xlim([-0.025*n,n*1.05]);
ylim([-0.025*ymax,ymax]);
plotnumbers(1:n,y,0,nid,1)
line([-0.025*n,n*1.05],repmat(max([cutoff,2.5]),1,2),'Color','r');
title([class]);

author: marianux

repo: ecg-kit

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


text:
classdef Tests < matlab.unittest.TestCase
    % SolverTest tests solutions to the quadratic equation
    % a * x^2 + b * x + c = 0

    methods (Test)
        function testVehicleSimple(testCase)
            % Simulation time
            T = 6;                      % Total simulation time [s]
            resol = 50;                 % Resolution
            TSPAN = 0:T/resol:T;        % Time span [s]

            %% Tire parameters
            % Chosen tire: <TirePacejka.html TirePacejka.m>.
            %

            TireModel = VehicleDynamicsLateral.TirePacejka();

            %% Vehicle parameters
            % Chosen Vehicle: <VehicleSimpleNonlinear.html VehicleSimpleNonlinear.m>.

            System = VehicleDynamicsLateral.VehicleSimpleNonlinear();
            System.tire = TireModel;

            % Initial conditions
            System.dPSI0 = 0.7;                % Initial yaw rate [rad/s]
            System.ALPHAT0 = -0.2;             % Initial side slip angle [rad]
            System.PSI0 = 0;                   % Initial yaw angle [rad]
            System.X0 = 0;                     % Initial CG horizontal position [m]
            System.Y0 = 0;                     % Initial CG vertical position [m]
            System.V0 = 20;                    % Initial CG velocity [m/s]

            simulator = VehicleDynamicsLateral.Simulator(System, TSPAN);
            simulator.Simulate();

            % Retrieving states
            dPSI = simulator.dPSI;
            ALPHAT = simulator.ALPHAT;
            PSI = simulator.PSI;
            XT = simulator.XT;
            YT = simulator.YT;
            VEL = simulator.VEL;
        end
    end
end

%% See Also
%
% <../../index.html Home>
%

author: andresmendes

repo: Vehicle-Dynamics-Lateral


text:
%--------------------------------------------------------------------------
%
% Example script for sinf_1D.m and sinf_3D.m
%
% Author        : dr.ir. Emanuel A.P. Habets
% Date          : 16-09-2010
%
% Related paper : E.A.P. Habets and S. Gannot, 'Generating sensor signals
%                 in isotropic noise fields', Submitted to the Journal
%                 of the Acoustical Society of America, May, 2007.
%
% Comment       :
%
%--------------------------------------------------------------------------

clear;
close all;

%% Initialization
fs = 8000;                       % Sample frequency
NFFT = 256;                      % Number of frequency bins (for analysis)
w = 2*pi*fs*(0:NFFT/2)/NFFT;
c = 340;                         % Speed of sound
L = 2^18;                        % Data length

[x1,y1,z1]=sph2cart(0,0,0.1);    % Sensor position 1
[x2,y2,z2]=sph2cart(0,0,0.2);    % Sensor position 2
P = [0 x1 x2; 0 y1 y2; 0 z1 z2]; % Construct position matrix
M = 3;                           % Number of sensors

% Calculate sensor distances w.r.t. sensor 1
d = zeros(1,M);
for m = 2:M
    d(m) = norm(P(:,m)-P(:,1),2);
end

%% Generate sensor signals
params.c = c;
params.fs = fs;

% 1D example
params.N_phi = 64;
z = sinf_1D(d,L,params);

% 3D example
% params.N = 256;
% z = sinf_3D(P,L,params);

%% Calculate spatial coherences
sc_sim = zeros(M-2,NFFT/2+1);
sc_theory = zeros(M-2,NFFT/2+1);
for m = 1:M-1
    [sc,F]=mycohere(z(1,:)',z(m+1,:)',NFFT,fs,hanning(NFFT),0.75*NFFT);
    sc_sim(m,:) = real(sc');

    sc_theory(m,:) = sinc(w*d(m+1)/c/pi);
end

%% Plot results
% Sensor pair 1-2
figure(1);
m=1;
plot(F/1000,sc_sim(m,:),'k')
hold on;
plot(F/1000,sc_theory(m,:),'--k')
hold off;
xlabel('Frequency [kHz]');
ylabel('Spatial Coherence');
title(sprintf('Distance %1.2f m',d(m+1)));
set(gca,'DataAspectRatio',[1 0.75 1]);
legend('Simulation','Theory');
grid on;

% Sensor pair 1-3
figure(2);
m=2;
plot(F/1000,sc_sim(m,:),'k')
hold on;
plot(F/1000,sc_theory(m,:),'--k')
hold off;
xlabel('Frequency [kHz]');
ylabel('Spatial Coherence');
title(sprintf('Distance %1.2f m',d(m+1)));
set(gca,'DataAspectRatio',[1 0.75 1]);
legend('Simulation','Theory');
grid on;

author: ZitengWang

repo: MASP

save_path: github-repos/MATLAB/ZitengWang-MASP


text:
%  WaveForm DataBase (WFDB) Toolbox
%
%
%This is a set of MATLAB functions and wrappers for reading, writing, and processing
%files in the formats used by PhysioBank databases (among others).
%The WFDB Toolbox has support for reading public PhysioNet databases directly from
%web. This feature allows your code to analyze a wide range of physiological
%signals available from PhysioBank without the need to download entire
%records and to store them locally. This toolbox is distributed under the LGPL
%license (see LICENSE.txt file in this directory). For more informationa about the
%toolbox please go to: http://www.physionet.org
%
%
%
%
%  Table of Contents (T0C)
% -----------------------
%   ann2rr          - Extract a list of intervals from an annotation file
%   bxb             - ANSI/AAMI-standard beat-by-beat annotation comparator
%   dfa             - Detrended Fluctuation Analysis
%   corrint         - Correlation Integral Analysis
%   edr             - Derives a respiration signal from an ECG signal
%   ecgpuwave       - Estimation of QRS and P waves from ECG signals
%   gqrs            - Estimation of QRS from ECG signals
%   lomb	        - Estimates power spectrum using the Lomb periodogram method
%   mat2wfdb        - Writes a MATLAB variable into a WDFB record file
%   mrgann          - Merge annotation files
%   msentropy       - Multi scale entropy estimation
%   physionetdb     - Get information about all of PhysioNet's available databases and signals
%   pbsearch        - Search PhysioBank from MATLAB's browser for records with specific features.
%   rdann           - Read annotation files for WFDB records
%   rdmat           - Reads a signal into the workspace from a *.mat file generated by WFDB2MAT
%   rdmimic2wave    - Searches MIMIC II matched waveform records within a clinical time range
%   rdsamp          - Read signal files of WFDB records
%   sortann         - Rearrange annotations in canonical order
%   snip            - Copy an excerpt of a WFDB record
%   sqrs            - Finds the QRS complexes of a WFDB ECG record signal
%   surrogate       - Generates amplitude adjusted phase shuffled surrogate time series
%   sumann	        - Summarize the contents of a WFDB annotation file
%   tach            - Calculates instantaneous heart rate of a WFDB ECG record signal
%   visgraph        - Visibility Graph Analysis
%   wfdb            - Prints this help information of the Toolbox
%   wfdb2mat        - Converts a record *.dat file into a *.mat file.
%   wfdbexec        - Executes a system call to any installed native WFDB command.
%   wabp	        - Arterial blood pressure (ABP) pulse detector
%   wfdbdemo        - Demonstration of the WFDB App Toolbox
%   wfdbdesc        - Return signal information for about a WFDB record
%   wfdbloadlib     - Load WFDB library and displays Toolbox configuration.
%   wfdbRecordViewer- GUI for visualizing WFDB records and annotations
%   wfdbtest        - Script to test installation of the Toolbox
%	wfdbtime        - Converts sample index to WFDB Time based on WFDB record information
%	wfdbtool        - Launches PhysioNet's Signal and Annotation Viewer ( LightWave )
%   woody           - Perform signal averaging with alignment
%   wqrs            - Finds the QRS complexes of a WFDB ECG record signal
%   wrann           - Writes annotations for WFDB records into annotation files
%   wrsamp          - Writes signal data into WFDB-compatible records
%
%
%
%   To credit this toolbox, please cite the following references in your work:
%
%   Silva, I, Moody, G.
%   "An Open-source Toolbox for Analysing and Processing PhysioNet Databases in MATLAB and Octave."
%   Journal of Open Research Software 2(1):e27
%   [http://dx.doi.org/10.5334/jors.bi] ;
%   2014 (September 24).
%
%
%   Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG, Mietus JE, Moody GB, Peng CK, Stanley HE.
%   "PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource for Complex
%   Physiologic Signals."
%   Circulation 101(23):e215-e220
%   [http://circ.ahajournals.org/cgi/content/full/101/23/e215];
%   2000 (June 13).
%   PMID: 10851218; doi: 10.1161/01.CIR.101.23.e215
%
%
%   In addition, some of these functions use binary executables compiled
%   from open-source third-party code contributed to PhysioNet. When using
%   these functions on your work, please look at the help for that function
%   in order find out how to credit the original paper and authors.
%
%   For questions, contributions, and feedback post at our community Forum:
%       http://groups.google.com/forum/#!forum/wfdb-app-toolbox
%
%
%   The source code for the native libraries used in this toolbox can be obtained from PhysioNet under
%   the GNU GPL agreement.
%
%   The original contributors of any open source native code that is available at PhysioNet
%   and is used by the Toolbox are credited in their respective MATLAB wrappers. In addition, the
%   following people have contributed to the development or testing of the MATLAB wrappers
%   and the JVM interface:
%
%   Sahar Alkhairy
%   Fernando Andreotti Lage
%   Joachim Behar
%   Eudald Bogatell
%   Jonas Carlson
%   Gari D. Clifford
%   Michael Craig
%   Mohammad Ghassemi
%   Alistair Johnson
%   Li-wei Lehman
%   Sara Mariani
%   Louis Mayaud
%   Blaž Merela
%   Benjamin Moody
%   George B. Moody
%   Shamin Nemati
%   Piotr Podziemski
%   Erina Katsumata
%   Aled Rowen
%   Kari Schoonbee
%   Daniel J. Scott
%   Ikaro Silva
%   Gabriel Squillace
%   Bryan Tripp
%   Chen Xie
%   Tingting Zhu
%
%
%Created by Ikaro Silva 2012

%endOfHelp

author: ikarosilva

repo: wfdb-app-toolbox


text:
classdef DirichletD
%%DIRICHLETD Functions to handle the Dirichlet distribution.
%Implemented methods are: mean, cov, PDF, rand
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)

function val=mean(alph)
%%MEAN Obtain the mean of the Dirichlet distribution for the given
%      concentration parameters.
%
%INPUTS: alph An NX1 vector of positive, real concentration parameters.
%             The length of alpha determines the length of the random
%             vector. If one wishes for the mean of all of the parameters
%             to be equal, then just set alph to all ones.
%
%OUTPUTS: val  The mean of the Dirichlet distribution under consideration.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=alph/sum(alph);
end

function covMat=cov(alph)
%%VAR Obtain the covariance matrix of the Dirichlet distribution for the
%     given concentration parameters.
%
%INPUTS: alph An NX1 vector of positive, real concentration parameters. The
%             length of alpha determines the length of the random vector.
%             If one wishes for the mean of all of the parameters to be
%             equal, then just set alph to all ones.
%
%OUTPUTS: covMat The covariance matrix of the Dirichlet distribution under
%                consideration.
%
%Note that since the elements of the Dirichlet random variable must sum to
%one (i.e. are somewhat dependent), the covariance matrix is singular.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    numVals=length(alph);
    alpha0=sum(alph);
    denom=alpha0^2*(alpha0+1);

    %First, fill in the diagonal elements.
    covMat=diag(alph.*(alpha0-alph)/denom);

    %Then, fill in the off-diagonal elements.
    for curRow=1:numVals
       for curCol=(curRow+1):numVals
           covMat(curRow,curCol)=-alph(curRow)*alph(curCol)/denom;
           covMat(curCol,curRow)=covMat(curRow,curCol);
       end
    end
end

function val=PDF(X,alph)
%%PDF Evaluate the PDF of the Dirichlet distribution at a particular vector
%     for the given concentration parameters.
%
%INPUTS: X An NX1 dirichlet random variable. Its elements must sum to one.
%     alph An NX1 vector of positive, real concentration parameters. If one
%          wishes for the mean of all of the parameters to be equal, then
%          just set alph to all ones.
%
%OUTPUTS: val  The value of the Dirichlet distribution at X given alph.
%
%The Dirichlet distribution is given in Appendix B of [1].
%
%REFERENCES:
%[1] C. M. Bishop, Pattern Recognition and Machine Learning. Cambridge,
%    United Kingdom: Springer, 2007.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    C=gamma(sum(alph))/prod(gamma(alph));
    val=C*prod(X.^(alph-1));
end

function val=rand(alph)
%%RAND Generate a Dirichlet random vector with concentration vector
%      parameter vector alpha.
%
%INPUTS: alph An NX1 vector of positive, real concentration parameters. The
%             length of alpha determines the length of the random vector
%             generated. If one wishes for the mean of all of the
%             parameters to be equal, then just set alph to all ones.
%
%OUTPUTS: val An NX1 Dirichlet random vector.
%
%Dirichlet random vectors have elements that sum to one and whose
%elements are bounded between 0 and 1. Thus, Dirichlet random vectors are
%useful for generating weights for random mixture distributions.
%
%The Dirichlet random vector is generated by transforming a set of gamma
%distributed random variables. The relationship between the Dirichlet
%distribution and the gamma distribution is discussed in [1].
%
%REFERENCES:
%[1] B. A. Frigyik, A. Kapila, and M. R. Gupta, "Introduction to the
%    Dirichlet Distribution and Related processes", University of
%    Washington Technical Report, No. UWEETR-20-0006, Dec. 2010.
%    https://www.ee.washington.edu/techsite/papers/documents/UWEETR-2010-0006.pdf
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    numVals=length(alph);
    gVals=zeros(numVals,1);
    for curVal=1:numVals
       gVals(curVal)=randGamma(1,alph(curVal),1);
    end

    val=gVals/sum(gVals);
end

end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%this is the code for AAAI 2016 paper: Top-N Recommender System via Matrix
%Completion by zhao kang,chong peng,qiang cheng.

%[hr,arhr,X] = mc_logdet(data,mu,rho)
function [X] = demo(data,mu,rho,toler,maxiter)
  [m,n]=size(data);
  M=data;
  id=find(M==0);
  ID=ones(m,n);
  ID(id)=0;
  [X] = MC_LogDet_v3(M,ID,mu,rho,toler,maxiter);
end

author: andrewssobral

repo: lrslibrary


text:
function P=perimeter(r,N)
% find perimeter of the polygon

n=1:N-1;
r1=r(:,n); % current points
r2=r(:,n+1); % next points
dr=r1-r2; % difference
dr2=sum(dr.^2); % squared lenghts
drl=sqrt(dr2); % lenghts
P=sum(drl); % perimeter

% last and first:
r1=r(:,N); % current points
r2=r(:,1); % next points
dr=r1-r2; % difference
dr2=sum(dr.^2); % squared lenghts
drl=sqrt(dr2); % lenghts
P=P+drl; % perimeter


text:
function jed = ymdf_to_jed_jelali ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_JED_JELALI converts a Jelali YMDF date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, real JED, the corresponding Julian Ephemeris Date.
%
  jed_epoch = epoch_to_jed_jelali ( );

  jed = jed_epoch + ( d - 1 ) + 30 * ( m - 1 ) + 365 * ( y - 1 ) ...
    + floor ( ( y - 1 ) / 4 ) + f;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x y] = BipolarCoordinates(Xi,Eta,a)

% a is radii that defies height of the domain
% a > 0

% Check the input arguments and set default values

if nargin == 2
    a = 1. ;
end
r = Xi ;
s = pi*(Eta-1/2) ;

x = (a*sinh(r))/(cosh(r)+cos(s)) ;
y = (a*sin(s))/(cosh(r)+cos(s)) ;


text:
function a = pascal2 ( n )

%*****************************************************************************80
%
%% PASCAL2 returns the PASCAL2 matrix.
%
%  Formula:
%
%    If ( I = 1 or J = 1 )
%      A(I,J) = 1
%    else
%      A(I,J) = A(I-1,J) + A(I,J-1)
%
%  Example:
%
%    N = 5
%
%    1 1  1  1  1
%    1 2  3  4  5
%    1 3  6 10 15
%    1 4 10 20 35
%    1 5 15 35 70
%
%  Properties:
%
%    A is a "chunk" of the Pascal binomial combinatorial triangle.
%
%    A is positive definite.
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    A is nonsingular.
%
%    det ( A ) = 1.
%
%    A is unimodular.
%
%    Eigenvalues of A occur in reciprocal pairs.
%
%    The condition number of A is approximately 16^N / ( N*PI ).
%
%    The elements of the inverse of A are integers.
%
%    A(I,J) = (I+J-2)! / ( (I-1)! * (J-1)! )
%
%    The Cholesky factor of A is a lower triangular matrix R,
%    such that A = R * R'.  The matrix R is a Pascal
%    matrix of the type generated by subroutine PASCAL.  In other
%    words, PASCAL2 = PASCAL * PASCAL'.
%
%    If the (N,N) entry of A is decreased by 1, the matrix is singular.
%
%    Gregory and Karney consider a generalization of this matrix as
%    their test matrix 3.7, in which every element is multiplied by a
%    nonzero constant K.  They point out that if K is the reciprocal of
%    an integer, then the inverse matrix has all integer entries.
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Brawer, Magnus Pirovino,
%    The linear algebra of the Pascal matrix,
%    Linear Algebra and Applications,
%    Volume 174, 1992, pages 13-23.
%
%    Robert Gregory, David Karney,
%    Example 3.7,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, 1969, page 32,
%    LC: QA263.G68.
%
%    Nicholas Higham,
%    Accuracy and Stability of Numerical Algorithms,
%    Society for Industrial and Applied Mathematics,
%    Philadelphia, PA, USA, 1996; section 26.4.
%
%    Sam Karlin,
%    Total Positivity, Volume 1,
%    Stanford University Press, 1968.
%
%    Morris Newman, John Todd,
%    The evaluation of matrix inversion programs,
%    Journal of the Society for Industrial and Applied Mathematics,
%    Volume 6, Number 4, pages 466-476, 1958.
%
%    Heinz Rutishauser,
%    On test matrices,
%    Programmation en Mathematiques Numeriques,
%    Centre National de la Recherche Scientifique,
%    1966, pages 349-365.
%
%    John Todd,
%    Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
%    Academic Press, 1977, page 172.
%
%    HW Turnbull,
%    The Theory of Determinants, Matrices, and Invariants,
%    Blackie, 1929.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n

      if ( i == 1 )
        a(i,j) = 1.0;
      elseif ( j == 1 )
        a(i,j) = 1.0;
      else
        a(i,j) = a(i,j-1) + a(i-1,j);
      end

    end
  end

  return
end


text:
function level = optiWarnLevel(warn)
%Return opti warning level in numerical form

%Get Warning Level
if(strcmpi(warn,'all'))
    level = 2;
elseif(strcmpi(warn,'critical'))
    level = 1;
else
    level = 0;
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function installSegmentation(query, keepAllNodes, filePaths, numGrayLayers)
%
% function installSegmentation(query, keepAllNodes, filePaths, numGrayLayers)
%
% Rebuild coords after a new segmentation.
%
% Note: could automatically rebuild Gray/corAnals by xforming from inplane, but we
% would have to do this for each dataType.
%
% 12/8/00 Wandell and Huk
% djh, 2/14/2001
%    Uses hiddenGray and hiddenFlat structures rather than opening windows.
%    Uses fullfile rather than cd'ing into the directory
% ras, dy, al 01/09: added option to keep all nodes from each hemisphere.
% ras 06/09: eliminated a redundant set of file requests. Streamlined the
% logic, so it doesn't call a large series of nested files. Condensed file
% selection to a single dialog, which guesses reasonable defaults.

if ~exist('query','var')
	% query only if Gray, Volume, or Flat files are found
	global HOMEDIR %#ok<TLEV> %
	w1 = dir( fullfile(HOMEDIR, 'Volume*') );
	w2 = dir( fullfile(HOMEDIR, 'Gray*') );
	w3 = dir( fullfile(HOMEDIR, 'Flat*') );
	if (~isempty(w1) || ~isempty(w2) || ~isempty(w3))
		query = 1;
	else
		query = 0;
	end
end

% if keepAllNodes is true, we store coordinates for all gray voxels. If it
% is false, then we only store gray voxels within the functional field of
% view. The latter is more resource efficient. The former has greater
% flexibility, for example to visualize ROIs on a mesh that extend beyond
% the field of view from a particular session.
if ~exist('keepAllNodes','var'), keepAllNodes=0; end

if notDefined('filePaths')
	% get a set of file paths in a dialog
	filePaths = segmentationFilesDialog;
    if isempty(filePaths), return; end
end

if ~exist('numGrayLayers', 'var')
    resp = inputdlg('number of gray layers?', 'enter a number', 1, {'3'});
    numGrayLayers = str2double(resp);
end

if query
    resp = questdlg('Delete Gray and all Flat data files and close all Gray and Flat windows?');
else
    resp = 'Yes';
end

% delete coords and rebuild
if ~strcmp(resp,'Yes')
	disp('Coords, corAnal, and parameter map files not deleted.');
    return;
end
closeAllGrayWindows;
closeAllFlatWindows;

% Deletes coords, corAnal, and other parameter map files
cleanGray;
cleanAllFlats;

% ras 06/09: initially the logic was different here. The function would
% call 'initHiddenGray', which in turn would call 'switch2Gray', which
% would call 'getGrayCoords', which would detect that you wanted to
% rebuild, then call 'buildGrayCoords'. This was a little confusing, but
% moreover didn't work if you wanted to keep all nodes (most of the
% intermediate functions don't allow such a parameter, and it wouldn't make
% much sense most of the time). So, I make it more direct here.

% build the gray coordinates fields (this saves the coords.mat file)
if keepAllNodes==1
	% the gray nodes/edges are either loaded or built in the call to
	% initHiddenGray above. It is awkward to pass a flag to this
	% function to load all nodes, so we recompute the nodes/edges here.
    if prefsVerboseCheck >= 1
        fprintf('[%s]: Rebuilding Gray coordinates with ALL nodes. \n', ...
            mfilename);
    end
else
    if prefsVerboseCheck >= 1
        fprintf(['[%s]: Rebuilding Gray coordinates with nodes only ' ...
            'where the inplanes have data. \n'], mfilename);
    end
end
buildGrayCoords([], [], keepAllNodes, filePaths, numGrayLayers);

disp('Rebuilt Gray coords. Deleted old corAnal and parameter map files.');

return;
% /-----------------------------------------------------------/ %

% /-----------------------------------------------------------/ %
function filePaths = segmentationFilesDialog
% dialog to get the file paths for installing a segmentation.

%% first, guess default paths
defaultClass = '';

try

    pattern = fullfile(pwd, '3DAnatomy', '*class*nii*');
	w = dir(pattern);
	if ~isempty(w),
		defaultClass = fullfile(pwd, '3DAnatomy', w(1).name);
	end

catch ME
    warning(ME.identifier, ME.mesage);
	% can't get it? don't worry...

end

%% now build the dialog structure
dlg(1).fieldName	= 'classFile';
dlg(end).style		= 'filename';
dlg(end).string		= 'Classification File (NIFTI)?';
dlg(end).value		= defaultClass;

%% put the dialog to the user
[resp, ok] = generalDialog(dlg, 'Install Gray Segmentation');
if ~ok
	filePaths = [];
    fprintf('User aborted.\n');
    return;
end

%% parse the response
filePaths = resp.classFile;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [g, c, tmp, m]=kMeansCluster(m,k,isRand)
%%%%%%%%%%%%%%%%
%
% kMeansCluster - Simple k means clustering algorithm
% Author: Kardi Teknomo, Ph.D.
%
% Purpose: classify the objects in data matrix based on the attributes
% Criteria: minimize Euclidean distance between centroids and object points
% For more explanation of the algorithm, see http://people.revoledu.com/kardi/tutorial/kMean/index.html
% Output: matrix data plus an additional column represent the group of each object
%
% Example: m = [ 1 1; 2 1; 4 3; 5 4]  or in a nice form
%          m = [ 1 1;
%                2 1;
%                4 3;
%                5 4]
%          k = 2
% kMeansCluster(m,k) produces m = [ 1 1 1;
%                                   2 1 1;
%                                   4 3 2;
%                                   5 4 2]
% Input:
%   m      - required, matrix data: objects in rows and attributes in columns
%   k      - optional, number of groups (default = 1)
%   isRand - optional, if using random initialization isRand=1, otherwise input any number (default)
%            it will assign the first k data as initial centroids
%
% Local Variables
%   f      - row number of data that belong to group i
%   c      - centroid coordinate size (1:k, 1:maxCol)
%   g      - current iteration group matrix size (1:maxRow)
%   i      - scalar iterator
%   maxCol - scalar number of rows in the data matrix m = number of attributes
%   maxRow - scalar number of columns in the data matrix m = number of objects
%   temp   - previous iteration group matrix size (1:maxRow)
%   z      - minimum value (not needed)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<3,        isRand=0;   end
if nargin<2,        k=1;        end
tmp = [];

[maxRow, maxCol]=size(m);
if maxRow<=k,
    y=[m, 1:maxRow];
else

	% initial value of centroid
    if isRand,
        p = randperm(size(m,1));      % random initialization
        for i=1:k
            c(i,:)=m(p(i),:)
    	end
    else
        for i=1:k
           c(i,:)=m(i,:);        % sequential initialization
    	end
    end

	temp=zeros(maxRow,1);   % initialize as zero vector

	while 1,
        d=DistMatrix(m,c);  % calculate objcets-centroid distances
        [z,g]=min(d,[],2);  % find group matrix g
        if g==temp,
            break;          % stop the iteration
        else
            temp=g;         % copy group matrix to temporary variable
        end
        for i=1:k
            f=find(g==i);
            if f            % only compute centroid if f is not empty
                c(i,:)=mean(m(find(g==i),:),1);
            end
        end
	end

	y=[m,g];

end

%The Matlab function kMeansCluster above call function DistMatrix as shown in the code below. It works for multi-dimensional Euclidean distance. Learn about other type of distance here.

    function d=DistMatrix(A,B)
                 %%%%%%%%%%%%%%%%%%%%%%%%%
                 % DISTMATRIX return distance matrix between points in A=[x1 y1 ... w1] and in B=[x2 y2 ... w2]
                 % Copyright (c) 2005 by Kardi Teknomo,  http://people.revoledu.com/kardi/
                 %
                 % Numbers of rows (represent points) in A and B are not necessarily the same.
                 % It can be use for distance-in-a-slice (Spacing) or distance-between-slice (Headway),
                 %
                 % A and B must contain the same number of columns (represent variables of n dimensions),
                 % first column is the X coordinates, second column is the Y coordinates, and so on.
                 % The distance matrix is distance between points in A as rows
                 % and points in B as columns.
                 % example: Spacing= dist(A,A)
                 % Headway = dist(A,B), with hA ~= hB or hA=hB
                 %          A=[1 2 3; 4 5 6; 2 4 6; 1 2 3]; B=[4 5 1; 6 2 0]
                 %          dist(A,B)= [ 4.69   5.83;
                 %                       5.00   7.00;
                 %                       5.48   7.48;
                 %                       4.69   5.83]
                 %
                 %          dist(B,A)= [ 4.69   5.00     5.48    4.69;
                 %                       5.83   7.00     7.48    5.83]
                 %%%%%%%%%%%%%%%%%%%%%%%%%%%
                 [hA,wA]=size(A);
                 [hB,wB]=size(B);
                 if wA ~= wB,  error(' second dimension of A and B must be the same'); end
                 for k=1:wA
                      C{k}= repmat(A(:,k),1,hB);
                      D{k}= repmat(B(:,k),1,hA);
                 end
                 S=zeros(hA,hB);
                 for k=1:wA
                      S=S+(C{k}-D{k}').^2;
                 end
                 d=sqrt(S);

author: PatternRecognition

repo: OpenBMI


text:
function [ triangle_num, triangle_node, triangle_neighbor ] = r8tris2 ( ...
  node_num, node_xy )

%*****************************************************************************80
%
%% R8TRIS2 constructs a Delaunay triangulation of 2D vertices.
%
%  Discussion:
%
%    The routine constructs the Delaunay triangulation of a set of 2D vertices
%    using an incremental approach and diagonal edge swaps.  Vertices are
%    first sorted in lexicographically increasing (X,Y) order, and
%    then are inserted one at a time from outside the convex hull.
%
%  Modified:
%
%    07 February 2005
%
%  Author:
%
%    Original FORTRAN77 version by Barry Joe,
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Barry Joe,
%    GEOMPACK - a software package for the generation of meshes
%    using geometric algorithms,
%    Advances in Engineering Software,
%    Volume 13, pages 325-331, 1991.
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of the nodes.
%
%    Output, integer TRIANGLE_NUM, the number of triangles in the triangulation;
%    TRIANGLE_NUM is equal to 2*NODE_NUM - BOUNDARY_NUM - 2, where BOUNDARY_NUM
%    is the number of boundary vertices.
%
%    Output, integer TRIANGLE_NODE(3,TRIANGLE_NUM), the nodes that make up each triangle.
%    The elements are indices of P.  The vertices of the triangles are
%    in counter clockwise order.
%
%    Output, integer TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), the triangle neighbor list.
%    Positive elements are indices of TIL; negative elements are used for links
%    of a counter clockwise linked list of boundary edges; LINK = -(3*I + J-1)
%    where I, J = triangle, edge index; TRIANGLE_NEIGHBOR(J,I) refers to
%    the neighbor along edge from vertex J to J+1 (mod 3).
%
  triangle_num = 0;
  triangle_node = [];
  triangle_neighbor = [];

  tol = 100.0 * eps;
%
%  Sort the vertices by increasing (x,y).
%
  indx = r82vec_sort_heap_index_a ( node_num, node_xy );

  node_xy = r82vec_permute ( node_num, node_xy, indx );
%
%  Make sure that the data points are "reasonably" distinct.
%
  m1 = 1;

  for i = 2 : node_num

    m = m1;
    m1 = i;

    k = 0;

    for j = 1 : 2

      cmax = max ( abs ( node_xy(j,m) ), abs ( node_xy(j,m1) ) );

      if ( tol * ( cmax + 1.0 ) < abs ( node_xy(j,m) - node_xy(j,m1) ) )
        k = j;
        break
      end

    end

    if ( k == 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8TRIS2 - Fatal error!\n' );
      fprintf ( 1, '  Fails for point number I = %d\n', i );
      fprintf ( 1, '  M = %d\n', m );
      fprintf ( 1, '  M1 = %d\n', m1 );
      fprintf ( 1, '  X,Y(M)  = %f  %f\n', node_xy(1,m), node_xy(2,m) );
      fprintf ( 1, '  X,Y(M1) = %f  %f\n', node_xy(1,m1), node_xy(2,m1) );
      error ( 'R8TRIS2 - Fatal error!' )
      return
    end

  end
%
%  Starting from points M1 and M2, search for a third point M that
%  makes a "healthy" triangle (M1,M2,M)
%
  m1 = 1;
  m2 = 2;
  j = 3;

  while ( 1 )

    if ( node_num < j )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8TRIS2 - Fatal error!\n' );
      error ( 'R8TRIS2 - Fatal error!' )
      return
    end

    m = j;

    lr = lrline ( node_xy(1,m), node_xy(2,m), node_xy(1,m1), node_xy(2,m1), ...
      node_xy(1,m2), node_xy(2,m2), 0.0 );

    if ( lr ~= 0 )
      break
    end

    j = j + 1;

  end
%
%  Set up the triangle information for (M1,M2,M), and for any other
%  triangles you created because points were collinear with M1, M2.
%
  triangle_num = j - 2;

  if ( lr == -1 )

    triangle_node(1,1) = m1;
    triangle_node(2,1) = m2;
    triangle_node(3,1) = m;
    triangle_neighbor(3,1) = -3;

    for i = 2 : triangle_num

      m1 = m2;
      m2 = i+1;
      triangle_node(1,i) = m1;
      triangle_node(2,i) = m2;
      triangle_node(3,i) = m;
      triangle_neighbor(1,i-1) = -3 * i;
      triangle_neighbor(2,i-1) = i;
      triangle_neighbor(3,i) = i - 1;

    end

    triangle_neighbor(1,triangle_num) = -3 * triangle_num - 1;
    triangle_neighbor(2,triangle_num) = -5;
    ledg = 2;
    ltri = triangle_num;

  else

    triangle_node(1,1) = m2;
    triangle_node(2,1) = m1;
    triangle_node(3,1) = m;
    triangle_neighbor(1,1) = -4;

    for i = 2 : triangle_num
      m1 = m2;
      m2 = i+1;
      triangle_node(1,i) = m2;
      triangle_node(2,i) = m1;
      triangle_node(3,i) = m;
      triangle_neighbor(3,i-1) = i;
      triangle_neighbor(1,i) = -3 * i - 3;
      triangle_neighbor(2,i) = i - 1;
    end

    triangle_neighbor(3,triangle_num) = -3 * triangle_num;
    triangle_neighbor(2,1) = -3 * triangle_num - 2;
    ledg = 2;
    ltri = 1;

  end
%
%  Insert the vertices one at a time from outside the convex hull,
%  determine visible boundary edges, and apply diagonal edge swaps until
%  Delaunay triangulation of vertices (so far) is obtained.
%
  top = 0;

  for i = j+1 : node_num

    m = i;
    m1 = triangle_node(ledg,ltri);

    if ( ledg <= 2 )
      m2 = triangle_node(ledg+1,ltri);
    else
      m2 = triangle_node(1,ltri);
    end

    lr = lrline ( node_xy(1,m), node_xy(2,m), node_xy(1,m1), node_xy(2,m1), ...
      node_xy(1,m2), node_xy(2,m2), 0.0 );

    if ( 0 < lr )
      rtri = ltri;
      redg = ledg;
      ltri = 0;
    else
      l = -triangle_neighbor(ledg,ltri);
      rtri = floor ( l / 3 );
      redg = mod(l,3) + 1;
    end

    [ ltri, ledg, rtri, redg ] = vbedg ( node_xy(1,m), node_xy(2,m), node_num, node_xy, ...
      triangle_num, triangle_node, triangle_neighbor, ltri, ledg, rtri, redg );

    n = triangle_num + 1;
    l = -triangle_neighbor(ledg,ltri);

    while ( 1 )

      t = floor ( l / 3 );
      e = mod ( l, 3 ) + 1;
      l = -triangle_neighbor(e,t);
      m2 = triangle_node(e,t);

      if ( e <= 2 )
        m1 = triangle_node(e+1,t);
      else
        m1 = triangle_node(1,t);
      end

      triangle_num = triangle_num + 1;
      triangle_neighbor(e,t) = triangle_num;
      triangle_node(1,triangle_num) = m1;
      triangle_node(2,triangle_num) = m2;
      triangle_node(3,triangle_num) = m;
      triangle_neighbor(1,triangle_num) = t;
      triangle_neighbor(2,triangle_num) = triangle_num - 1;
      triangle_neighbor(3,triangle_num) = triangle_num + 1;
      top = top + 1;

      if ( node_num < top )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'R8TRIS2 - Fatal error!\n' );
        fprintf ( 1, '  Stack overflow.\n' );
        error ( 'R8TRIS2 - Fatal error!' )
      end

      work(top) = triangle_num;

      if ( t == rtri & e == redg )
        break
      end

    end

    triangle_neighbor(ledg,ltri) = -3 * n - 1;
    triangle_neighbor(2,n) = -3 * triangle_num - 2;
    triangle_neighbor(3,triangle_num) = -l;
    ltri = n;
    ledg = 2;

    [ top, ltri, ledg, triangle_node, triangle_neighbor ] = swapec ( ...
      m, top, ltri, ledg, node_num, node_xy, triangle_num, triangle_node, ...
      triangle_neighbor, work );

  end
%
%  Now account for the sorting that we did.
%
  for i = 1 : 3
    for j = 1 : triangle_num
      triangle_node(i,j) = indx ( triangle_node(i,j) );
    end
  end

  indx = perm_inverse ( node_num, indx );

  node_xy = r82vec_permute ( node_num, node_xy, indx );

  return
end


text:
function order = r82row_order_type ( n, a )

%*****************************************************************************80
%
%% R82ROW_ORDER_TYPE finds if an R82ROW is (non)strictly ascending/descending.
%
%  Discussion:
%
%    An R82ROW is a (2,N) array of R8's.
%
%    The dictionary or lexicographic ordering is used.
%
%    (X1,Y1) < (X2,Y2)  <=>  X1 < X2 or ( X1 = X2 and Y1 < Y2).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries of the array.
%
%    Input, real A(2,N), the array to be checked.
%
%    Output, integer ORDER, order indicator:
%    -1, no discernable order;
%    0, all entries are equal;
%    1, ascending order;
%    2, strictly ascending order;
%    3, descending order;
%    4, strictly descending order.
%
  dim_num = 2;
%
%  Search for the first value not equal to A(1,1).
%
  i = 1;

  while ( 1 )

    i = i + 1;

    if ( n < i )
      order = 0;
      return
    end

    if ( ...
         a(1,1) <  a(1,i) || ...
       ( a(1,1) == a(1,i) && a(2,1) < a(2,i) ) )

      if ( i == 2 )
        order = 2;
      else
        order = 1;
      end

      break

    elseif ( ...
        a(1,i) <  a(1,1) || ...
      ( a(1,i) == a(1,1) && a(2,i) < a(2,1) ) )

      if ( i == 2 )
        order = 4;
      else
        order = 3;
      end

      break

    end

  end
%
%  Now we have a "direction".  Examine subsequent entries.
%
  while ( 1 )

    i = i + 1;
    if ( n < i )
      break
    end

    if ( order == 1 )

      if ( ...
          a(1,i) <  a(1,i-1) || ...
        ( a(1,i) == a(1,i-1) && a(2,i) < a(2,i-1) ) )
        order = -1;
        break
      end

    elseif ( order == 2 )

      if ( ...
          a(1,i) <  a(1,i-1) || ...
        ( a(1,i) == a(1,i-1) && a(2,i) < a(2,i-1) ) )
        order = -1;
        break
      elseif ( ...
         a(1,i) == a(1,i-1) && a(2,i) == a(2,i-1) )
        order = 1;
      end

    elseif ( order == 3 )

      if ( ...
          a(1,i-1) <  a(1,i) || ...
        ( a(1,i-1) == a(1,i) && a(2,i-1) < a(2,i) ) )
        order = -1;
        break
      end

    elseif ( order == 4 )

      if ( ...
          a(1,i-1) <  a(1,i) || ...
        ( a(1,i-1) == a(1,i) && a(2,i-1) < a(2,i) ) )
        order = -1;
        break
      elseif ( a(1,i) == a(1,i-1) && a(2,i) == a(2,i-1) )
        order = 3;
      end

    end

  end

  return
end


text:
report_this_filefun(mfilename('fullpath'));

di = km2deg(red(:,7));
az = red(:,8);
[la,lo] =reckon(az*0+r(2),az*0+r(1),di,az);
s = [lo la red(:,4) red(:,5) red(:,3)-red(:,2) red(:,6)];
cd /Seis/obelix/stefan/DEM
save redaz.dat s -ascii

author: CelsoReyes


text:
% Copyright (C) 1994-2015 John W. Eaton
%
% This file is part of Octave.
%
% Octave 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.
%
% Octave 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 Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

% -*- texinfo -*-
% @deftypefn  {Function File} {} fftfilt (@var{b}, @var{x})
% @deftypefnx {Function File} {} fftfilt (@var{b}, @var{x}, @var{n})
% Filter @var{x} with the FIR filter @var{b} using the FFT.
%
% If @var{x} is a matrix, filter each column of the matrix.
%
% Given the optional third argument, @var{n}, @code{fftfilt} uses the
% overlap-add method to filter @var{x} with @var{b} using an N-point FFT@.
% The FFT size must be an even power of 2 and must be greater than or equal to
% the length of @var{b}.  If the specified @var{n} does not meet these
% criteria, it is automatically adjusted to the nearest value that does.
%
% @seealso{filter, filter2}
% @end deftypefn

% Author: Kurt Hornik <Kurt.Hornik@wu-wien.ac.at>
% Created: 3 September 1994
% Adapted-By: jwe

function y = oc_fftfilt (b, x, n)

  % If N is not specified explicitly, we do not use the overlap-add
  % method at all because loops are really slow.  Otherwise, we only
  % ensure that the number of points in the FFT is the smallest power
  % of two larger than N and length(b).  This could result in length
  % one blocks, but if the user knows better ...

  if (nargin < 2 || nargin > 3)
    print_usage ();
  end

  transpose = (size(x,1) == 1);

  if (transpose)
    x = x.';
  end

  [r_x, c_x] = size (x);
  [r_b, c_b] = size (b);

  if (~isvector (b))
    error ('fftfilt: B must be a vector');
  end

  if (ndims (x) ~= 2)
    error ('fftfilt: X must be a 1-D or 2-D array');
  end

  l_b = r_b * c_b;
  b = reshape (b, l_b, 1);

  if (nargin == 2)
    % Use FFT with the smallest power of 2 which is >= length (x) +
    % length (b) - 1 as number of points ...
    n = 2 ^ nextpow2 (r_x + l_b - 1);
    B = fft (b, n);
    y = ifft (fft (x, n) .* B(:, ones (1, c_x)));
  else
    % Use overlap-add method ...
    if (~(isscalar (n)))
      error ('fftfilt: N has to be a scalar');
    end
    n = 2 ^ nextpow2 (max ([n, l_b]));
    L = n - l_b + 1;
    B = fft (b, n);
    B = B(:, ones (c_x,1));
    R = ceil (r_x / L);
    y = zeros (r_x, c_x);
    for r = 1:R;
      lo = (r - 1) * L + 1;
      hi = min (r * L, r_x);
      tmp = zeros (n, c_x);
      tmp(1:(hi-lo+1),:) = x(lo:hi,:);
      tmp = ifft (fft (tmp) .* B);
      hi  = min (lo+n-1, r_x);
      y(lo:hi,:) = y(lo:hi,:) + tmp(1:(hi-lo+1),:);
    end
  end

  y = y(1:r_x, :);

  % Final cleanups:

  % - If both b and x are real, y should be real.
  % - If b is real and x is imaginary, y should be imaginary.
  % - If b is imaginary and x is real, y should be imaginary.
  % - If both b and x are imaginary, y should be real.
  xisreal = all (imag (x) == 0);
  xisimag = all (real (x) == 0);

  if (all (imag (b) == 0))
    y (:,xisreal) = real (y (:,xisreal));
    y (:,xisimag) = complex (real (y (:,xisimag)) * 0, imag (y (:,xisimag)));
  elseif (all (real (b) == 0))
    y (:,xisreal) = complex (real (y (:,xisreal)) * 0, imag (y (:,xisreal)));
    y (:,xisimag) = real (y (:,xisimag));
  end

  % - If both x and b are integer in both real and imaginary
  %   components, y should be integer.
  if (~any(b - fix (b)))
    idx = find (~any(x - fix (x)));
    y (:, idx) = round (y (:, idx));
  end

  % Transpose after cleanup, otherwise rounding fails.
  if (transpose)
    y = y.';
  end

end

author: brainstorm-tools


text:
function M=axisAng2RotMat(u,theta,handed)
%%AXISANG2ROTMAT Get a rotation matrix to rotate a 3D vector an angle of
%                theta counterlockwise (right-handed) or clockwise
%                (left-handed) about an axis given by the unit vector u.
%
%INPUTS: u A unit vector representing the axis about which a vector should
%          be rotated using the rotation matrix.
%    theta The angle in radians by which a vector is to be rotated is to be
%          rotated about the axis u when multiplied by the rotation matrix.
%   handed The handedness of the rotation angle. If omitted, it is assumed
%          that the rotation is right-handed (the standard). Possible
%          values are
%          'right' The default if omitted. The rotation is right-handed.
%           'left' The rotation is left-handed. The rotation angle is
%                  counterclockwise when one is looking in the same
%                  direction that the rotation axis points.
%
%OUTPUTS: M The rotation matrix such that M*v rotates the vector v by the
%           desired rotation.
%
%To rotate a 3D vector x an angle of theta about the axis u, simply
%evaluate M*x.
%
%The rotation matrix is obtained by obtaining the unit quaternion for the
%rotation and turning it into the corresponding rotation matrix as
%described in  [1].
%
%REFERENCES:
%[1] David F. Crouse , "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%December 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(handed))
    handed='right';
end

q=axisAng2Quat(u,theta);
M=quat2RotMat(q,handed);

end

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


text:
function [ sout, slist ] = cvx_solver( sname )

%CVX_SOLVER    CVX solver selection.
%   CVX_SOLVER <solver_name> or CVX_SOLVER('<solver_name>')
%   selects the named solver the CVX uses to solve models. The solver name
%   is case-insensitive; so, for example, both 'SeDuMi' and 'sedumi' will
%   select the same solver.
%
%   When CVX is first installed, the solver SDPT3 is selected as the
%   default. For most problems, this will be a good choice; nevertheless,
%   no solver is perfect, so if you encounter issues you may wish to
%   experiment with other solvers.
%
%   There are two ways to use the CVX_SOLVER command. If you use it within
%   a model---that is, between the statements CVX_BEGIN and CVX_END---then
%   the new solver selection will apply only to that particular model. For
%   instance, if the default solver is SDPT3, then the following structure
%   will solve a single model using SeDuMi instead:
%       cvx_begin
%           cvx_solver sedumi
%           variables ...
%           ...
%       cvx_end
%   On the other hand, if CVX_SOLVER is called *outside* of a model, then
%   the change will apply for all subsequent models, or until you call
%   CVX_SOLVER once again.
%
%   [ SOLVER, SOLVER_LIST ] = CVX_SOLVER returns the name of the current
%   solver, and a cell array containing the names of all available choices.
%
%   Calling CVX_SOLVER with no input or output arguments produces a listing
%   of the solvers that CVX currently recognized, and an indication of the
%   current solver selection and/or the default.

global cvx___
cvx_global
if nargin,
    if isempty( sname ),
        sname = 'default';
    elseif ~ischar( sname ) || size( sname, 1 ) ~= 1,
        cvx_throw( 'Argument must be a string.' );
    end
    try
        snumber = cvx___.solvers.map.(lower(sname));
    catch
        cvx_throw( 'Unknown, unusable, or missing solver: %s', sname );
    end
    if ~isempty( cvx___.solvers.list(snumber).error ),
        cvx_throw( 'Solver unusable due to prior errors: %s', sname );
    end
    if isempty( cvx___.problems ),
        cvx___.solvers.selected = snumber;
    elseif ~isa( evalin( 'caller', 'cvx_problem', '[]' ), 'cvxprob' ),
        cvx_throw( 'The global CVX solver selection cannot be changed while a model is being constructed.' );
    elseif cvx___.problems(end).solver ~= snumber,
        cvx___.problems(end).solver = snumber;
        cvx___.problems(end).settings = cvx___.solvers.list.settings;
    end
    if cvx___.solvers.active,
        cvx_setspath;
    end
elseif nargout == 0,
    solvers = cvx___.solvers;
    snames = solvers.names;
    statvec = [ 0, solvers.map.default, solvers.active ];
    statstr = { 'selected', 'default', 'active' };
    if ~isempty( cvx___.problems ),
        statvec(1) = cvx___.problems(end).solver;
    else
        statvec(1) = solvers.selected;
    end
    fprintf( '\n' );
    dash = '-';
    solvers = solvers.list;
    nsolv = length( solvers );
    lens = [4,6,7,8];
    for k = 1 : nsolv,
        if isempty( solvers(k).error ),
            nstat = statstr(k==statvec);
            if ~isempty(nstat),
                nstat = sprintf( '%s,', nstat{:} );
                nstat = nstat(1:end-1);
            end
        else
            nstat = 'disabled';
        end
        lens = max( lens, [ length(snames{k}), length(nstat), length(solvers(k).version), length(solvers(k).location) ] );
    end
    fmt = sprintf( '   %%-%ds   %%-%ds   %%-%ds   %%s\\n', lens(1), lens(2), lens(3) );
    fprintf( fmt, 'Name', 'Status', 'Version', 'Location' ); %#ok
    fprintf( '%s\n', dash(ones(1,sum(lens)+15)) );
    for k = 1 : nsolv,
        if isempty( solvers(k).error ),
            nstat = statstr(k==statvec);
            if ~isempty(nstat),
                nstat = sprintf( '%s,', nstat{:} );
                nstat = nstat(1:end-1);
            else
                nstat = '';
            end
        else
            nstat = 'disabled';
        end
        fprintf( fmt, snames{k}, nstat, solvers(k).version, solvers(k).location );
    end
    fprintf( '\n' );
end
if nargout > 0,
    sout = cvx___.solvers.list(cvx___.solvers.selected).name;
    slist = cvx___.solvers.names;
end

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

author: yu-jiang

repo: radpbook

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


text:
% HEALPIX Border

clear

n = 2;
n = 4;
n = 4;

% sampling interval on the unit sphere
interval = pi / 90;

hold on
axis equal

% polar cap area

for k = 1:n
    start_phi = pi * k / (2 * n);
    end_phi = pi / 2;
    S = HealpixBorderLinePC(n, k, interval, start_phi, end_phi, 0);
    for m = 0:3
        % Northern Hemisphere
        C = SphToCart([S(:, 1), S(:, 2) + pi * m / 2]);
        plot3(C(:,1),C(:,2),C(:,3), 'k-')
        % Southern Hemisphere
        C = SphToCart([pi - S(:, 1), S(:, 2) + pi * m / 2]);
        plot3(C(:,1),C(:,2),C(:,3), 'k-')
    end

    start_phi = 0;
    end_phi = pi / 2 - pi * k / (2 * n);
    S = HealpixBorderLinePC(n, k, interval, start_phi, end_phi, 1);
    for m = 0:3
        % Northern Hemisphere
        C = SphToCart([S(:, 1), S(:, 2) + pi * m / 2]);
        plot3(C(:,1),C(:,2),C(:,3), 'k-')
        % Southern Hemisphere
        C = SphToCart([pi - S(:, 1), S(:, 2) + pi * m / 2]);
        plot3(C(:,1),C(:,2),C(:,3), 'k-')
    end
end

for k = 0:3
    % Northern Hemisphere
    PHI = 0:interval:acos(2 / 3);
    if PHI(end) ~= acos(2 / 3)  % if the end edge point was not included
        PHI = [PHI acos(2 / 3)];    % add the end edge
    end
    S = [PHI.', ones(size(PHI, 2), 1) * k * pi / 2];
    C = SphToCart(S);
    plot3(C(:,1),C(:,2),C(:,3), 'k-')

    % Southern Hemisphere
    PHI = acos(-2 / 3):interval:pi;
    if PHI(end) ~= pi  % if the end edge point was not included
        PHI = [PHI pi];    % add the end edge
    end
    S = [PHI.', ones(size(PHI, 2), 1) * k * pi / 2];
    C = SphToCart(S);
    plot3(C(:,1),C(:,2),C(:,3), 'k-')
end

% equatorial belt area

for k = (-3 * n):(n - 1)
    start_theta = acos(-2 / 3);
    start_phi = (-4 / 3 + 4 * k / (3 * n)) * 3 * pi / 8;
    end_phi = 4 * k / (3 * n) * 3 * pi / 8;
    S = HealpixBorderLineEB(n, k, interval, start_phi, end_phi, start_theta);
    C = SphToCart(S);
    plot3(C(:,1),C(:,2),C(:,3), 'k-')

    start_theta = acos(2 / 3);
    temp = -start_phi;
    start_phi = -end_phi;
    end_phi = temp;
    S = HealpixBorderLineEB(n, k, interval, start_phi, end_phi, start_theta);
    C = SphToCart(S);
    plot3(C(:,1),C(:,2),C(:,3), 'k-')
end

k = 5;
n = 2^k-1;
[x,y,z] = sphere(n);
c = zeros(32);
scale = 0.995;
surf(x*scale,y*scale,z*scale,c);
colormap([1  1  1])
shading flat
view(20,30);

author: mtex-toolbox

repo: mtex

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


text:
function r8_erf_inverse_test ( )

%*****************************************************************************80
%
%% R8_ERF_INVERSE_TEST tests R8_ERF_INVERSE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 August 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_ERF_INVERSE_TEST:\n' );
  fprintf ( 1, '  R8_ERF_INVERSE inverts the error function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    FX            X    R8_ERF_INVERSE(FX)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x1, fx ] = erf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    x2 = r8_erf_inverse ( fx );

    fprintf ( 1, '  %6f  %12f  %12f\n', fx, x1, x2 );

  end

  return
end


text:
function sgmgg_print ( gg_ma, gg_mi, gg_mo, gg_na, gg_nd, gg_ni, gg_no, ...
  gg_a, gg_b, gg_f, gg_g, gg_i, gg_o, gg_s )

%*****************************************************************************80
%
%% SGMGG_PRINT prints out an SGMGG data structure.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer GG_MA, the maximum dimension for GG_A.
%
%    Input, integer GG_MI, the maximum dimension for GG_I.
%
%    Input, integer GG_MO, the maximum dimension for GG_O.
%
%    Input, integer GG_NA, the current dimension for GG_A.
%
%    Input, integer GG_ND, the spatial dimension.
%
%    Input, integer GG_NI, the current dimension for GG_I.
%
%    Input, integer GG_NO, the current dimension for GG_O.
%
%    Input, integer GG_A(GG_MA), the active indices.
%
%    Input, integer GG_B(GG_ND,GG_MI), the forward neighbors.
%
%    Input, integer GG_F(GG_ND,GG_MI), the backward neighbors.
%
%    Input, real GG_G(GG_MA), the error estimators.
%
%    Input, integer GG_I(GG_ND,GG_MI), the index set.
%
%    Input, integer GG_O(GG_MO), the old indices.
%
%    Input, integer GG_S(GG_MI), 0 if index I is old, 1
%    if it is active.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GG DATA STRUCTURE:\n' );
  fprintf ( 1, '  ND = %d\n', gg_nd );
  fprintf ( 1, '  NI = %d\n', gg_ni );
  fprintf ( 1, '  NO = %d\n', gg_no );
  fprintf ( 1, '  NA = %d\n', gg_na );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Indices:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I       A/O    G(I)        Index values\n' );
  fprintf ( 1, '\n' );
  for i = 1 : gg_ni
    if ( gg_s(i) == 0 )
      s = 'o';
    elseif ( gg_s(i) == 1 )
      s = 'a';
    else
      s = '?';
    end
    fprintf ( 1, '  %4d        %c  %14g', i, s, gg_g(i) );
    for j = 1 : gg_nd
      fprintf ( 1, '  %4d', gg_i(j,i) );
    end
    fprintf ( 1, '\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Backward neighbors:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : gg_ni
    fprintf ( 1, '  %4d      ', i );
    for j = 1 : gg_nd
      fprintf ( 1, '  %4d', gg_b(j,i) );
    end
    fprintf ( 1, '\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Forward neighbors:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : gg_ni
    fprintf ( 1, '  %4d      ', i );
    for j = 1 : gg_nd
      fprintf ( 1, '  %4d', gg_f(j,i) );
    end
    fprintf ( 1, '\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Active Heap:\n' );
  fprintf ( 1, '     I     A      G\n' );
  fprintf ( 1, '\n' );

  for i = 1 : gg_na
    fprintf ( 1, '  %4d  %4d  %14g\n', i, gg_a(i), gg_g(gg_a(i)) );
  end

  return
end


text:
function [p sres sres_ns] = ResidScan(res, FWHM)
% Calculates P(M>=t) where M is the max value of the smoothed residuals.
% In this implementation the residuals are smoothed using a Gaussian
% kernel.
%
% :Usage:
% ::
%
%     function [p sres sres_ns] = ResidScan(res, FWHM)
%
% :Inputs:
%
%   **res:**
%        residual time course
%
%   **FWHM:**
%        Full Width Half Maximum (in time units)
%
% :Outputs:
%
%   **p:**
%        pvalues
%
%   **sres:**
%        smoothed residuals
%
%   **sres_ns:**
%        smoothed residuals (non standardized)
%
% ..
%    By Martin Lindquist & Ji-Meng Loh, July 2007
%
%    Edited by ML on 10/02/09
% ..

res_ns = res;
res = res./std(res);
len = length(res);

% Create Gaussian Kernel
sig = ceil(FWHM/(2*sqrt(2*log(2))));
klen = 3*sig;
kern = normpdf((-klen:klen),0,sig);
kern = kern./sqrt(sum(kern.^2));

% Convolve
x = conv(res,kern);
sres = x((klen + 1):(end-klen));

x = conv(res_ns,kern/sum(kern));
sres_ns = x((klen + 1):(end-klen));

% Find Max value
[a,location] = max(abs(sres));

% Find p-values using Gaussian Random Field theory
z = Euler_p(1, a, len, FWHM);
z = 2*z;        %Two-sided test
p = min(1, z);

end

% END MAIN FUNCTION

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Subfunctions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function pval = Euler_p(myDim, value, N, fwhm)
% function z = Euler_p(myDim, value, N, fwhm)
%
% Finds the p value using the expected Euler characteristic.
%
% This function returns P(M \ge value) using the approximation
% \sum_{d=0}^D R_d(V) \rho_d(value) following Worsley et al's "A Unified
% Statistical Approach for Determining Significant Signals in Images of
% Cerebral Activation".
%
% INPUTS:
%
% myDim - the number of dimensions in the data
% value - the value of the maximum.
% N     - the number of (time) points in that 1 dimension
% fwhm  - the full width half maximum
%
% OUTPUTS:
%
% pval - the p-value

% NOTE: CURRENTLY THIS FUNCTION IS ONLY IMPLEMENTED FOR THE 1D CASE

  % Constants
  myfactor = 4*log(2);
  pi2 = 2*pi;
  exptsq = exp(-(value^2)/2);

  % Euler Characteristc Densties
  rho = zeros(5,1);
  rho(1) = 1-normcdf(value);
  rho(2) = myfactor^(0.5)*exptsq/pi2;
  rho(3) = myfactor * exptsq * value / (pi2 ^ (1.5));
  rho(4) = myfactor ^ (1.5) * exptsq * (value^2-1) / (pi2 ^2);
  rho(5) = myfactor ^2 * exptsq * (value^3-3*value) / (pi2 ^ (5/2));

  % Resel Count
  R0 = 1;
  R1 = N/fwhm;

  % P-value
  pval = R0 * rho(1) + R1 * rho(2);

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% driver for statistical and systems heterogeneity experiments

%% load dataset
datarepo = 'data/'; % location of data folder
dataset = 'small'; % small test dataset
load([datarepo dataset]); % load data

%% set parameters
addpath('opt/'); addpath('util/'); % add helper functions
training_percent = 0.75; % percentage of data for training
[Xtrain, Ytrain, Xtest, Ytest] = split_data(X, Y, training_percent);
opts.obj='C'; % classification
opts.avg = 1; % compute average error
opts.sys_het = 0; % run systems (1) or stats heterogeneity exps (0)
opts.top = 1.0; % highest number of rounds
opts.bottom = 0.1; % lowest number of rounds
lambda = 1e-4;

%% mocha [need to tune sdca_frac]
opts.mocha_outer_iters = 1;
opts.mocha_inner_iters = 2000;
opts.mocha_sdca_frac = 1.0;
opts.w_update = 1; % just do a single w-update
[rmse_mocha, primal_mocha, dual_mocha] = run_mocha(Xtrain, Ytrain, Xtest, Ytest, lambda, opts);

%% cocoa [need to tune sdca_frac]
opts.cocoa_outer_iters = 1;
opts.cocoa_inner_iters = 500;
opts.theta = 0.5;
[rmse_cocoa, primal_cocoa, dual_cocoa, max_its] = run_cocoa(Xtrain, Ytrain, Xtest, Ytest, lambda, opts);

%% mbsdca [need to tune sdca_frac, scaling]
opts.mbsdca_outer_iters = 1;
opts.mbsdca_inner_iters = 5000;
opts.mbsdca_sdca_frac = 0.5;
opts.mbsdca_scaling = 10;
[rmse_mbsdca, primal_mbsdca, dual_mbsdca] = run_mbsdca(Xtrain, Ytrain, Xtest, Ytest, lambda, opts);

%% mbsgd [need to tune sgd_frac, scaling]
opts.mbsgd_outer_iters = 1;
opts.mbsgd_inner_iters = 5000;
opts.mbsgd_sgd_frac = 0.5;
opts.mbsgd_scaling = 0.1;
[rmse_mbsgd, primal_mbsgd] = run_mbsgd(Xtrain, Ytrain, Xtest, Ytest, lambda, opts);

if(opts.sys_het)
    %% plot systems heterogeneity
    % note: ensure methods have reached global optimal or enter manually
    optimal = min([primal_mocha; primal_cocoa; primal_mbsdca; primal_mbsgd]);

    %% calculate estimated time
    comm_cost = 100; % communication cost: Wifi=10, LTE=100, 3G=1000
    train_n = 0;
    for t=1:length(Xtrain)
         train_n = train_n + size(Xtrain{t},1);
    end

    %% calculate time based on flops and communication cost
    local_mocha_time = 8 * opts.top * train_n + comm_cost;
    local_mbsdca_time = 6 * opts.top * train_n + comm_cost;
    local_mbsgd_time = 4 * opts.top * train_n + comm_cost;
    mocha_time = 1:local_mocha_time:local_mocha_time*length(primal_mocha);
    cocoa_time = cumsum(max_its .* (8 * length(Xtrain)) + comm_cost);
    mbsdca_time = 1:local_mbsdca_time:local_mbsdca_time*length(primal_mbsdca);
    mbsgd_time = 1:local_mbsgd_time:local_mbsgd_time*length(primal_mbsgd);

    %% plot results
    figure;
    step = 100;
    semilogy(mocha_time(1:step:end), primal_mocha(1:step:end) - optimal, 'LineWidth', 6)
    hold on;
    semilogy(cocoa_time(1:step:end), primal_cocoa(1:step:end) - optimal, 'LineWidth', 6)
    hold on;
    semilogy(mbsdca_time(1:step:end), primal_mbsdca(1:step:end) - optimal, 'LineWidth', 6)
    hold on;
    semilogy(mbsgd_time(1:step:end), primal_mbsgd(1:step:end) - optimal, 'LineWidth', 6)
    title([dataset ': Systems Heterogeneity'])
    xlabel('Estimated Time')
    ylabel('Primal Sub-Optimality')
    set(gca, 'fontsize', 16)
    legend({'MOCHA', 'CoCoA', 'Mb-SDCA', 'Mb-SGD'})
    axis([0 8000000 .001 100]) % set manually

else
    %% plot statistical heterogeneity
    % note: ensure methods have reached global optimal or enter manually
    optimal = min([primal_mocha; primal_cocoa; primal_mbsdca; primal_mbsgd]);

    %% calculate time based on flops and communication cost
    train_n = 0;
    for t=1:length(Xtrain)
         train_n = train_n + size(Xtrain{t},1);
    end
    comm_cost = 100; % communication cost: Wifi=10, LTE=100, 3G=1000
    local_mocha_time = 8 * opts.mocha_sdca_frac * train_n + comm_cost;
    local_mbsdca_time = 6 * opts.mbsdca_sdca_frac * train_n + comm_cost;
    local_mbsgd_time = 4 * opts.mbsgd_sgd_frac * train_n + comm_cost;
    cocoa_time = cumsum(max_its .* (8 * length(Xtrain)) + comm_cost);
    mocha_time = 1:local_mocha_time:local_mocha_time*length(primal_mocha);
    mbsdca_time = 1:local_mbsdca_time:local_mbsdca_time*length(primal_mbsdca);
    mbsgd_time = 1:local_mbsgd_time:local_mbsgd_time*length(primal_mbsgd);

    %% plot results
    figure;
    step = 100;
    semilogy(mocha_time(1:step:end), primal_mocha(1:step:end) - optimal, 'LineWidth', 6)
    hold on;
    semilogy(cocoa_time, primal_cocoa - optimal, 'LineWidth', 6)
    hold on;
    semilogy(mbsdca_time(1:step:end), primal_mbsdca(1:step:end) - optimal, 'LineWidth', 6)
    hold on;
    semilogy(mbsgd_time(1:step:end), primal_mbsgd(1:step:end) - optimal, 'LineWidth', 6)
    title([dataset ': Statistical Heterogeneity'])
    xlabel('Estimated Time')
    ylabel('Primal Sub-Optimality')
    set(gca, 'fontsize', 16)
    legend({'MOCHA', 'CoCoA', 'Mb-SDCA', 'Mb-SGD'})
    axis([0 7000000 .001 100]) % set manually

end

author: gingsmith

repo: fmtl

save_path: github-repos/MATLAB/gingsmith-fmtl