text:
function X = create2dbf(g1,g2,N)

n1 = length(g1);
n2 = length(g2);

X1 = zeros(n1,N,2);
for i = 1:N
    X1(:,i,1) = cos((i-1).*pi.*[0:(n1-1)]./(n1-1));
    X1(:,i,2) = sin((i-1).*pi.*[0:(n1-1)]./(n1-1));
end
X2 = zeros(n2,N,2);
for i = 1:N
    X2(:,i,1) = cos((i-1).*pi.*[0:(n2-1)]./(n2-1));
    X2(:,i,2) = sin((i-1).*pi.*[0:(n2-1)]./(n2-1));
end

X = zeros(n1,n2,N^2);
ij = 1;
for i=1:N
    for j=1:N
        for i1=1:2
            for i2=1:2
                X(:,:,ij) = repmat(X1(:,i,i1),1,n2).*repmat(X2(:,j,i2)',n1,1);
                ij = ij +1;
            end
        end
    end
end

author: MBB-team

repo: VBA-toolbox


text:
function testvlfeatsift()
 dir='\\File.ceegs.ohio-state.edu\SPIN\UAV\August 08, 2013\Nikon D800\109ND800\';
fn1=[dir,'DSC_1023.jpg'];
fn2=[dir,'DSC_1024.jpg'];
im1=imread(fn1);
im2=imread(fn2);
Nimg1 = gaussdownsample(im1,4,[5 5],1.6);
Nimg2 = gaussdownsample(im2,4,[5 5],1.6);

figure; imshow(Nimg1);

sift_mosaic(Nimg1, Nimg2);

function Nimg = gaussdownsample(img,levels,gsize,sigma);
%
%  Function to generate a gaussian-pyramid for the given input image
% image is resize to width/2^(levels-1)*height/2^(levels-1)
% Input:
%        img: input image-matrix grayscale
%        levels:  number of levels of the pyramid, we only return the level
%        of the least size
%        gsize: size of the gaussian kernel [w h] ([5 5] normally provides a smooth output)
%        sigma:  sigma for gaussian kernel
% Output:
%        Nimg:  the smoothed image of the smallest size
% Usage:
%      im = imread('cameraman.tif');
%      Nimg = Gscale(im,3,[5 5],1.6);
%      figure; imshow(Nimg);
Nimg=img;
if(levels==1)
    return;
end
%guassian filter  with a sigma=1.6
g = fspecial('gaussian',gsize,sigma);
%pyramid
for i = 2:levels
        %perform guassian filtering
        im = imfilter(Nimg,g,'conv');
        %perform downsampling (horizontal)
        im1 = im(:,1:2:size(Nimg,2),:);
        %vertical
        im2 = im1(1:2:size(Nimg,1),:,:);
        %store it in a struct format
        Nimg = im2;
end

%End

repo: ekfmonoslam


text:
function lbeta_test ( )

%*****************************************************************************80
%
%% LBETA_TEST tests R4_LBETA and R8_LBETA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LBETA_TEST:\n' );
  fprintf ( 1, '  Test BETA_LOG_VALUES, R4_LBETA, R8_LBETA.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             A               B     LBETA(A,B)\n' );
  fprintf ( 1, '                                R4_LBETA(A,B)         Diff\n' );
  fprintf ( 1, '                                R8_LBETA(A,B)         Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, a, b, fx1 ] = beta_log_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_lbeta ( single ( a ), single ( b ) );
    fx3 = r8_lbeta ( a, b );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.4g  %14.6g\n', a, b, fx1 );
    fprintf ( 1, '                                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ know, x ] = p18_sol ( n )

%*****************************************************************************80
%
%% P18_SOL returns the solution for problem 18.
%
%  Discussion:
%
%    The solution values are taken from Brent.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the problem.  This value
%    is only needed for those problems with variable N.
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Output, real X(N), the solution, if known.
%
  if ( n == 2 )
    know = 1;
    x = [ 0.2113249, 0.7886751 ]';
  elseif ( n == 4 )
    know = 1;
    x = [ 0.1026728, 0.4062037, 0.5937963, 0.8973272 ]';
  elseif ( n == 6 )
    know = 1;
    x = [ 0.066877, 0.288741, 0.366682, 0.633318, ...
      0.711259, 0.933123 ]';
  elseif ( n == 8 )
    know = 1;
    x = [ 0.043153, 0.193091, 0.266329, 0.500000, ...
      0.500000, 0.733671, 0.806910, 0.956847 ]';
  else
    know = 0;
    x = zeros ( n, 1 );
  end

  return
end


text:
function [tSnrImageArray, fileTsnrArray, ...
            tSnrRatioImageArray, fileTsnrRatioArray ] = ...
    tapas_physio_compute_tsnr_gains(physio, SPM, indexContrastForSnrRatio, ...
    namesPhysContrasts)
% Computes tSNR gains through physiological noise correction for all
% confound regressor sets modelled, after estimation of the SPM-GLM
%
%   [fileTsnrGainArray, fileTsnrArray] = ...
%           tapas_physio_compute_tsnr_gains(physio, SPM, doSave);
%
%
% IN
%   physio      physio-structure (or filename with structure), after estimating multiple_regressors.txt
%   SPM         SPM variable (or filename with saved structure)
%               after estimation of general linear model (GLM)
%   indexContrastForSnrRatio
%               contrast id for comparison to get relative tSNR values
%               default: 0 (raw time series after filtering/pre-whitening)
% OUT
%   fileTsnrGainArray
%               cell(nPhysioSets,1) of nii-filenames holding tSNR gain
%               images for all physiological regressor sets (in percent?)
%
%   fileTsnrArray
%               cell(nPhysioSets+1,1) of nii-filenames for tSNR images of
%               for all physiological regressor sets
%               and, as last element, raw tSNR (after preprocessing)
% EXAMPLE
%   tapas_physio_compute_tsnr_gains
%
%   See also tapas_physio_compute_tsnr_spm spm_write_residuals

% Author: Lars Kasper
% Created: 2015-07-03
% Copyright (C) 2015 TNU, Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public
% License (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   0. Check variable structure, cast for convenience
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 3
    indexContrastForSnrRatio = 0; % contrast for comparison
end

if nargin < 4
    namesPhysContrasts = tapas_physio_get_contrast_names_default();
end

% load physio-variable, if filename given
if ~isstruct(physio)
    % load SPM variable from file
    if iscell(physio)
        filePhysio = physio{1};
    else
        filePhysio = physio;
    end
    load(filePhysio, 'physio');
end

% load SPM-variable, if filename given
if ~isstruct(SPM)
    % load SPM variable from file
    if iscell(SPM)
        fileSpm = SPM{1};
    else
        fileSpm = SPM;
    end
    load(fileSpm, 'SPM');
end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Determine which physiological contrasts could be created with existing
% model and create the missing ones
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

SPM = tapas_physio_create_missing_physio_contrasts(SPM, physio.model, ...
    namesPhysContrasts);

nContrasts = numel(namesPhysContrasts);

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. Compute tSNR maps via tapas_physio_compute_tsnr_spm, as follows
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tSnrImageArray = cell(nContrasts,1);
fileTsnrArray = cell(nContrasts,1);
tSnrRatioImageArray = cell(nContrasts,1);
fileTsnrRatioArray = cell(nContrasts,1);

for c = 1:nContrasts
    % First try, whether inverted contrast already exists, then no
    % additional computation has to be performed!
    indContrast = tapas_physio_check_get_xcon_index(SPM, ...
        ['All but: ' namesPhysContrasts{c}]);
    doComputeInvertedContrast = indContrast == 0; % not found, do invert
    doSaveNewContrasts = doComputeInvertedContrast; % and save contrasts for later
    if doComputeInvertedContrast % find contrast before inversion
        indContrast = tapas_physio_check_get_xcon_index(SPM, ...
            namesPhysContrasts{c});
    end

    if indContrast > 0
        % if contrast exist, compute tSNR and save it!
        [tSnrImageArray{c}, fileTsnrArray{c}, ...
            tSnrRatioImageArray{c}, fileTsnrRatioArray{c}] = ...
            tapas_physio_compute_tsnr_spm(SPM, indContrast, ...
            indexContrastForSnrRatio, doComputeInvertedContrast, doSaveNewContrasts);

        % rename files for better human readability
        filenameSfx = regexprep(namesPhysContrasts{c}, ' ', '_');
        fileTsnrOld = fileTsnrArray{c};
        fileTsnrArray{c} = regexprep(fileTsnrArray{c}, ...
            sprintf('con%04d', indContrast), filenameSfx);

        fileTsnrRatioOld = fileTsnrRatioArray{c};
        fileTsnrRatioArray{c} = regexprep(fileTsnrRatioArray{c}, ...
            sprintf('con%04dvs%04d', indContrast, indexContrastForSnrRatio), ...
            [filenameSfx '_vs_Raw']);

        movefile(fileTsnrOld, fileTsnrArray{c});
        movefile(fileTsnrRatioOld, fileTsnrRatioArray{c});
    end
end

author: translationalneuromodeling

repo: tapas


text:
% HIST2 - draw superimposed histograms
%
% Usage:
%   >> hist2(data1, data2, bins);
%
% Inputs:
%   data1   - data to plot first process
%   data2   - data to plot second process
%
% Optional inputs:
%   bins    - vector of bin center
%
% Author: Arnaud Delorme (SCCN, UCSD)

% Copyright (C) 2003 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% 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.

% draw superimposed histograms
% ---------------
function hist2(data1, data2, bins);

if nargin < 1
    help hist2;
    return;
end
if nargin < 3
    bins = linspace(min(min(data1), min(data2)), max(max(data1), max(data2)), 100);
elseif length(bins) == 1
    bins = linspace(min(min(data1), min(data2)), max(max(data1), max(data2)), bins);
end

hist(data1, bins);
hold on; hist(data2, bins);
%figure; hist( [ measure{:,5} ], 20);
%hold on; hist([ measure{:,2} ], 20);
c = get(gca, 'children');
numfaces = size(get(c(1), 'Vertices'),1);
set(c(1), 'FaceVertexCData', repmat([1 0 0], [numfaces 1]), 'Cdatamapping', 'direct', 'facealpha', 0.5, 'edgecolor', 'none');
numfaces = size(get(c(2), 'Vertices'),1);
set(c(2), 'FaceVertexCData', repmat([0 0 1], [numfaces 1]), 'Cdatamapping', 'direct', 'facealpha', 0.5, 'edgecolor', 'none');
ylabel('Number of values');
xlim([bins(1) bins(end)]);

yl = ylim;
xl = xlim;
line([xl(1) xl(1)]+(xl(2)-xl(1))/2000, yl, 'color', 'k');
line(xl, [yl(1) yl(1)]+(yl(2)-yl(1))/2000, 'color', 'k');

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [ xtab, weight ] = laguerre_compute ( norder, alpha )

%*****************************************************************************80
%
%% LAGUERRE_COMPUTE computes a Gauss-Laguerre quadrature rule.
%
%  Discussion:
%
%    In the simplest case, ALPHA is 0, and we are approximating the
%    integral from 0 to +oo of EXP(-X) * F(X).  When this is so,
%    it is easy to modify the rule to approximate the integral from
%    A to +oo as well.
%
%    If ALPHA is nonzero, then there is no simple way to extend the
%    rule to approximate the integral from A to +oo.  The simplest
%    procedures would be to approximate the integral from 0 to A.
%
%    The integration interval is [ A, +oo ) or [ 0, +oo ).
%
%    The weight function w(x) = EXP ( - X ) or EXP ( - X ) * X^ALPHA.
%
%
%    If the integral to approximate is:
%
%        Integral ( A <= X < +oo ) EXP ( - X ) * F(X) dX
%      or
%        Integral ( 0 <= X < +oo ) EXP ( - X ) * X^ALPHA * F(X) dX
%
%    then the quadrature rule is:
%
%      EXP ( - A ) * Sum ( 1 <= I <= ORDER ) WEIGHT(I) * F ( A+XTAB(I) )
%    or
%      sum ( 1 <= I <= ORDER ) WEIGHT(I) * F ( XTAB(I) )
%
%
%    If the integral to approximate is:
%
%        Integral ( A <= X < +oo ) F(X) dX
%      or
%        Integral ( 0 <= X < +oo ) X^ALPHA * F(X) dX
%
%    then the quadrature rule is:
%
%      EXP ( - A ) * sum ( 1 <= I <= ORDER )
%        WEIGHT(I) * EXP(A+XTAB(I)) * F ( A+XTAB(I) )
%    or
%      sum ( 1 <= I <= ORDER ) WEIGHT(I) * EXP(XTAB(I)) * F ( XTAB(I) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 October 2005
%
%  Author:
%
%    Original FORTRAN77 version by Arthur Stroud, Don Secrest.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%  Parameters:
%
%    Input, integer NORDER, the order of the quadrature rule to be computed.
%    NORDER must be at least 1.
%
%    Input, real ALPHA, the exponent of the X factor.
%    Set ALPHA = 0.0 for the simplest rule.
%    ALPHA must be nonnegative.
%
%    Output, real XTAB(NORDER), the Gauss-Laguerre abscissas.
%
%    Output, real WEIGHT(NORDER), the Gauss-Laguerre weights.
%

%
%  Set the recursion coefficients.
%
  for i = 1 : norder
    b(i) = ( alpha + 2 * i - 1 );
  end

  for i = 1 : norder
    c(i) = ( i - 1 ) * ( alpha + i - 1 );
  end

  cc = gamma ( alpha + 1.0 ) * prod ( c(2:norder) );

  for i = 1 : norder
%
%  Compute an estimate for the root.
%
    if ( i == 1 )

      x = ( 1.0 + alpha ) * ( 3.0 + 0.92 * alpha ) ...
        / ( 1.0 + 2.4 * norder + 1.8 * alpha );

    elseif ( i == 2 )

      x = x + ( 15.0 + 6.25 * alpha ) ...
        / ( 1.0 + 0.9 * alpha + 2.5 * norder );

    else

      r1 = ( 1.0 + 2.55 * ( i - 2 ) ) / ( 1.9 * ( i - 2 ) );

      r2 = 1.26 * ( i - 2 ) * alpha / ( 1.0 + 3.5 * ( i - 2 ) );

      ratio = ( r1 + r2 ) / ( 1.0 + 0.3 * alpha );

      x = x + ratio * ( x - xtab(i-2) );

    end
%
%  Use iteration to find the root.
%
    [ x, dp2, p1 ] = laguerre_root ( x, norder, alpha, b, c );
%
%  Set the abscissa and weight.
%
    xtab(i) = x;
    weight(i) = ( cc / dp2 ) / p1;

  end

  return
end


text:
function [output] = volumefillholes(input, along)

% VOLUMEFILLHOLES is a helper function for segmentations
%
% See also VOLUMETHRESHOLD, VOLUMESMOOTH

% ensure that SPM is available, needed for spm_bwlabel
hasspm = ft_hastoolbox('spm8up', 3) || ft_hastoolbox('spm2', 1);

if nargin<2
  inflate = false(size(input)+2);                   % grow the edges along each dimension
  inflate(2:end-1, 2:end-1, 2:end-1) = (input~=0);  % insert the original volume
  [lab, num] = spm_bwlabel(double(~inflate), 18);   % note that 18 is consistent with imfill, 26 is not
  if num>1
    inflate(lab~=lab(1)) = true;
    output  = inflate(2:end-1, 2:end-1, 2:end-1);   % trim the edges
  else
    output = input;
  end

else
  output = input;
  dim    = size(input);
  switch along
    case 1
      for i=1:dim(1)
        slice=reshape(input(i,:,:),dim([2 3]));
        im = imfill(slice,8,'holes');
        output(i,:,:) = im;
      end

    case 2
      for i=1:dim(2)
        slice=reshape(input(:,i,:),dim([1 3]));
        im = imfill(slice,8,'holes');
        output(:,i,:) = im;
      end

    case 3
      for i=1:dim(3)
        slice=reshape(input(:,:,i),dim([1 2]));
        im = imfill(slice,8,'holes');
        output(:,:,i) = im;
      end

    otherwise
      ft_error('invalid dimension along which to slice the volume');
  end % switch
end % if nargin

author: spm


text:
% Color code for each image:

if ~exist('n_ima')|~exist('fc'),
    fprintf(1,'No calibration data available.\n');
    return;
end;

check_active_images;

if n_ima ~=0,
if ~exist(['ex_' num2str(ind_active(1)) ]),
    fprintf(1,'Need to calibrate before analysing reprojection error. Maybe need to load Calib_Results.mat file.\n');
    return;
end;
end;

%if ~exist('no_grid'),
no_grid = 0;
%end;

colors = 'brgkcm';

figure(5);

for kk = 1:n_ima,
    if exist(['y_' num2str(kk)]),
        if active_images(kk) & eval(['~isnan(y_' num2str(kk) '(1,1))']),

            if ~no_grid,
                eval(['XX_kk = X_' num2str(kk) ';']);
                N_kk = size(XX_kk,2);

                if ~exist(['n_sq_x_' num2str(kk)]),
                    no_grid = 1;
                end;

                if ~no_grid,
                    eval(['n_sq_x = n_sq_x_' num2str(kk) ';']);
                    eval(['n_sq_y = n_sq_y_' num2str(kk) ';']);
                    if (N_kk ~= ((n_sq_x+1)*(n_sq_y+1))),
                        no_grid = 1;
                    end;
                end;
            end;

            eval(['plot(ex_' num2str(kk) '(1,:),ex_' num2str(kk) '(2,:),''' colors(rem(kk-1,6)+1) '+'');']);

            hold on;
        end;
    end;
end;

hold off;
axis('equal');
if 1, %~no_grid,
    title('Reprojection error (in pixel) - To exit: right button');
else
    title('Reprojection error (in pixel)');
end;
xlabel('x');
ylabel('y');

set(5,'color',[1 1 1]);
set(5,'Name','error','NumberTitle','off');

if n_ima == 0,

        text(.5,.5,'No image data available','fontsize',24,'horizontalalignment' ,'center');

else

err_std = std(ex')';

fprintf(1,'Pixel error:          err = [ %3.5f   %3.5f] (all active images)\n\n',err_std);

b = 1;

while b==1,

    [xp,yp,b] = ginput3(1);

    if b==1,
        ddd = (ex(1,:)-xp).^2 + (ex(2,:)-yp).^2;

        [mind,indmin] = min(ddd);

        done = 0;
        kk_ima = 1;
        while (~done)&(kk_ima<=n_ima),
            %fprintf(1,'%d...',kk_ima);
            eval(['ex_kk = ex_' num2str(kk_ima) ';']);
            sol_kk = find((ex_kk(1,:) == ex(1,indmin))&(ex_kk(2,:) == ex(2,indmin)));
            if isempty(sol_kk),
                kk_ima = kk_ima + 1;
            else
                done = 1;
            end;
        end;

        eval(['x_kk = x_' num2str(kk_ima) ';']);
        xpt = x_kk(:,sol_kk);

        if ~no_grid,

            eval(['n_sq_x = n_sq_x_' num2str(kk_ima) ';']);
            eval(['n_sq_y = n_sq_y_' num2str(kk_ima) ';']);

            Nx = n_sq_x+1;
            Ny = n_sq_y+1;

            y1 = floor((sol_kk-1)./Nx);
            x1 = sol_kk - 1 - Nx*y1; %rem(sol_kk-1,Nx);

            y1 = (n_sq_y+1) - y1;
            x1 = x1 + 1;

            fprintf(1,'\n');
            fprintf(1,'Selected image: %d\n',kk_ima);
            fprintf(1,'Selected point index: %d\n',sol_kk);
            fprintf(1,'Pattern coordinates (in units of (dX,dY)): (X,Y)=(%d,%d)\n',[x1-1 y1-1]);
            fprintf(1,'Image coordinates (in pixel): (%3.2f,%3.2f)\n',[xpt']);
            fprintf(1,'Pixel error = (%3.5f,%3.5f)\n',[ex(1,indmin) ex(2,indmin)]);

        else

            fprintf(1,'\n');
            fprintf(1,'Selected image: %d\n',kk_ima);
            fprintf(1,'Selected point index: %d\n',sol_kk);
            fprintf(1,'Image coordinates (in pixel): (%3.2f,%3.2f)\n',[xpt']);
            fprintf(1,'Pixel error = (%3.5f,%3.5f)\n',[ex(1,indmin) ex(2,indmin)]);

        end;

        if exist(['wintx_' num2str(kk_ima)]),

            eval(['wintx = wintx_' num2str(kk_ima) ';']);
            eval(['winty = winty_' num2str(kk_ima) ';']);

            fprintf(1,'Window size: (wintx,winty) = (%d,%d)\n',[wintx winty]);
        end;

    end;

end;

disp('done');

end;

author: strawlab

repo: MultiCamSelfCal


text:
function V = vhs(im,varargin)
  % VHS Apply a VHS filter to an image
  %
  % V = vhs(im)
  % V = vhs(im,'ParameterName',ParameterValue,...)
  %
  % Inputs:
  %   im  h by w by c image
  %   Optional:
  %     'VerticalLoop' followed by whether to loop the bent strip vertically
  %     over time and output a sequence of images.
  % Output:
  %   V  h by w by c by f image of (f long sequence of images)
  %

  % must use double
  im = im2double(im);

  looping = false;
  % Map of parameter names to variable names
  params_to_variables = containers.Map( {'VerticalLoop'},{'looping'});
  v = 1;
  iter = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  if looping
    strip_top = 1;
  else
    strip_top = ceil(0.4*size(im,1));
  end

  first = 1;
  while true

    % http://mikeyjam.buzznet.com/user/journal/12237761/tutorial-getting-vhs-tv-effect/
    A = im;
    B = im;
    C = im;
    D = im;
    % exclusion blend like photoshop
    % http://www.deepskycolors.com/archive/2010/04/21/formulas-for-Photoshop-blending-modes.html
    ex = @(T,B) 0.5 - 2.*bsxfun(@times,T-0.5,B-0.5);
    % kill color channels
    A(:,:,1) = 0;
    B(:,:,2) = 0;
    C(:,:,3) = 0;
    % Shift color layers
    nudge = @(f) ceil(f*rand(1)*size(im,2));
    A = A(:,mod(nudge(0.02)+(1:end)-1,end)+1,:);
    B = B(:,mod(nudge(0.02)+(1:end)-1,end)+1,:);
    C = C(:,mod(nudge(0.02)+(1:end)-1,end)+1,:);
    A = A(mod(nudge(0.005)+(1:end)-1,end)+1,:,:);
    B = B(mod(nudge(0.005)+(1:end)-1,end)+1,:,:);
    C = C(mod(nudge(0.005)+(1:end)-1,end)+1,:,:);
    % exclusion blend colored layers and alpha blend with original

    F = D+0.3*(ex(ex(C,B),A)-D);
    N = rand(size(im));

    % inverse mapping function
    bend_w = (2*rand(1)-1)*5;
    bend = @(x,u) [mod(x(:,1)+(1-x(:,2)/max(x(:,2))).^2*bend_w,max(x(:,1))) x(:,2)];
    % maketform arguments
    ndims_in = 2;
    ndims_out = 2;
    tdata = [];
    tform = maketform('custom', 2,2, [], bend, tdata);

    % Bend strip
    strip_h = ceil(1/4*size(im,1));
    strip = mod(strip_top+(1:strip_h)-1,size(im,1))+1;
    F(strip,:,:) = imtransform(F(strip,:,:), tform);

    % overlay gray line
    ol= @(T,B) (T>0.5).*(1-(1-2.*(T-0.5)).*(1-B))+(T<=0.5).*((2.*T).*B);
    G = repmat(0.75,[numel(-1:1) size(F,2) size(F,3)]);
    F(mod(strip_top+(-1:1)-1,size(F,1))+1,:,:) = ...
      ol( F(mod(strip_top+(-1:1)-1,size(F,1))+1,:,:),G);

    % overlay horizontal lines
    L = zeros(size(F));
    L(1:4:end,:,:) = 1;
    L(2:4:end,:,:) = 1;
    L = imfilter(L,fspecial('gaussian',[5 5],1.5),'replicate');
    F = ol(F,L);

    % Fade in random color gradient
    m = rand(1,2)*2-1;
    R = bsxfun(@plus,m(1)*(1:size(F,2)),m(2)*(1:size(F,1))');
    R = (R-min(R(:)))./(max(R(:))-min(R(:)));
    R = gray2rgb(R,colormap(jet(255)));
    F = F+0.05*(ol(F,R)-F);

    % soft light
    sl = @(T,B) (B>0.5).*(1-(1-T).*(1-(B-0.5))) + (B<=0.5).*(T.*(B+0.5));
    F = F + 0.15*(sl(F,N)-F);

    % sharpen
    V(:,:,:,iter) = imsharpen(F);

    %imshow(V(:,:,:,iter));
    %drawnow;

    strip_top = strip_top+round(0.057*size(im,1));
    iter = iter+1;
    if ~looping || strip_top > size(im,1)
      break;
    end

  end
end

author: alecjacobson

repo: gptoolbox


text:
%% RATE OF CONVERGENCE OF BILINEAR FINITE ELEMENT METHOD
%
% This example is to show the rate of convergence of bilinear finite
% element approximation of the Poisson equation on the unit square with the
% following boundary conditions:
%
% - Non-empty Dirichlet boundary condition.
% - Pure Neumann boundary condition.
% - Robin boundary condition.
%
% The basis, data structure and numerical test is summarized in <a
% href="matlab:ifem PoissonQ1femrate">PoissonQ1femrate</a>.
%
% See also Poissonfemrate, Poissonafemrate
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

clear variable
%% Setting
[node,elem] = squarequadmesh([0,1,0,1],1/2^3);
mesh  = struct('node',node,'elem',elem);
option.L0 = 2;
option.maxIt = 4;
option.printlevel = 1;
option.plotflag = 0;
option.elemType = 'Q1';

%% Non-empty Dirichlet boundary condition.
pde = sincosdata;
mesh.bdFlag = setboundary(node,elem,'Dirichlet','~(x==0)','Neumann','x==0');
femPoisson(mesh,pde,option);

%% Pure Neumann boundary condition.
pde = sincosNeumanndata;
mesh.bdFlag = setboundary(node,elem,'Neumann');
femPoisson(mesh,pde,option);

%% Pure Robin boundary condition.
option.plotflag = 0;
pde = sincosRobindata;
mesh.bdFlag = setboundary(node,elem,'Robin');
femPoisson(mesh,pde,option);

%% Conclusion
% To do:
% 1. Fix the getH1errorQ1. The H1 error is computed wrong.
% 2. For pure Neuman boundary condition, the rate is not quite right.

repo: ifem


text:
function [rs, s] = fracdims(s, kmin, kmax, Nref, gstart, gend, past, steps)

%tstoolbox/@signal/fracdims
%   Syntax:
%     * rs = fracdims(s, kmin, kmax, Nref, gstart, gend, past, steps)
%     * rs = fracdims(s, kmin, kmax, Nref, gstart, gend, past)
%     * rs = fracdims(s, kmin, kmax, Nref, gstart, gend)
%
%   Input arguments:
%     * kmin - minimal number of neighbors for each reference point
%     * kmax - maximal number of neighbors for each reference point
%     * Nref - number of randomly chosen reference points (n == -1 means :
%       use all points)
%     * gstart - starting value for moments
%     * gend - end value for moments
%     * past - (optional) number of samples to exclude before and after
%       each reference index, default is 0
%     * steps - (optional) number of moments to calculate, default is 32
%
%   Compute fractal dimension spectrum D(q) using moments of neighbor
%   distances for time-delay reconstructed timeseries s.
%
%   Do the main job - computing nearest neighbors for reference points.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(6,8)

if nargin < 7
	past = 0
end
if nargin < 8
	steps = 32;
end

if steps < 1
	error('Number of steps must be positive')
end

if (gend-gstart > 0) & (steps > 1)
	gammas = linspace(gstart, gend, steps);
else
	gammas = gstart;
end

N = dlens(s,1);

try
	atria = optparams(s, 1);
	names = fieldnames(atria);
catch
	atria = nn_prepare(data(s), 'euclidian');
	s = setoptparams(s, 1, atria);
end

% Do the main job - computing nearest neighbors for reference points
[nn, dist] = nn_search(data(s), atria, randref(1,N,Nref), kmax, past);

out = gendimest(dist, gammas, kmin, kmin, kmax);

rs = signal(core(out(:)), s);

a = achse(unit, 1 - (gammas(:) ./ out));
a = setname(a, 'q');
rs = setaxis(rs, 1, a);
rs = addhistory(rs,  ['Computed fractal dimension spectrum']);
rs = addcommandlines(rs, 's = fracdims(s', kmin, kmax, Nref, gstart, gend, past, steps);
rs = setyname(rs, 'D(q)');
rs = setlabel(rs, 'Renyi dimension spectrum');

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [y,e] = powspec(x, sr, wintime, steptime, dither)
%[y,e] = powspec(x, sr, wintime, steptime, sumlin, dither)
%
% compute the powerspectrum and frame energy of the input signal.
% basically outputs a power spectrogram
%
% each column represents a power spectrum for a given frame
% each row represents a frequency
%
% default values:
% sr = 8000Hz
% wintime = 25ms (200 samps)
% steptime = 10ms (80 samps)
% which means use 256 point fft
% hamming window
%
% $Header: /Users/dpwe/matlab/rastamat/RCS/powspec.m,v 1.3 2012/09/03 14:02:01 dpwe Exp dpwe $

% for sr = 8000
%NFFT = 256;
%NOVERLAP = 120;
%SAMPRATE = 8000;
%WINDOW = hamming(200);

if nargin < 2
  sr = 8000;
end
if nargin < 3
  wintime = 0.025;
end
if nargin < 4
  steptime = 0.010;
end
if nargin < 5
  dither = 1;
end

winpts = round(wintime*sr);
steppts = round(steptime*sr);

NFFT = 2^(ceil(log(winpts)/log(2)));
%WINDOW = hamming(winpts);
%WINDOW = [0,hanning(winpts)'];
WINDOW = [hanning(winpts)'];
% hanning gives much less noisy sidelobes
NOVERLAP = winpts - steppts;
SAMPRATE = sr;

% Values coming out of rasta treat samples as integers,
% not range -1..1, hence scale up here to match (approx)
y = abs(specgram(x*32768,NFFT,SAMPRATE,WINDOW,NOVERLAP)).^2;

% imagine we had random dither that had a variance of 1 sample
% step and a white spectrum.  That's like (in expectation, anyway)
% adding a constant value to every bin (to avoid digital zero)
if (dither)
  y = y + winpts;
end
% ignoring the hamming window, total power would be = #pts
% I think this doesn't quite make sense, but it's what rasta/powspec.c does

% that's all she wrote

% 2012-09-03 Calculate log energy - after windowing, by parseval
e = log(sum(y));

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
function [] = getGtDepthMap(fnames,statesDir,annotationsDir,depthMapsDir,cadModels)
%GETGTDEPTHMAP Summary of this function goes here
%   Detailed explanation goes here

parfor i=1:length(fnames)
    warning off;
    disp(i);
    state=load(fullfile(statesDir,fnames{i}));
    state=state.state;
    dmap = zeros(size(state.mask));
    voc_id = fnames{i}(1:end-6);
    recordFile = fullfile(annotationsDir,[voc_id '.mat']);
    dmapFile = fullfile(depthMapsDir,fnames{i});
    if(~exist(dmapFile,'file'))
        if(exist(recordFile,'file'))
            record = load(recordFile);
            record = record.record;
            index = 0;
            for j=1:length(record.objects)
                bb = record.objects(j).bbox;
                bb2 = state.modelbbox;
                bb2(3:4)=bb2(3:4)+bb(1:2);
                if(sum(abs(bb-bb2))<=4)
                    index = j;
                end
            end
            if(index)
                vertex = cadModels(record.objects(index).cad_index).vertices;
                face = cadModels(record.objects(index).cad_index).faces;
                [x2d,D] = project_3d(vertex, record.objects(index));
                if(size(D,1)>0)
                    dmap = meshToDepth([round(x2d),-D],face,size(state.mask));
                end
            end
        end
        savefunc(dmapFile,dmap);
    end
end
end

function savefunc(file,dmap)
    save(file,'dmap');
end

repo: CategoryShapes


text:
function geneCopyFactor = childCopyGivenFreqsFactor(alleleFreqs, geneCopyVar)
% This function creates a factor whose values are the frequencies of each
% allele in the population.
%
% Input:
%   alleleFreqs: A list of the frequencies of the alleles in the population
%   genotypeVar: The variable number for the genotype
%
% Output:
%   geneCopyFactor: A factor for the prior probability of genotypeVar (note
%   that this is the FULL CPD with no evidence observed)

numAlleles = length(alleleFreqs);
geneCopyFactor = struct('var', [], 'card', [], 'val', []);
geneCopyFactor.var(1) = geneCopyVar;
geneCopyFactor.card(1) = numAlleles;
geneCopyFactor.val = alleleFreqs';

author: anhncs

repo: Probabilistic-Graphical-Models


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

%*****************************************************************************80
%
%% FRAC_BORROW_JULIAN borrows fractions from days in a Julian YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, D, real F,
%    a YMDF date.
%
  while ( f < 0.0 )

    f = f + 1.0;

    d = d - 1;

  end

  [ y, m, d ] = day_borrow_julian ( y, m, d );

  return
end


text:
function [refMesh, distListIters, refVerticesIters] = averageMesh(meshList, varargin)
%AVERAGEMESH Compute average mesh from a list of meshes.
%
%   AVG = averageMesh(MESHLIST)
%
%   Example
%   averageMesh
%
%   See also
%     meshes3d

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2020-01-31, using Matlab 9.7.0.1247435 (R2019b) Update 2
% Copyright 2020-2022 INRAE - BIA Research Unit - BIBS Platform (Nantes)

%% Parse input values

% default values
nIters = 10;
verbose = false;

% parse input arguments
while length(varargin) > 1
    name = varargin{1};
    if ~ischar(name)
        error('require parameter name-value pairs');
    end

    if strcmpi(name, 'verbose')
        verbose = varargin{2};
    elseif strcmpi(name, 'nIters')
        nIters = varargin{2};
    else
        error(['Unknown parameter name: ' name]);
    end
    varargin(1:2) = [];
end

%% Initialisations

nMeshes = length(meshList);

% initialize kd-trees to accelerate nearest-neighbor searches
treeList = cell(nMeshes, 1);
for iMesh = 1:nMeshes
    treeList{iMesh} = KDTreeSearcher(meshList{iMesh}.vertices);
end

% choose arbitrary initial mesh
refMesh = struct('vertices', meshList{1}.vertices, 'faces', meshList{1}.faces);

refVerticesIters = cell(1, nIters);
distListIters = cell(1, nIters);

%% Main iteration

for iIter = 1:nIters
    if verbose
        fprintf('iter %d/%d\n', iIter, nIters);
    end
    refMesh = smoothMesh(refMesh);

    % create new array for average vertices
    newVerts = zeros(size(refMesh.vertices));
    distList = zeros(size(refMesh.vertices, 1), 1);

    % iterate over all meshes
    for iMesh = 1:nMeshes
        if verbose
            fprintf('    mesh %d/%d\n', iMesh, nMeshes);
        end

%         mesh = meshList{iMesh};
        inds = knnsearch(treeList{iMesh}, refMesh.vertices);

%         closest = mesh.vertices(inds,:);
        closest = treeList{iMesh}.X(inds,:);
        newVerts = newVerts + closest;
        distList = distList + sum((closest - refMesh.vertices).^2, 2);
    end

    % update new vertices
    newVerts = newVerts / nMeshes;
    refVerticesIters{iIter} = newVerts;
    refMesh.vertices = newVerts;

    % keep list of distances
    distList = sqrt(distList / nMeshes);
    distListIters{iIter} = distList;
end

% figure; drawMesh(refMesh, 'lineStyle', 'none', 'faceColor', [.5 .5 .5])
% axis equal; view(3); hold on; axis([-2.5 2.5 -2 2 -3.5 3.5]); light;
% lighting gouraud
% title('Average mesh');
% print(gcf, 'averageMesh_initial.png', '-dpng');

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function test_suite = testIsTearDownString
%testIsTearDownString Unit tests for isTearDownString

%   Steven L. Eddins
%   Copyright 2008 The MathWorks, Inc.

initTestSuite;

function testOneStringIs
assertTrue(xunit.utils.isTearDownString('teardownfoobar'));
assertTrue(xunit.utils.isTearDownString('TearDown_foobar'));

function testOneStringIsNot
assertFalse(xunit.utils.isTearDownString('tEardown'));

function testCellArray
strs = {'teardown', 'tearup'};
assertEqual(xunit.utils.isTearDownString(strs), [true false]);
assertEqual(xunit.utils.isTearDownString(strs'), [true; false]);

author: ricedsp


text:
function [cl,cu] = matRad_getConstraintBounds(optiProb,cst)
% matRad IPOPT get constraint bounds wrapper function
%
% call
%   [cl,cu] = matRad_getConstraintBounds(optiProb,cst)
%
% input
%   cst:            matRad cst struct
%
% output
%   cl: lower bounds on constraints
%   cu: lower bounds on constraints
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2016 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

BPtype = class(optiProb.BP);
isEffectBP = strcmp(BPtype,'matRad_EffectProjection');

% Initialize bounds
cl = [];
cu = [];

% compute objective function for every VOI.
for  i = 1:size(cst,1)

    % Only take OAR or target VOI.
    if ~isempty(cst{i,4}) && ( isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') )

        % loop over the number of constraints for the current VOI
        for j = 1:numel(cst{i,6})

            optiFunc = cst{i,6}{j};

            % only perform computations for constraints
            %if ~isempty(strfind(cst{i,6}{j}.type,'constraint'))
            if isa(optiFunc,'DoseConstraints.matRad_DoseConstraint')

                if isEffectBP
                    doses = optiFunc.getDoseParameters();

                    effect = cst{i,5}.alphaX*doses + cst{i,5}.betaX*doses.^2;

                    optiFunc = optiFunc.setDoseParameters(effect);
                end

                 cl = [cl;optiFunc.lowerBounds(numel(cst{i,4}{1}))];
                 cu = [cu;optiFunc.upperBounds(numel(cst{i,4}{1}))];

                %end
            end

        end % over all objectives of structure

    end % if structure not empty and target or oar

end % over all structures

repo: matRad


text:
function [ll,lls,Ht] = gogarch_likelihood(parameters,data,p,q,gjrType,P,L,isOgarch,isInference)
% Log-likelihood for use in estimation GOGARCH and OGARCH models
%
% USAGE:
%  [LL,LLS,HT] = gogarch_likelihood(PARAMETERS,DATA,P,Q,GJRTYPE,P,L,ISOGARCH,ISINFERNCE)
%
% INPUTS:
%   PARAMETERS - K*(K-1)/2 + sum(P) + sum(Q) by 1 vector of parameters (only sum(P) + sum(Q) if ISOGARCH)
%   DATA       - K by K by Tarray of data (either realized measures or outer-products of daily data)
%   P          - K by 1 vector of positive, scalar integer representing the number of symmetric innovations
%   Q          - K by 1 vector of non-negative, scalar integer representing the number of conditional covariance lags
%   GJRTYPE    - K by 1 vector containing either 1 (TARCH/AVGARCH) or 2 (GJRGARCH/GARCH/ARCH)
%   P          - Eigenvector of unconditional covariance matrix of the data
%   L          - Diagonal matrix containing the eigenvalues of unconditional covariance matrix of the data
%   ISOGARCH   - Boolean indicating that the model is OGARCH (otherwise GOGARCH)
%   ISINFERNCE - Boolean indicating the likelihood is being used for inference, so that the first
%                  K(K+1)/2 parameters are ivech(S)
%
% OUTPUTS:
%   LL         - The log likelihood computed at PARAMETERS
%   LLS        - T by 1 vector of log-likelihoods
%   HT         - K by K by T vector of conditional covariances
%
% COMMENTS:

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 4/15/2012

[k,~,T] = size(data);
if size(parameters,1)>size(parameters,2)
    parameters = parameters';
end

offset = 0;
if isInference
    S =  ivech(parameters(1:k*(k+1)/2));
    [P,L] = eig(S);
    P = P';
    offset = offset + k*(k+1)/2;
end

if ~isOgarch
    phi = parameters(offset + (1:k*(k-1)/2));
    U = phi2u(phi);
    offset = offset + k*(k-1)/2;
else
    U = eye(k);
end
Z = P*L^(0.5)*U;
Zinv = U'*L^(-0.5)*P';
stdData = zeros(k,k,T);
for t=1:T
    stdData(:,:,t) = Zinv*data(:,:,t)*Zinv';
end
% Univariate GARCH models
V = zeros(T,k);
w = .06 * .94.^(0:sqrt(T));
w = w/sum(w);
likData = zeros(T,k);
for i=1:k
    count = p(i) + q(i);
    volParameters = parameters(offset + (1:count));
    volParameters = max(volParameters,0);
    volParameters = [1-sum(volParameters) volParameters]; %#ok<AGROW>
    offset = offset + count;
    likData(:,i) = squeeze(stdData(i,i,:));
    volData = sqrt(likData(:,i));
    backCast = w*volData(1:length(w)).^2;
    v = tarch_core_simple(volData,volParameters,backCast,0,p(i),0,q(i),gjrType(i));
    V(:,i) = v;
end

likConst = k*log(2*pi);
logdetZZp = log(det(Z*Z'));
lls = 0.5*(likConst + logdetZZp + sum(log(V),2) + sum(likData./V,2));
ll = sum(lls);

if ~isreal(ll) || isnan(ll) || isinf(ll)
    ll = 1e7;
end

if nargout>2
    Ht = zeros(k,k,T);
    for t=1:T
        Ht(:,:,t) = Z*diag(V(t,:))*Z';
    end
end

author: bashtage

repo: mfe-toolbox


text:
function xy_io_test ( )

%*****************************************************************************80
%
%% XY_IO_TEST tests the XY_IO library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 January 2009
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'XY_IO_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the XY_IO library.\n' );

  xy_io_test01 ( );
  xy_io_test02 ( );
  xy_io_test03 ( );
  xy_io_test04 ( );
  xy_io_test05 ( );
  xy_io_test06 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'XY_IO_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function A10 = synA10min(A11, A00, cmax)
%-----------------------------------------------------------------------------
%
% For each point of colour 10 this function assigns the minimum value at the
% neighbouring gridpoints of colours 11 and 00.
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: December 7, 2001.
% (c) 1998-2002 Stichting CWI, Amsterdam
%-----------------------------------------------------------------------------
[n00, m00]=size(A00);
[n11, m11]=size(A11);
n10=n11;
m10=m00;
%[n10, m10]=size(A10);
if     m10 == m11
  S=min(A11, stripR(extL(A11, cmax)));
elseif m10 == m11+1
  S=min(extL(A11, cmax), extR(A11, cmax));
else
  disp([' size A11 = ' int2str(size(A11)) ' size A00 = ' int2str(size(A00))]);
  error(' synA10min - A11 and A00 do not match ');
end
if     n10 == n00
  T=min(A00, stripU(extD(A00, cmax)));
elseif n10 == n00-1
  T=min(stripD(A00), stripU(A00));
else
  disp([' size A11 = ' int2str(size(A11)) ' size A00 = ' int2str(size(A00))]);
  error(' synA10min - A11 and A00 do not match ');
end
%Note: all(size(S) == size(T)) & all(size(S) == [n10 m10]) always holds.
A10=min(S, T);
%-----------------------------------------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function coeff = reg_als(basis_cr, y, r, n)
%REG_ALS Summary of this function goes here
%   d - dimension of the regression function
%   N - number of regression points
%   y - 1 x N array of values of the function to approximate
%   n - d x 1 array, stores number of basis functions in each dimension
%   r - d+1 x 1 array, stores TT-ranks of the tensor of coefficients of the
%       regression function
%   basis_cr - array that stores values of basis functions at regression
%   points
%   fixed - d x 1 cell of aggregated fixed dimensions of the regression
%   function at k-th step of the ALS algorithm
%   coeff - TT-tensor of coefficient of the regression function

N = size(y, 2);
d = numel(n);
coeff = tt_rand(n, d, r);
coeff_ps = coeff.ps;
coeff_cr = coeff.core;
basis_ps = cumsum([1 ; n*N]);

% assemble fixed part of the TT-function
fixed_ps = cumsum([1 ; N ; r(2:d)*N ; N]);
fixed_cr = zeros(fixed_ps(d+2) - fixed_ps(1), 1);
fixed_cr(fixed_ps(d+1): fixed_ps(d+2)-1) = ones(N, r(d+1));
for dim = d: -1: 2
    c = reshape(coeff_cr(coeff_ps(dim): coeff_ps(dim+1)-1), [r(dim) n(dim) r(dim+1)]);
    c = permute(c, [1 3 2]);
    b = reshape(basis_cr(basis_ps(dim): basis_ps(dim+1)-1), [n(dim) N]);
    t1 = ten_conv(c, 3, b); % r(dim) x r(dim+1) x N
    t2 = reshape(fixed_cr(fixed_ps(dim+1): fixed_ps(dim+2)-1), [r(dim+1), N]);
    t1 = reshape(t1, [r(dim), r(dim+1)*N]);
    t2 = kron(ones(r(dim),1), reshape(t2, [1, r(dim+1)*N]));
    t3 = t1.*t2;
    t3 = reshape(t3, [r(dim) r(dim+1) N]);
	t3 = sum(t3, 2);
    fixed_cr(fixed_ps(dim): fixed_ps(dim+1)-1) = reshape(t3, [r(dim)*N, 1]);
end
fixed_cr(fixed_ps(1): fixed_ps(2)-1) = ones(N, r(1));

% initial sum of squares
VAR = var(y);

n_swp = 0;
R2 = 0;
swp_tp = 'f';
dim = 1;
while R2 < 0.9999 && n_swp < 10
    n1 = r(dim);
    n2 = n(dim);
    n3 = r(dim+1);

    u1 = reshape(fixed_cr(fixed_ps(dim): fixed_ps(dim+1)-1), [N n1]);
    u2 = reshape(basis_cr(basis_ps(dim): basis_ps(dim+1)-1), [n2 N]);
    u3 = reshape(fixed_cr(fixed_ps(dim+1): fixed_ps(dim+2)-1), [n3 N]);

    t1 = kron(kron(ones(n3,1), ones(n2,1)), u1');
    t2 = kron(kron(ones(n3,1), u2), ones(n1,1));
    t3 = kron(kron(u3, ones(n2,1)), ones(n1,1));
    t = t1.*t2.*t3;

    A = t*t';
    b = t*y';
    c = A\b;
    coeff_cr(coeff_ps(dim):coeff_ps(dim+1)-1) = c;

    if (swp_tp == 'f' && dim < d) || (swp_tp == 'b' && dim == 1)
        c = permute(reshape(c, [n1 n2 n3]), [1 3 2]);
        t1 = ten_conv(c, 3, u2);
        t1 = reshape(t1, [n1*n3, N]);
        u1 = kron(ones(1,n3), u1);
        t2 = u1.*t1';
        t2 = reshape(t2, [N n1 n3]);
        t2 = sum(t2, 2);
        t2 = reshape(t2, [N, n3]);
        fixed_cr(fixed_ps(dim+1): fixed_ps(dim+2)-1) = reshape(t2, [N*n3, 1]);

        %R^2
        if swp_tp == 'b'
            MSE = sum((sum(t2'.*u3,1) - y).^2)/N;
            R2 = 1 - MSE/VAR;
        end
    elseif (swp_tp == 'b' && dim > 1) || (swp_tp == 'f' && dim == d)
        c = permute(reshape(c, [n1 n2 n3]), [1 3 2]);
        t1 = ten_conv(c, 3, u2);
        t1 = reshape(t1, [n1, n3*N]);
        u3 = kron(ones(n1,1), reshape(u3, [1, n3*N]));
        t2 = t1.*u3;
        t2 = reshape(t2, [n1 n3 N]);
        t2 = sum(t2, 2);
        t2 = reshape(t2, [n1 N]);
        fixed_cr(fixed_ps(dim): fixed_ps(dim+1)-1) = reshape(t2, [n1*N, 1]);

        %R^2
        if swp_tp == 'f'
            MSE = sum((sum(u1'.*t2,1) - y).^2)/N;
            R2 = 1 - MSE/VAR;
        end
    end

    if swp_tp == 'f'
        if dim < d
            dim = dim + 1;
        elseif dim == d
            swp_tp = 'b';
            dim = dim - 1;
            n_swp = n_swp + 1;
            fprintf('=tt_reg_als= Sweep %d, R^2: %.2f%%, MSE: %1.2e\n', n_swp, 100*R2, MSE);
        end
    elseif swp_tp == 'b'
        if dim > 1
            dim = dim - 1;
        elseif dim == 1
            swp_tp = 'f';
            dim = dim + 1;
            n_swp = n_swp + 1;
            fprintf('=tt_reg_als= Sweep %d, R^2: %.2f%%, MSE: %1.2e\n', n_swp, 100*R2, MSE);
        end
    end
end

coeff.core = coeff_cr;

end

author: oseledets

repo: TT-Toolbox


text:
function [Pgood,Pbad] = NDS(Population,K)
% Sort the population based on non-dominated sorting and crowding distance

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    FrontNo  = NDSort(Population.objs,inf);
    CrowdDis = CrowdingDistance(Population.objs,FrontNo);
    [~,rank] = sortrows([FrontNo;-CrowdDis]');
    Pgood    = Population(rank(1:K));
    Pbad     = Population(rank(end-K+1:end));
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [pt_errs_phys, pts_moved_pix, TRE_phys, TREstd_phys] = DIR_movepoints_v2(pts_mov_pix, ...
                        pts_fix_pix, spc_orig, spc_disp, D)
    pts_moved_pix = pts_fix_pix;
    D_pix = conv_3d_T_from_phys_to_pix(D, spc_disp);
    for i = 1 : size(pts_mov_pix, 1)
        x = fl( pts_fix_pix(i, :) - 1).*spc_orig(:)./spc_disp + 1;
        d1 = lerp_eval(x, D_pix(:,:,:,1));
        d2 = lerp_eval(x, D_pix(:,:,:,2));
        d3 = lerp_eval(x, D_pix(:,:,:,3));

        pts_moved_pix(i, :) =  pts_moved_pix(i,:) + [d1, d2, d3] .*spc_disp./spc_orig;
    end
    pts_moved_pix = round(pts_moved_pix);
    koef = repmat(spc_orig, [size(pts_moved_pix, 1), 1]);
    pt_errs_phys = sqrt( sum((  (pts_moved_pix - pts_mov_pix).*koef  ).^2, 2) );
    TRE_phys = mean(pt_errs_phys);
    TREstd_phys = std(sqrt( sum((  (pts_moved_pix - pts_mov_pix).*koef  ).^2, 2) ));
end

repo: pTVreg


text:
function zz=lpcss2zz(ss)
%LPCSS2ZZ Convert s-place poles to z-plane poles ZZ=(SS)
%the s-plane is in units of Normalized Hz and so the imaginary part
% of each ss() value is in the range +-0.5
%
% If you multiply ss by the sample frequency, a formant with
% frequency f and bandwidth b will give an s-plane pole-pair
% of approximately -b/2 +- j f

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: lpcss2zz.m,v 1.4 2007/05/04 07:01:39 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

zz=exp(2*pi*ss);

author: decouples


text:
function innerpro = cinnerprodgeneral(x, y)
% Computes the Euclidean inner product between x and y in the complex case
%
% function innerpro = cinnerprodgeneral(x, y)
%
% The input x and y are numeric data structures which can be defined
% recursively by arrays, structs and cells. Each part of x and y should
% be a struct which contains the fields real and iamg which indicate
% the real and imaginary part of the stored complex numbers. The inner
% product between x and y is defined as sum(real(conj(x(:)).*y(:))).
% The return is the sum of the inner products over each part of x and y.
% In case that x and y are structs with different fields, the inner products
% are computed only for the common fields.
%
% Note: Operations between dlarrays containing complex numbers have been
% introduced in Matlab R2021b or later. This file is only useful for Matlab
% R2021a or earlier. It will be discarded when Matlab R2021b is stable.
%
% See also: innerprodgeneral, manoptADhelp

% This file is part of Manopt: www.manopt.org.
% Original author: Xiaowen Jiang, Aug. 31, 2021.
% Contributors: Nicolas Boumal
% Change log:

    if ~((isstruct(x) && isstruct(y)) || (iscell(x) && iscell(y))...,
            || (isnumeric(x) && isnumeric(y)) || (isstruct(x) && isnumeric(y)))

        up = MException('manopt:autodiff:cinnerprodgeneral' ,...
            'cinnerprodgeneral should only accept structs, cells or arrays.');
        throw(up);

    end
    % recursively compute the inner product
    if isstruct(x) && isstruct(y) && (~isfield(x,'real')) && (~isfield(y,'real'))
        innerpro  = cinnerprodgeneral_struct(x,y);
    elseif iscell(x) && iscell(y)
        innerpro = cinnerprodgeneral_cell(x,y);
    else
        xconj = cconj(x);
        product = cdottimes(xconj,y);
        innerpro = sum(creal(product),'all');
        % slower
        % xcol = cmat2col(x);
        % innerpro = creal(cprod(ctransp(xcol),xcol));
    end

    % struct case
    function innerpro = cinnerprodgeneral_struct(x,y)
        innerpro = 0;
        elemsx = fieldnames(x);
        elemsy = fieldnames(y);
        % find the common fields
        [elems,ix,iy] = intersect(elemsx,elemsy, 'stable');
        nelems = numel(elems);
        for ii = 1:nelems
            if isstruct(x.(elemsx{ix(ii)})) && (~isfield(x.(elemsx{ix(ii)}),'real'))...,
                    && (~isfield(y.(elemsy{iy(ii)}),'real'))
                innerpro = innerpro + cinnerprodgeneral_struct(...,
                    x.(elemsx{ix(ii)}),y.(elemsy{iy(ii)}));
            elseif iscell(x.(elemsx{ix(ii)}))
                innerpro = innerpro + cinnerprodgeneral_cell(...,
                    x.(elemsx{ix(ii)}),y.(elemsy{iy(ii)}));
            else
                xconj = cconj(x.(elemsx{ix(ii)}));
                product = cdottimes(xconj, y.(elemsy{iy(ii)}));
                innerpro = innerpro + sum(creal(product), 'all');
            end
        end
    end

    % cell case
    function innerpro = cinnerprodgeneral_cell(x,y)
        innerpro = 0;
        ncell = length(x);
        for ii = 1:ncell
            if isstruct(x{ii}) && (~isfield(x{ii},'real')) && (~isfield(y{ii},'real'))
                innerpro = innerpro + cinnerprodgeneral_struct(...,
                    x{ii},y{ii});
            elseif iscell(x{ii})
                innerpro = innerpro + cinnerprodgeneral_cell(...,
                    x{ii},y{ii});
            else
                xconj = cconj(x{ii});
                product = cdottimes(xconj, y{ii});
                innerpro = innerpro + sum(creal(product), 'all');
            end
        end
    end

end

author: NicolasBoumal

repo: manopt


text:
function pass = test_ultrapts(pref)

% Choose a tolerance:
tol = 1e-14;

% Test a small n (using REC)
n = 42;
lambda = .3;
x = ultrapts(n, lambda);
pass(1) = all(size(x) == [n, 1]);
[x, w, v] = ultrapts(n, lambda);
pass(2) = all(size(x) == [n, 1]) && all(size(w) == [1, n]) && all(size(v) == [n, 1]);
pass(3) = abs(w*x) < tol && abs(w*x.^2 - 0.8843414686338345) < tol;
pass(4) = abs(x(37) - 9.131896381993957E-01) < tol;
pass(5) = abs(w(37) - 4.332670514309510E-02) < tol;
pass(6) = abs(v(37) - 3.115587460502451E-01) < tol;

% Test on [0, 10]:
[x, w, v] = ultrapts(n, lambda, [0, 10]);
pass(7) = abs(w*x - 30.1957169274024) < 31*tol && abs(w*x.^2 - 209.0472710358626) < 300*tol;
pass(8) = abs(x(38) - 9.704505777068543E+00) < tol;
pass(9) = abs(w(38) - 1.018229378664342E-01) < tol;
pass(10) = abs(v(38) + 2.449177929215358E-01) < tol;

% Test a larger n
lambda = .8;
n = 251;
x = ultrapts(n, lambda);
pass(11) = all(size(x) == [n, 1]);
[x, w, v] = ultrapts(n, lambda);
pass(12) = all(size(x) == [n, 1]) && all(size(w) == [1, n]) && all(size(v) == [n, 1]);
pass(13) = abs(w*x) < tol && abs(w*x.^2 - 0.4744211549960596) < tol;
pass(14) = abs(x(37) + 8.958806879214126E-01) < tol;
pass(15) = abs(w(37) - 3.406945649865882E-03) < tol;
pass(16) = abs(v(37) - 2.321704534650446E-01) < tol;

% Test on [0, 10]:
[x, w, v] = ultrapts(n, lambda, [0, 10]);
pass(17) = abs(w*x - 112.1472319135050) < 100*tol && abs(w*x.^2 - 716.4962038918374) < 100*tol;
pass(18) = abs(x(38) - 5.486606034997460E-01 ) < tol;
pass(19) = abs(w(38) - 4.655102134393607E-02) < tol;
pass(20) = abs(v(38) + 2.427562206703888E-01 ) < tol;

% Test a larger n with a larger lambda
lambda = 7;
n = 551;
x = ultrapts(n, lambda,'asy');
pass(21) = all(size(x) == [n, 1]);
[x, w, v] = ultrapts(n, lambda,'asy');
pass(22) = all(size(x) == [n, 1]) && all(size(w) == [1, n]) && all(size(v) == [n, 1]);
pass(23) = abs(w*x) < tol && abs(w*x.^2 - (429/32768)*pi) < tol;
pass(24) = abs(x(37) + 9.748144265829347E-01) < tol;
pass(25) = abs(w(37) - 4.244751593204416E-12) < tol;
pass(26) = abs(v(37) - 6.123401324799126E-06) < tol;

% Test a larger n with a larger lambda
lambda = 25.5;
n = 5e3;
x = ultrapts(n, lambda,'asy');
pass(27) = all(size(x) == [n, 1]);
[x, w, v] = ultrapts(n, lambda,'asy');
pass(28) = all(size(x) == [n, 1]) && all(size(w) == [1, n]) && all(size(v) == [n, 1]);
pass(29) = abs(w*x) < tol && abs(w*x.^2 - 281474976710656/42710983650155457) < tol;
pass(30) = abs(sum(w) - 281474976710656/805867616040669 ) < tol;

% Test n = 1:
lambda = .6;
[x, w] = ultrapts(1, lambda);
pass(31) = (x == 0);
pass(32) = abs(w - 1.887181162535959) < tol;
[x, w] = ultrapts(1, lambda, [-10, 3]);
pass(33) = abs( (x + 10) + (x-3) ) < tol ; % midpoint
pass(34) = abs( w*x + 62.42774750787926 ) < 2*tol;

% Test n = 2:
[x, w] = ultrapts(2, lambda);
pass(35) = all(abs(x - [-sqrt(5)/4; sqrt(5)/4]) < tol);
pass(36) = all(abs(w - [.5*gamma(lambda+.5)*sqrt(pi)/gamma(lambda+1),...
    .5*gamma(lambda+.5)*sqrt(pi)/gamma(lambda+1)]) < tol);
[x, w] = ultrapts(2, lambda, [-10, 3]);
pass(37) = abs(sum(w) - 17.83649928796550) < tol;
pass(38) = abs( w*x + 62.42774750787926 ) < 10*tol;
pass(39) = abs( w*x.^2 - 453.9946459389969 ) < 100*tol;
pass(40) = abs( w*x.^3 + 3237.463968416426 ) < 100*tol;
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function out = runtests(varargin)
%runtests Run unit tests
%   runtests runs all the test cases that can be found in the current directory
%   and summarizes the results in the Command Window.
%
%   Test cases can be found in the following places in the current directory:
%
%       * An M-file function whose name starts or ends with "test" or
%         "Test" and that returns no output arguments.
%
%       * An M-file function whose name starts or ends with "test" or
%         "Test" and that contains subfunction tests and uses the
%         initTestSuite script to return a TestSuite object.
%
%       * An M-file defining a subclass of TestCase.
%
%   runtests(dirname) runs all the test cases found in the specified directory.
%
%   runtests(packagename) runs all the test cases found in the specified
%   package. (This option requires R2009a or later).
%
%   runtests(mfilename) runs test cases found in the specified function or class
%   name. The function or class needs to be in the current directory or on the
%   MATLAB path.
%
%   runtests('mfilename:testname') runs the specific test case named 'testname'
%   found in the function or class 'name'.
%
%   Multiple directories or file names can be specified by passing multiple
%   names to runtests, as in runtests(name1, name2, ...) or
%   runtests({name1, name2, ...}, ...)
%
%   runtests(..., '-verbose') displays the name and result, result, and time
%   taken for each test case to the Command Window.
%
%   runtests(..., '-logfile', filename) directs the output of runtests to
%   the specified log file instead of to the Command Window.
%
%   out = runtests(...) returns a logical value that is true if all the
%   tests passed.
%
%   Examples
%   --------
%   Find and run all the test cases in the current directory.
%
%       runtests
%
%   Find and run all the test cases in the current directory. Display more
%   detailed information to the Command Window as the test cases are run.
%
%       runtests -verbose
%
%   Save verbose runtests output to a log file.
%
%       runtests -verbose -logfile my_test_log.txt
%
%   Find and run all the test cases contained in the M-file myfunc.
%
%       runtests myfunc
%
%   Find and run all the test cases contained in the TestCase subclass
%   MyTestCase.
%
%       runtests MyTestCase
%
%   Run the test case named 'testFeature' contained in the M-file myfunc.
%
%       runtests myfunc:testFeature
%
%   Run all the tests in a specific directory.
%
%       runtests c:\Work\MyProject\tests
%
%   Run all the tests in two directories.
%
%       runtests c:\Work\MyProject\tests c:\Work\Book\tests

%   Steven L. Eddins
%   Copyright 2009-2010 The MathWorks, Inc.

verbose = false;
logfile = '';
if nargin < 1
    suite = TestSuite.fromPwd();
else
    [name_list, verbose, logfile] = getInputNames(varargin{:});
    if numel(name_list) == 0
        suite = TestSuite.fromPwd();
    elseif numel(name_list) == 1
        suite = TestSuite.fromName(name_list{1});
    else
        suite = TestSuite();
        for k = 1:numel(name_list)
            suite.add(TestSuite.fromName(name_list{k}));
        end
    end
end

if isempty(suite.TestComponents)
    error('xunit:runtests:noTestCasesFound', 'No test cases found.');
end

if isempty(logfile)
    logfile_handle = 1; % File handle corresponding to Command Window
else
    logfile_handle = fopen(logfile, 'w');
    if logfile_handle < 0
        error('xunit:runtests:FileOpenFailed', ...
            'Could not open "%s" for writing.', logfile);
    else
        cleanup = onCleanup(@() fclose(logfile_handle));
    end
end

fprintf(logfile_handle, 'Test suite: %s\n', suite.Name);
if ~strcmp(suite.Name, suite.Location)
    fprintf(logfile_handle, 'Test suite location: %s\n', suite.Location);
end
fprintf(logfile_handle, '%s\n\n', datestr(now));

if verbose
    monitor = VerboseTestRunDisplay(logfile_handle);
else
    monitor = TestRunDisplay(logfile_handle);
end
did_pass = suite.run(monitor);

if nargout > 0
    out = did_pass;
end

function [name_list, verbose, logfile] = getInputNames(varargin)
name_list = {};
verbose = false;
logfile = '';
k = 1;
while k <= numel(varargin)
    arg = varargin{k};
    if iscell(arg)
        name_list = [name_list; arg];
    elseif ~isempty(arg) && (arg(1) == '-')
        if strcmp(arg, '-verbose')
            verbose = true;
        elseif strcmp(arg, '-logfile')
            if k == numel(varargin)
                error('xunit:runtests:MissingLogfile', ...
                    'The option -logfile must be followed by a filename.');
            else
                logfile = varargin{k+1};
                k = k + 1;
            end
        else
            warning('runtests:unrecognizedOption', 'Unrecognized option: %s', arg);
        end
    else
        name_list{end+1} = arg;
    end
    k = k + 1;
end


text:
classdef Constant_LineSearchInitiator < LineSearchInitiator

    properties (Access = private)
        value
    end

    methods (Access = public)

        function obj = Constant_LineSearchInitiator(cParams)
            obj.init(cParams);
            obj.value = cParams.value;
        end

        function initStep = compute(obj,lastStep)
            initStep = obj.value;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [cols] = cubicle2col(vol,scale)
% convert a 3D volume to columns, i.e, each column is a small cube in the
% 3D volume. (An extension of im2col).
% vol: 3D volume.
% scale: the size of the small cube. The size of the volume should be
% divisible by scale.

if ndims(vol) < 3
    error('input data should be a volume');
end

a_vol = vol;
square_size = size(a_vol);
block_size = square_size(1:3) ./ scale;
assert(all(block_size == floor(block_size)));

cols = zeros([scale^2, numel(vol)./scale^3], class(vol));
for c = 1 : scale ^ 3
    x_off = mod(c-1 , scale); y_off = mod(floor((c-1) / scale), scale); z_off = floor((c-1) / scale^2 );
    per_channel = vol((x_off + 1) : scale : end, (y_off + 1) : scale : end , (z_off + 1) : scale : end, :);
    cols(c, :) = per_channel(:);
end

author: zhirongw


text:
function Z = sparsemultipart(Z1,Z2,info)

% Find the elements in Z1 that are in the convex hull of Z2, where Z2 is bipartite

% This file is part of SOSTOOLS - Sum of Squares Toolbox ver 3.00.
%
% Copyright (C)2002, 2004, 2013  A. Papachristodoulou (1), J. Anderson (1),
%                                G. Valmorbida (1), S. Prajna (2),
%                                P. Seiler (3), P. A. Parrilo (4)
% (1) Department of Engineering Science, University of Oxford, Oxford, U.K.
% (2) Control and Dynamical Systems - California Institute of Technology,
%     Pasadena, CA 91125, USA.
% (3) Aerospace and Engineering Mechanics Department, University of
%     Minnesota, Minneapolis, MN 55455-0153, USA.
% (4) Laboratory for Information and Decision Systems, M.I.T.,
%     Massachusetts, MA 02139-4307
%
% Send bug reports and feedback to: sostools@cds.caltech.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.

% AP Feb 01 0

[Y,I] = find(sum(Z2,1)==0);
sizeinf = length(info);

if sizeinf == 1
   error('Error in sparsemultipart option - at least two sets are required');
end

for i = 1:sizeinf
  Z3{i} = inconvhull(Z1(:,info{i}),Z2(:,info{i}));
end

for i = 1:sizeinf-1%GVcomment for the variables introducing homogeneous monomials
  Z3{i+1} = [repmat(Z3{i},size(Z3{i+1},1),1),kron(Z3{i+1},ones(size(Z3{i},1),1))];
end
Z3 = Z3{sizeinf};
Z = zeros(size(Z3,1),size(Z2,2));
lgth = 0;
for i = 1:sizeinf
   Z(:,info{i})= Z3(:,(1+lgth):(lgth+length(info{i})));
   lgth = length(info{i})+lgth;
end
Z(:,I) = zeros(size(Z,1),size(I,2));

author: yu-jiang

repo: radpbook

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


text:
function out = cfg_example_run_sum(job)
% Example function that returns the sum of an vector given in job.a in out.
% The output is referenced as out(1), this is defined in
% cfg_example_vout_sum.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_example_run_sum.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

out = sum(job.a);

author: spm


text:
function [mi entropy fd_bins] = mutualinformationx(x,y,fd_bins,permtest)
% MUTUALINFORMATIONX   Compute mutual information between two vectors
%
%   Inputs:
%       x,y     :  data matrices of equal size
%
%   Optional inputs:
%       bins    :  number of bins to use for distribution discretization
%      permtest :  perform permutation test and return mi in standard-Z values
%
%   Outputs:
%       mi      :  mutual information in bits
%       entropy :  entropy of x, y, and joint
%       nbins   :  number of bins used for discretization
%                  (based on Freedman-Diaconis rule)
%
%  Mike X Cohen (mikexcohen@gmail.com)

if nargin<2, error('Specify two inputs.'); end
if length(x)~=length(y), error('X and Y must have equal length'); end

%% determine the optimal number of bins for each variable

% vectorize in the case of matrices
x=x(:); y=y(:);

if nargin<3 || isempty(fd_bins)
    n            = length(x);
    maxmin_range = max(x)-min(x);
    fd_bins1     = ceil(maxmin_range/(2.0*iqr(x)*n^(-1/3))); % Freedman-Diaconis

    n            = length(y);
    maxmin_range = max(y)-min(y);
    fd_bins2     = ceil(maxmin_range/(2.0*iqr(y)*n^(-1/3)));

    % and use the average...
    fd_bins = ceil((fd_bins1+fd_bins2)/2);
end

%% bin data

edges = linspace(min(x),max(x),fd_bins+1);
[nPerBin1,bins1] = histc(x,edges);

edges = linspace(min(y),max(y),fd_bins+1);
[nPerBin2,bins2] = histc(y,edges);

%% compute entropies

% recompute entropy with optimal bins for comparison
hdat1 = hist(x,fd_bins);
hdat1 = hdat1./sum(hdat1);
hdat2 = hist(y,fd_bins);
hdat2 = hdat2./sum(hdat2);

% convert histograms to probability values
for i=1:2
    eval([ 'entropy(' num2str(i) ') = -sum(hdat' num2str(i) '.*log2(hdat' num2str(i) '+eps));' ]);
end

%% compute joint probabilities

jointprobs = zeros(fd_bins);
for i1=1:fd_bins
    for i2=1:fd_bins
        jointprobs(i1,i2) = sum(bins1==i1 & bins2==i2);
    end
end
jointprobs=jointprobs./sum(jointprobs(:));

entropy(3) = -sum(jointprobs(:).*log2(jointprobs(:)+eps));

%% mutual information

mi = sum(entropy(1:2)) - entropy(3);

%% optional permutation testing

if nargin==4

    npermutes = 500;
    n = length(bins2);

    perm_mi = zeros(1,npermutes);

    for permi=1:npermutes

        jointprobs = zeros(fd_bins);

        % shuffle bins
        binbreak = randsample(round(n*.8),1,1)+round(n*.1);
        switch mod(permi,4)
            case 0, bins2 = [ bins2(binbreak:end);    bins2(1:binbreak-1) ];
            case 1, bins2 = [ bins2(end:-1:binbreak); bins2(1:binbreak-1) ];
            case 2, bins2 = [ bins2(binbreak:end);    bins2(binbreak-1:-1:1) ];
            case 3, bins2 = [ bins2(end:-1:binbreak); bins2(binbreak-1:-1:1) ];
        end

        for i1=1:fd_bins
            for i2=1:fd_bins
                jointprobs(i1,i2) = sum(bins1==i1 & bins2==i2);
            end
        end
        jointprobs=jointprobs./sum(jointprobs(:));

        perm_jentropy = -sum(jointprobs(:).*log2(jointprobs(:)+eps));

        % mutual information
        perm_mi(permi) = sum(entropy(1:2)) - perm_jentropy;
    end

    mi = (mi-mean(perm_mi))/std(perm_mi);
end

%%
% simplified replacement for randsample
function y = randsample(x,n,junk)
y=randperm(x);
y=y(1:n);

author: mikexcohen

repo: AnalyzingNeuralTimeSeries


text:
% Descriptor for simple shapes (e.g letters).
% Computes the descriptor twice in order to describe horizontally and vertically with the same accuracy.
%
% Input: image (2-D, binary or gray)
%        depth of recursion (default 7),
%        plotFlag(if on - plot Illustration)
% Output:
%        vec - descriptor: Division locations, values around zero. Of length 2*(2^depth - 2)
%        (levels: the depth of the division locations.)
%
% Shahar Armon. 3/1/2012

function  [vec levels] = hierarchicalCentroid(im, depth, plotFlag)
    if nargin < 2
        depth = 7;
        plotFlag = 0;
    elseif nargin < 3
        plotFlag = 0;
    end

    im(im < 0) = 0;

    [vec1 levels1] = hierarchicalCentroid1(im  ,depth, plotFlag);
    [vec2 levels2] = hierarchicalCentroid1(im' ,depth, 0);  % transpose iamge

    vec = [vec1, vec2];
    levels = [levels1,-levels2];  % Minus is only to mark the transpose.
end

% The descriptor. Includes Normalization.
function  [vec levels] = hierarchicalCentroid1(im, d, plotFlag)

    p = hierarchicalCentroidRec(im, 1, d);

    if plotFlag  % Illustration
        showLines(im, p);
    end

    meany = ([1:size(im,1)]*sum(im,2))/sum(sum(im));
    indVer = (mod(p(:,2),2) == 1);
    % Normalization for location:
    p(indVer,1) = p(indVer,1) - p(1,1);
    p(~indVer,1) = p(~indVer,1) - meany;

    % Normalization for size (keeping aspect ratio):
    p(:,1) = p(:,1)/size(im,1);
    % Normalization for size:
    % p(indVer,1) = p(indVer,1)/size(im,2);
    % p(~indVer,1) = p(~indVer,1)/size(im,1);

    % Illustration of the lines after normalization without the image
    %     showLines([], p);

    vec = p(2:end,1)';
    levels = p(2:end,2)';
end

function p = hierarchicalCentroidRec(im, depth, maxDepth)
    if depth > maxDepth
        p = [];
    else
        area = sum(sum(im));
        [rows,cols] = size(im);
        % compute the centroid-x
        if cols == 1
            centroid = 0.5;
        elseif area == 0
            centroid = cols/2;
        elseif rows == 1
            centroid = (im*[1:cols]')/area;
        else
            centroid = sum(im)*[1:cols]'/area;
        end

        leftIm = im(:,1:floor(centroid));
        rightIm = im(:,ceil(centroid):end);

        pLeft  = hierarchicalCentroidRec(leftIm'  , depth+1, maxDepth);
        pRight = hierarchicalCentroidRec(rightIm' , depth+1, maxDepth);

        %  Updates the distances so they will be in relation to the complete image
        if size(pRight,1) > 1
            ind = find((1-mod(depth,2))-mod(pRight(2:end,2),2)) + 1;
            pRight(ind,1) = pRight(ind,1) + ceil(centroid) - 1;
            pRight(1,1) = pRight(1,1) - 1;
        end

        p = [centroid, depth; pLeft; pRight];
    end
end

function showLines(im, p)
    if max(p(:,1))<3
        hold on;
        axis([-1 1 -1 1])
        axis ij;
        temp = mod(p(:,2),2);
        indVer = find(temp);
        indHor = find(~temp);
        showLinesRec(min(p(indVer,1)),max(p(indVer,1)),min(p(indHor,1)),max(p(indHor,1)), p);
    else
        imshow(im);
        hold on;
        showLinesRec(0,size(im,2),0,size(im,1), p);
    end
    hold off
end

function p = showLinesRec(x1,x2,y1,y2, p)
    p1 = p(1,:);
    w = 1+max(p(:,2))-p1(1,2);
    t = p1(1);
    p = p(2:end,:);
    if mod(p1(1,2),2)
        plot([t t],[y1 y2],'-m','LineWidth',w);
    else
        plot([x1 x2],[t t],'-r','LineWidth',w);
    end
    if (p1(1,2) < max(p(:,2)))
        if mod(p1(1,2),2)
            p = showLinesRec(x1,t,y1,y2, p);
            p = showLinesRec(t,x2,y1,y2, p);
        else
            p = showLinesRec(x1,x2,y1,t, p);
            p = showLinesRec(x1,x2,t,y2, p);
        end
    end

end


text:
% Online PRMF (Wang et al. 2012)
% process_video('RPCA', 'OPRMF', 'dataset/demo.avi', 'output/demo_OPRMF.avi');
X = normalize(M);
rk = 2;
lambdaU = 1;
lambdaV = 1;
tol = 1e-2;
mask = ones(size(X));
[~, ~, L] = onlineRPMF(X, rk, lambdaU, lambdaV, tol, mask);
S = X - L;

author: andrewssobral

repo: lrslibrary


text:
function [dat] =  test(algo,dat,lossType)
%
%               [res]=test(algo,data,loss)
%   An algorihm algo is trained on the training set data. If a loss
%   function is supplied, it will be calculated.
%   It returns a data object containing the results.
%   If no loss function is supplied, the data object will contain the loss
%   on the training set. If not, it usually contain the label estimates in
%   the X part and the true labels in the Y part of the data object.
%
%   It is also possible to call test(algo), where the emtpy data set is passed
%   into the algorithm to test it. This is useful for methods that generate their
%   own data.
%
%   Note:
%               test(algo) <=> test(algo,[]).

%% Programming note:
% It is used to call testing.m function which doesn't include loss
% function calculations in a child object

if nargin==1
    dat=[];% <--- data is optional
end;

if nargin<3
    lossType=[];
end;

if iscell(dat)
    dat=group(dat);
end;
%%<<----test data---->>
if ~isempty(dat)
    if ~am_i_data(dat)
        dat=test(dat); %% <--- i.e generate data
    end
end

%%%% SGE support (note sightly different behaviour from algorithm/train)
if isdeferred(algo)
	if submitted(algo.deferred)
		[algo jobfailed] = waitcollect(algo);
		if jobfailed, return, end
	else
		dat.deferred = qsub(algo.deferred, max(nargout,1), mfilename, algo, dat, lossType);
		return
	end
end
if ~isempty(dat)
if isdeferred(dat)
	[dat jobfailed] = waitcollect(dat);
	if jobfailed, return, end
end
end
%%%%

e=struct(dat);
%%<<------ if there are multiple datasets as input ---->>
if isa(dat,'group')  & strcmp(e.group,'separate')
  dat=e.child; % dat seen as @algorithm object => so it is not possible to access child directly
  dat=make_cell(dat);
  res=[];
  for i=1:length(dat)
      [r]=test(algo,dat{i},lossType);
      res{i}=r;
  end
  dat=group(res);  %% <--- return data/res as set of group
else
    [dat]=testing(algo,dat);
    if ~isempty(lossType)
      if isa(lossType,'loss')
	dat=train(lossType,dat);
      else
	dat=loss(dat,lossType);
      end
    end;
    if iscell(dat)
        dat=group(dat);
    end;  %%  <--- return data/res as set of group
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function W = SimGraph_NearestNeighbors(M, k, Type, sigma)
% SIMGRAPH_NEARESTNEIGHBORS Returns kNN similarity graph
%   Returns adjacency matrix for an k-Nearest Neighbors
%   similarity graph
%
%   'M' - A d-by-n matrix containing n d-dimensional data points
%   'k' - Number of neighbors
%   'Type' - Type if kNN Graph
%      1 - Normal
%      2 - Mutual
%   'sigma' - Parameter for Gaussian similarity function. Set
%      this to 0 for an unweighted graph. Default is 1.
%
%   Author: Ingo Buerk
%   Year  : 2011/2012
%   Bachelor Thesis

if nargin < 3
   ME = MException('InvalidCall:NotEnoughArguments', ...
       'Function called with too few arguments');
   throw(ME);
end

if ~any(Type == (1:2))
   ME = MException('InvalidCall:UnknownType', ...
       'Unknown similarity graph type');
   throw(ME);
end

n = size(M, 2);

% Preallocate memory
indi = zeros(1, k * n);
indj = zeros(1, k * n);
inds = zeros(1, k * n);

for ii = 1:n
    % Compute i-th column of distance matrix
    dist = distEuclidean(repmat(M(:, ii), 1, n), M);

    % Sort row by distance
    [s, O] = sort(dist, 'ascend');

    % Save indices and value of the k
    indi(1, (ii-1)*k+1:ii*k) = ii;
    indj(1, (ii-1)*k+1:ii*k) = O(1:k);
    inds(1, (ii-1)*k+1:ii*k) = s(1:k);
end

% Create sparse matrix
W = sparse(indi, indj, inds, n, n);

clear indi indj inds dist s O;

% Construct either normal or mutual graph
if Type == 1
    % Normal
    W = max(W, W');
else
    % Mutual
    W = min(W, W');
end

if nargin < 4 || isempty(sigma)
    sigma = 1;
end

% Unweighted graph
if sigma == 0
    W = (W ~= 0);

% Gaussian similarity function
elseif isnumeric(sigma)
    W = spfun(@(W) (simGaussian(W, sigma)), W);

else
    ME = MException('InvalidArgument:NotANumber', ...
        'Parameter epsilon is not numeric');
    throw(ME);
end

end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function cout = surrogate3(cin)

%tstoolbox/@core/surrogate3
%   Syntax:
%     * cout = surrogate3(cin)
%
%   Input Arguments:
%     * cin - core object
%
%   create surrogate data for a scalar time series by permuting samples
%   randomly
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

if ndim(cin) > 1
	error('not a scalar time series')
end

y = data(cin);

len = dlens(cin,1);

cout = core(y(randperm(len)));

author: benfulcher

repo: hctsa


text:
function event = read_micromed_event(trcfile)

% reads the events of the Micromed TRC format files

fid=fopen_or_error(trcfile,'r');

%------------------reading patient & recording info----------
fseek(fid,64,-1);
surname=char(fread(fid,22,'char'))';
name=char(fread(fid,20,'char'))';

fseek(fid,128,-1);
day=fread(fid,1,'char');
if length(num2str(day))<2
    day=['0' num2str(day)];
else
    day=num2str(day);
end
month=fread(fid,1,'char');
switch month
case 1
    month='JAN';
case 2
    month='FEB';
case 3
    month='MAR';
case 4
    month='APR';
case 5
    month='MAY';
case 6
    month='JUN';
case 7
    month='JUL';
case 8
    month='AUG';
case 9
    month='SEP';
case 10
    month='OCT';
case 11
    month='NOV';
case 12
    month='DEC';
end
year=num2str(fread(fid,1,'char')+1900);

%------------------ Reading Header Info ---------

fseek(fid,175,-1);
Header_Type=fread(fid,1,'char');
if Header_Type ~= 4
    ft_error('*.trc file is not Micromed System98 Header type 4')
end

fseek(fid,138,-1);
Data_Start_Offset=fread(fid,1,'uint32');
Num_Chan=fread(fid,1,'uint16');
Multiplexer=fread(fid,1,'uint16');
Rate_Min=fread(fid,1,'uint16');
Bytes=fread(fid,1,'uint16');
fseek(fid,176+8,-1);
Code_Area=fread(fid,1,'uint32');
Code_Area_Length=fread(fid,1,'uint32');
fseek(fid,192+8,-1);
Electrode_Area=fread(fid,1,'uint32');
Electrode_Area_Length=fread(fid,1,'uint32');

fseek(fid,400+8,-1);
Trigger_Area=fread(fid,1,'uint32');
Tigger_Area_Length=fread(fid,1,'uint32');

%----------------- Read Trace Data ----------

fseek(fid,Data_Start_Offset,-1);
switch Bytes
case 1
    trace=fread(fid,'uint8');
case 2
    trace=fread(fid,'uint16');
case 4
    trace=fread(fid,'uint32');
end
m=length(trace);
if rem(m,Num_Chan)~=0
    roundata=floor(m/Num_Chan);
    trace=trace(1:roundata*Num_Chan);
    m=length(trace);
end
trace=reshape(trace,Num_Chan,m/Num_Chan);

%---------------- Reading Trigger Data ----------
fseek(fid,Trigger_Area,-1);
for l=1:Tigger_Area_Length/6
    trigger(1,l)=fread(fid,1,'uint32');
    trigger(2,l)=fread(fid,1,'uint16');
end

first_trigger=trigger(1,1);
m=length(trace);
tl=length(trigger);
NoTrig=0;
for tr=1:tl
    if ((trigger(1,tr) <= m) && (trigger(1,tr) >= first_trigger))
        NoTrig=NoTrig+1;
    end
end
if NoTrig > 0
   	trigger=trigger(:,1:NoTrig);
else
	trigger=[];
	first_trigger=[];
end

fclose(fid);

event = [];

if ~isempty(trigger)
  for E=1:length(trigger)
    event(E).type    = 'MARKER';
    event(E).sample  = trigger(1,E)+1;
    event(E).value   = trigger(2,E);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
filename='Cantilever_triangle_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'MMA';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.4;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;


text:
function [coeff, score, latent, tsquare] = princomp_largedata(x,econFlag)
% Principal Components Analysis
%
% This is a version created for large data sets by Matthew Davidson
%
% The default Matlab PRINCOMP is naive to large data sets. Out-of-memory
% errors are easily obtained on imaging data. The solution is to replace
% concise but inefficient calls to repmat with loops and to eliminate
% large, unused variables.
%
% Tested on random data sets and produces identical output as original
% PRINCOMP. Speed penalty is drastic for small data. 50% slower on
% 50x1000 element data set, but on a 50x10000, only 2% slower.
%
% COEFF = PRINCOMP(X) performs principal components analysis on the N-by-P
% data matrix X, and returns the principal component coefficients, also
% known as loadings.  Rows of X correspond to observations, columns to
% variables.  COEFF is a P-by-P matrix, each column containing coefficients
% for one principal component.  The columns are in order of decreasing
% component variance.
%
% PRINCOMP centers X by subtracting off column means, but does not
% rescale the columns of X.  To perform PCA with standardized variables,
% i.e., based on correlations, use PRINCOMP(ZSCORE(X)).  To perform PCA
% directly on a covariance or correlation matrix, use PCACOV.
%
% [COEFF, SCORE] = PRINCOMP(X) returns the principal component scores,
% i.e., the representation of X in the principal component space.  Rows
% of SCORE correspond to observations, columns to components.
%
% [COEFF, SCORE, LATENT] = PRINCOMP(X) returns the principal component
% variances, i.e., the eigenvalues of the covariance matrix of X, in
% LATENT.
%
% [COEFF, SCORE, LATENT, TSQUARED] = PRINCOMP(X) returns Hotelling's
% T-squared statistic for each observation in X.
%
% When N <= P, SCORE(:,N:P) and LATENT(N:P) are necessarily zero, and the
% columns of COEFF(:,N:P) define directions that are orthogonal to X.
%
% [...] = PRINCOMP(X,'econ') returns only the elements of LATENT that are
% not necessarily zero, i.e., when N <= P, only the first N-1, and the
% corresponding columns of COEFF and SCORE.  This can be significantly
% faster when P >> N.
%
% :See Also: BARTTEST, BIPLOT, CANONCORR, FACTORAN, PCACOV, PCARES, ROTATEFACTORS.
%
% :References:
%   1. Jackson, J.E., A User's Guide to Principal Components,
%      Wiley, 1988.
%   2. Jolliffe, I.T. Principal Component Analysis, 2nd ed.,
%      Springer, 2002.
%   3. Krzanowski, W.J., Principles of Multivariate Analysis,
%      Oxford University Press, 1988.
%   4. Seber, G.A.F., Multivariate Observations, Wiley, 1984.
%
% ..
%   Copyright 1993-2005 The MathWorks, Inc.
%   $Revision: 2.9.2.9 $  $Date: 2006/10/02 16:35:01 $
% ..

% ..
%    When X has more variables than observations, the default behavior is to
%    return all the pc's, even those that have zero variance.  When econFlag
%    is 'econ', those will not be returned.
% ..

if nargin < 2, econFlag = 0; end

[n,p] = size(x);
if isempty(x)
    pOrZero = ~isequal(econFlag, 'econ') * p;
    coeff = zeros(p,pOrZero); coeff(1:p+1:end) = 1;
    score = zeros(n,pOrZero);
    latent = zeros(pOrZero,1);
    tsquare = zeros(n,1);
    return
end

% Center X by subtracting off column means
%x0 = x - repmat(mean(x,1),n,1);
for i=1:size(x, 2)
    x(:,i) = x(:,i) - mean(x(:,i));
end
r = min(n-1,p); % max possible rank of X0

% The principal component coefficients are the eigenvectors of
% S = X0'*X0./(n-1), but computed using SVD.
[U,sigma,coeff] = svd(x,econFlag); % put in 1/sqrt(n-1) later

if nargout < 2
    % When econFlag is 'econ', only (n-1) components should be returned.
    % See comment below.
    if (n <= p) && isequal(econFlag, 'econ')
        coeff(:,n) = [];
    end

else
    % Project X0 onto the principal component axes to get the scores.
    if n == 1 % sigma might have only 1 row
        sigma = sigma(1);
    else
        sigma = diag(sigma);
    end
    score = U .* repmat(sigma',n,1); % == x*coeff
    sigma = sigma ./ sqrt(n-1);

    % When X has at least as many variables as observations, eigenvalues
    % n:p of S are exactly zero.
    if n <= p
        % When econFlag is 'econ', nothing corresponding to the zero
        % eigenvalues should be returned.  svd(,'econ') won't have
        % returned anything corresponding to components (n+1):p, so we
        % just have to cut off the n-th component.
        if isequal(econFlag, 'econ')
            sigma(n,:) = []; % make sure this shrinks as a column
            coeff(:,n) = [];
            score(:,n) = [];

        % Otherwise, set those eigenvalues and the corresponding scores to
        % exactly zero.  svd(,0) won't have returned columns of U
        % corresponding to components (n+1):p, need to fill those out.
        else
            sigma(n:p,1) = 0; % make sure this extends as a column
            score(:,n:p) = 0;
        end
    end

    % The variances of the pc's are the eigenvalues of S = X0'*X0./(n-1).
    latent = sigma.^2;

    % Hotelling's T-squared statistic is the sum of squares of the
    % standardized scores, i.e., Mahalanobis distances.  When X appears to
    % have column rank < r, ignore components that are orthogonal to the
    % data.
    if nargout == 4
        if n > 1
            q = sum(sigma > max(n,p).*eps(sigma(1)));
            if q < r
                warning('stats:princomp:colRankDefX', ...
                        ['Columns of X are linearly dependent to within machine precision.\n' ...
                         'Using only the first %d components to compute TSQUARED.'],q);
            end
        else
            q = 0;
        end
        tsquare = (n-1) .* sum(U(:,1:q).^2,2); % == sum((score*diag(1./sigma)).^2,2)
    end
end


text:
function [y,ty]=v_correlogram(x,inc,nw,nlag,m,fs)
%V_CORRELOGRAM calculate correlogram [y,ty]=(x,inc,nw,nlag,m,fs)
% Usage:
%        do_env=1; do_hp=1;                            % flags to control options
%        [b,a,fx,bx,gd]=v_gammabank(25,fs,'',[80 5000]); % determine v_filterbank
%        y=v_filterbank(b,a,s,gd);                       % filter signal s
%        if do_env
%            [bl,al]=butter(3,2*800/fs);
%            y=filter(bl,al,v_teager(y,1),[],1);           % low pass envelope @ 800 Hz
%        end
%        if do_hp
%            y=filter(fir1(round(16e-3*fs),2*64/fs,'high'),1,y,[],1);  % HP filter @ 64 Hz
%        end
%        v_correlogram(y,round(10e-3*fs),round(16e-3*fs),round(12e-3*fs),'',fs);
%
% Inputs:
%        x(*,chan)  is the output of a filterbank (e.g. v_filterbank)
%                   with one column per filter channel
%        inc        frame increment in samples
%        nw         window length in samples [or window function]
%        nlag       number of lags to calculate
%        m          mode string:
%               'h' = Hamming window
%        fs         sample freq (affects only plots)
%
% Outputs:
%        y(lag,chan,frame) is v_correlogram. Lags are 1:nlag samples
%        ty                is time of the window energy centre for each frame
%                            x(n) is at time n
%
% For each channel, the calculated correlation for frame n comprises
%       y(t+1,*,n+1)=(win.*x(n*inc+(1:nw))'*x(n*inc+t+(1:nw))/sqrt(win'*abs(x(n*inc+(1:nw))).^2 * win'*abs(x(n*inc+t+(1:nw))).^2)
% This corresponds to the expression in (1.7) of [1] but incorporating the normalization from (1) of [2].
%
% Future planned mode options:
%       'd' = subtract DC component
%       'n' = unnormalized
%       'v' = variable analysis window proportional to lag
%       'p' = output the peaks only
%
% Refs:
% [1]	D. Wang and G. J. Brown. Fundamentals of computational auditory scene analysis.
%       In D. Wang and G. Brown, editors, Computational Auditory Scene Analysis: Principles,
%       Algorithms, and Applications, chapter 1. Wiley, Oct. 2006. doi: 10.1109/9780470043387.ch1
% [2]	S. Granqvist and B. Hammarberg. The correlogram: a visual display of periodicity. J. Acoust. Soc. Amer., 114: 2934, 2003.

%      Copyright (C) Mike Brookes 2011-2018
%      Version: $Id: v_correlogram.m 10867 2018-09-21 17:35:59Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
memsize=v_voicebox('memsize'); 	% set memory size to use
[nx,nc]=size(x);                % number of sampes and channels
if nargin<6
    fs=1;                       % default sample frequency is 1 Hz
    if nargin<5
        m='h';                  % default analysis window is Hamming
        if nargin<4
            nlag=[];
            if nargin<3
                nw=[];
                if nargin<2
                    inc=[];
                end
            end
        end
    end
end
if ~numel(inc)
    inc=128;                    % default frame hop is 128 samples
end
if ~numel(nw)
    nw=inc;                     % default analysis window length is the fame increment
end
nwin=length(nw);
if nwin>1                       % nw specifies the window function explicitly
    win=nw(:);
else                            % nw gives the window length
    nwin=nw;
    if any(m=='h')
        win=v_windows(3,nwin)'; % Hamming window
    else
        win=ones(nwin,1);       % window function
    end
end
if ~numel(nlag)
    nlag=nwin;
end
nwl=nwin+nlag-1;
nt=pow2(nextpow2(nwl));         % transform length
nf=floor((nx-nwl+inc)/inc);     % number of frames
i1=repmat((1:nwl)',1,nc)+repmat(0:nx:nx*nc-1,nwl,1);
nb=min(nf,max(1,floor(memsize/(16*nwl*nc))));    % chunk size for calculating
nl=ceil(nf/nb);                  % number of chunks
jx0=nf-(nl-1)*nb;                % size of first chunk in frames
wincg=(1:nwin)*win.^2/(win'*win); % determine window centre of energy
fwin=conj(fft(win,nt,1));       % conjugate fft of window
y=zeros(nlag,nc,nf);
% first do partial chunk
jx=jx0;
x2=zeros(nwl,nc*jx);
x2(:)=x(repmat(i1(:),1,jx)+repmat((0:jx-1)*inc,nwl*nc,1));
% the next line was previously: v=ifft(conj(fft(x2(1:nwin,:),nt,1)).*fft(x2,nt,1));
v=ifft(conj(fft(x2(1:nwin,:).*repmat(win(:),1,nc*jx),nt,1)).*fft(x2,nt,1));                 % v(1:nlag,:) contains second half of xcorr(x2) output
w=max(real(ifft(fwin(:,ones(1,nc*jx)).*fft(x2.*conj(x2),nt,1))),0); % v(1:nlag,:) contains second half of xcorr(|x2|^2,win) output
w=sqrt(w(1:nlag,:).*w(ones(nlag,1),:));
if isreal(x)
    y(:,:,1:jx)=reshape(real(v(1:nlag,:))./w,nlag,nc,jx); % note: some values may be NaN is x=0 throughout the window
else
    y(:,:,1:jx)=reshape(conj(v(1:nlag,:))./w,nlag,nc,jx);
end
% now do remaining chunks
x2=zeros(nwl,nc*nb);
for il=2:nl
    ix=jx+1;            % first frame in this chunk
    jx=jx+nb;           % last frame in this chunk
    x2(:)=x(repmat(i1(:),1,nb)+repmat((ix-1:jx-1)*inc,nwl*nc,1));
    % the next line was previously: v=ifft(conj(fft(x2(1:nwin,:),nt,1)).*fft(x2,nt,1));
    v=ifft(conj(fft(x2(1:nwin,:).*repmat(win(:),1,nc*nb),nt,1)).*fft(x2,nt,1));
    w=max(real(ifft(fwin(:,ones(1,nc*nb)).*fft(x2.*conj(x2),nt,1))),0);
    w=sqrt(w(1:nlag,:).*w(ones(nlag,1),:));
    if isreal(x)
        y(:,:,ix:jx)=reshape(real(v(1:nlag,:))./w,nlag,nc,nb);
    else
        y(:,:,ix:jx)=reshape(conj(v(1:nlag,:))./w,nlag,nc,nb);
    end
end
ty=(0:nf-1)'*inc+wincg;       % calculate times of window centres
if ~nargout
    imagesc(ty/fs,(1:nlag)/fs,squeeze(mean(abs(y),2)));
    if nargin<6
        us='samp';
    else
        us='s';
    end
    xlabel(['Time (' v_xticksi us ')']);
    ylabel(['Lag (' v_yticksi us ')']);
    axis 'xy';
    v_colormap('v_thermliny');
    colorbar;
    v_cblabel('Mean Correlation');
    title('Summary Correlogram');
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
% Matlab-Script: seisshift.m
% Script to evaluate increase and decrease of seismicity
% Author: J. Woessner
% woessner@seismo.ifg.ethz.ch
% last update: 04.09.02

% Track changes:
% 04.09.02 : Changed fcumulsum to calc_cumulsum

report_this_filefun(mfilename('fullpath'));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END Header %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Load data file %%%%%%%%
%%  Data is stored in catalog matrix a in the ZMAP form
%%% lon lat year month day mag depth hour min
mCatalog = ga; % Catalog initialazation
fSplitTime = mean(mCatalog(:,3));
fLastDate = max(mCatalog(:,3));
fFirstDate = min(mCatalog(:,3));

% [mCatalog1, mCatalog2, fPeriodExact1, fPeriodExact2, fPeriod, fPeriod2] = ...
%     ex_SplitCatalog(mCatalog, fSplitTime, 0, 0, 0, 0);
[mCatalog1, mCatalog2, fPeriodExact1, fPeriodExact2, fPeriod, fPeriod2] = ...
     ex_SplitCatalog(ga, 1984.7, 1, 1, 1,1);
maxmag = max(mCatalog1(:,6))
maxmag2 = max(mCatalog2(:,6))

[ev_val mags ev_valsum ev_valsum_rev,  mags_rev] =calc_cumulsum(mCatalog1);
[ev_val2 mags2 ev_valsum2 ev_valsum_rev2,  mags_rev2] =calc_cumulsum(mCatalog2);

% Figures
if exist('h_mag_fig','var') & ishandle(h_mag_fig)
    set(0,'Currentfigure',h_mag_fig);
    disp('Figure exists');
else
    h_mag_fig=figure_w_normalized_uicontrolunits('tag','ev_hist','Name','Magnitude histogram','Units','normalized','Nextplot','add',...
        'Numbertitle','off');
    h_mag_axs=axes('tag','ax_ev_hist','Nextplot','add','box','on');
end

xmag=max([maxmag maxmag2]); % Define length of x-axis

% Figure: Magnitude historgrams

subplot(2,1,1);
set(gca,'tag','ax_ev_hist1','Nextplot','replace','box','on');
axs1=findobj('tag','ax_ev_hist1');
axes(axs1(1));
histogram(mCatalog1(:,6),mags)
axis([0 xmag 0 max(ev_val)]);
ylabel('Events');
subplot(2,1,2);
set(gca,'tag','ax_ev_hist2','Nextplot','replace','box','on');
axs2=findobj('tag','ax_ev_hist2');
axes(axs2(1));
histogram(mCatalog2(:,6),mags2);
fig_hd1 = findobj(gca,'Type','patch');
set(fig_hd1,'FaceColor',[0 0.5 0])
axis([0 xmag 0 max(ev_val2)]);
ylabel('Events');
xlabel('Magnitude')

if exist('cum_mag_fig','var') &  ishandle(cum_mag_fig)
   set(0,'Currentfigure',cum_mag_fig);
   disp('Figure exists');
else
    cum_mag_fig=figure_w_normalized_uicontrolunits('tag','cumFMD','Name','Cumulative FMD','Units','normalized','Nextplot','add','Numbertitle','off');
    cum_mag_axs=axes('tag','ax_cumFMD','Nextplot','add','box','on');
end
% Figure: Bar and cumulative sums

subplot(2,1,1);
set(gca,'tag','ax_cumFMD1','Nextplot','replace','box','on');
axs3=findobj('tag','ax_cumFMD1')
axes(axs3(1));
plot(mags,ev_val,'-o',mags2,ev_val2,'-*')
ylabel('Non-cumulative number of events ');
subplot(2,1,2);
set(gca,'tag','ax_cumFMD2','Nextplot','replace','box','on','Yscale','log');
axs4=findobj('tag','ax_cumFMD2')
axes(axs4(1));
semilogy(mags,ev_valsum,'-o','Color',[0 0 1])
hold on;
semilogy(mags_rev,ev_valsum_rev,'-o','Color',[0 0 1])
semilogy(mags2,ev_valsum2,'-*','Color',[0 0.5 0])
semilogy(mags_rev2,ev_valsum_rev2,'-*','Color',[0 0.5 0])
ylabel('Cumulative sum');
xlabel('Magnitude')
hold off;

%% Difference illustrations
[params] = calc_totdiff(mCatalog, fSplitTime);

params.dNdiffsum
params.dNdiffsumYearVal
params.dNdiffsumMonthVal

if exist('hDiff_fig','var') &  ishandle(hDiff_fig)
   set(0,'Currentfigure',hDiff_fig);
   disp('Figure exists');
else
    hDiff_fig=figure_w_normalized_uicontrolunits('tag','hDiff','Name','Normalized differences','Units','normalized','Nextplot','add','Numbertitle','off');
    hDiff_axs=axes('tag','ax_hDiff','Nextplot','add','box','on');
end

subplot(3,1,1);
set(gca,'tag','ax_hDiff1','Nextplot','replace','box','on');
axs5=findobj('tag','ax_hDiff1')
axes(axs5(1));
bar(params.dMags,params.dNdiff);
ylabel('Seismicity difference');

subplot(3,1,2);
set(gca,'tag','ax_hDiff2','Nextplot','replace','box','on');
axs6=findobj('tag','ax_hDiff2')
axes(axs6(1));
bar(params.dMags,params.dNdiffYear)
ylabel(' Normalized per year');

subplot(3,1,3);
set(gca,'tag','ax_hDiff3','Nextplot','replace','box','on');
axs7=findobj('tag','ax_hDiff3')
axes(axs7(1));
bar(params.dMags,params.dNdiffMonth)
ylabel(' Normalized per month');

author: CelsoReyes


text:
EulerShockDriver2D
figure; PlotField2D(N, x, y, Q(:,:,4)); drawnow; pause(.1);

repo: Aero-matlab