text:
function out = isreal(f)
%ISREAL   Test if a CHEBFUN2V object F is real-valued.
%   ISREAL(F) returns logical true if F does not have any nonzero imaginary
%   part and false otherwise.
%
%   (This is slightly different from the Matlab convention, where isreal(x)
%   is false if x is a complex number whose imaginary part is 0.)
%
% See also CHEBFUN/ISREAL, CHEBFUN2/ISREAL, and CHEBFUN3/ISREAL.

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

if ( isempty(f) )
    out = true;
    return
end

for ii = 1:f.nComponents
    out(ii) = isreal(f.components{ii});
end
out = all(out);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function ye = pbd(t, x, y, side)
%PBD	Prescribes pressure p at outflow boundaries
% Possible values for side: 'lower', 'upper',  'left', 'right'

global geval Re

ye = 0;
if geval == 1|geval == 7		% Horizontal Poiseuille flow
  if strcmp(side, 'right') == 1		% Outflow
    ye = 0;
  else
    ye = 0;
  end
elseif geval == 2	% Vertical Poiseuille flow
  if strcmp(side, 'upper') == 1		% Outflow
    ye = 0;
  else
    ye = 0;
  end
elseif geval == 3	% Backward facing step
    ye = 0;
elseif geval == 4	% Driven cavity
    ye = 0;
elseif geval == 5|geval == 6	% Uniform flow under angle alpha
  ye = 1;
else
  error('Wrong value in input for parameter geval')
end

repo: Aero-matlab


text:
function [L,filters] = gaussSmooth( I, sigmas, shape, radius )
% Applies Gaussian smoothing to a (multidimensional) image.
%
% Smooths the n-dimensional array I with a n-dimensional gaussian with
% standard deviations specified by sigmas.  This operation in linearly
% seperable and is implemented as such.
%
% USAGE
%  [L,filters] = gaussSmooth( I, sigmas, [shape], [radius] )
%
% INPUTS
%  I       - input image
%  sigmas  - either n dimensional or 1 dimensional vector of standard devs
%            if sigmas(n)<=.3 then does not smooth along that dimension
%  shape   - ['full'] shape flag 'valid', 'full', 'same', or 'smooth'
%  radius  - [2.25] radius in units of standard deviation
%
% OUTPUTS
%  L       - smoothed image
%  filters - actual filters used, cell array of length n
%
% EXAMPLE
%  load trees; I=ind2gray(X,map);
%  I2 = gaussSmooth( I, 1, 'same' );
%  figure(1); im(I); figure(2); im(I2);
%
% See also FILTERGAUSS
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2008 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Lesser GPL [see external/lgpl.txt]

nd = ndims(I);  if(length(sigmas)==1); sigmas=repmat(sigmas,[1,nd]); end
if( nd > length(sigmas)); error('Incorrect # of simgas specified'); end
sigmas = sigmas(1:nd);

if( isa( I, 'uint8' ) ); I = double(I); end
if( nargin<3 || isempty(shape) ); shape='full'; end
if( nargin<4 || isempty(radius) ); radius=2.25; end

% create and apply 1D gaussian masks along each dimension
L = I;  filters = cell(1,nd);
for i=1:nd
  if (sigmas(i)>.3)
    r = ceil( sigmas(i)*radius );
    f = filterGauss( 2*r+1, [], sigmas(i)^2 );
    f = permute( f, circshift(1:nd,[1,i-1]) );
    filters{i} = f;
    L = convnFast( L, f, shape );
  else
    filters{i} = 1;
  end
end

author: kristinbranson

repo: JAABA


text:
function outputImage = kfilter(this, filterType, applicationDimensions, varargin)
% Filters by multiplication of a specified window function in k-space
%
%   Y = MrDataNd()
%   outputImage = Y.kfilter(filterType, applicationDimensions)
%
% This is a method of class MrDataNd.
%
% IN
%   filterType  string of filter to be applied, possible values are
%               'hamming' (default)
%               'hanning'
%               'raised_cosine'
%   applicationDimensions
%               '2D' or '3D'
%                   default: '2D' for single slices, '3D' otherwise
%               '2D' performs the filter slice-wise with the same filter,
%               '3D' performs the filter for a 3D symmetric version of the
%                    filter
%   varargin
%               extra filter parameters, depending on the chosen filter
%               'raised_cosine'
%               kfilter('raised_cosine', '2D' or '3D', 'fractionFOV', 0.5, ...
%               'beta', 0.5)
%                   fractionFOV  - fraction of FOV (1-dim!) where filter
%                                  reaches half Maximum
%                                  default: 0.5
%                   beta         - roll-off factor between 0 and 1 for the
%                                  raised-cosine window
%                                  (0 giving a box-car function,
%                                  and 1 a cosine without plateau)
%                                  default: 0.5
%   doPlotFilter    true of false (default)
%                   if true, an extra plot is generated, showing the
%                   filter response in k-space alongside central x- and
%                   y-line profiles of the image
%
% OUT
%
% EXAMPLE
%   kfilter
%
%   See also MrDataNd

% Author:   Lars Kasper, based on code by Johanna Vannesjo for 1D filtering
% Created:  2018-11-07
% Copyright (C) 2018 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC 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/>.

defaults.fractionFOV = 0.5;
defaults.beta = 0.5;
defaults.doPlotFilter = false;

args = tapas_uniqc_propval(varargin, defaults);
tapas_uniqc_strip_fields(args);

is3D = ndims(this) >= 2;

if nargin < 2
    filterType = 'hamming';
end

if nargin < 3
    if is3D
        applicationDimensions = '3D';
    else
        applicationDimensions = '2D';
    end
end

switch filterType
    case 'raised_cosine'
        % from J. Vannesjo, utils/general/raised_cosine.m Recon5-6, IBT
        funFilter = @(x) tapas_uniqc_raised_cosine((1:x) - floor(x/2), ...
            1/(fractionFOV*x), beta);
    otherwise
        funFilter = str2func(filterType);
end

dimInfoFilter = this.dimInfo.copyobj();

% column vector * row vector = matrix coordinate-wise product
filterMatrix = reshape(funFilter(this.dimInfo.nSamples(1)), [],1)*...
    reshape(funFilter(this.dimInfo.nSamples(2)), 1, []);

if doPlotFilter
    filterProfile = reshape(funFilter(this.dimInfo.nSamples(1)),[],1);
    figure('Name', 'k-filter profile');
    plot(filterProfile);
    xlim([1,this.dimInfo.nSamples(1)]);
    hold all;
    kDataProfile = this.image2k.abs.data(:,round(this.dimInfo.nSamples(2)/2), ...
        round(this.dimInfo.nSamples(3)/2));
    plot(kDataProfile/max(kDataProfile));
    plot(kDataProfile/max(kDataProfile).*filterProfile);
    legend('kfilter', 'unfiltered kx-profile of central slice', 'kfiltered kx-profile');
end

% replicate same filter for all slices
if is3D
    filterMatrix = repmat(filterMatrix, 1, 1, this.dimInfo.nSamples(3));
    dimInfoFilter.remove_dims(4:dimInfoFilter.nDims);

    switch applicationDimensions
        case '3D'
            % create the filter in 3rd dimension by replicating in other 2
            % dims and multiplying with slice-replicated 2D-filter
            filterMatrixThirdDim = reshape(funFilter(this.dimInfo.nSamples(3)), 1, 1, []);
            filterMatrix = filterMatrix.*repmat(filterMatrixThirdDim,...
                dimInfoFilter.nSamples(1), dimInfoFilter.nSamples(2), 1);
        case '2D'
            % everything fine, we just replicated for all slices
    end

else
    dimInfoFilter.remove_dims(3:dimInfoFilter.nDims);
end

filterImage = MrImage(filterMatrix, 'dimInfo', dimInfoFilter);
outputImage = k2image(image2k(this, applicationDimensions).*filterImage, ...
    applicationDimensions);

author: translationalneuromodeling

repo: tapas


text:
function PDP=ieee802_11_model(sigma_tau,Ts)
% IEEE 802.11 channel model PDP generator
% Input:
%       sigma_tau  : RMS delay spread
%       Ts         : Sampling time
% Output:
%       PDP        : Power delay profile

lmax = ceil(10*sigma_tau/Ts);
sigma02=(1-exp(-Ts/sigma_tau))/(1-exp(-(lmax+1)*Ts/sigma_tau)); % (2.9)
l=0:lmax;
PDP = sigma02*exp(-l*Ts/sigma_tau); % (2.8)

author: LyricYang


text:
function test_suite = test_createRay
%testCreateRay  One-line description here, please.
%   output = testCreateRay(input)
%
%   Example
%   testCreateRay
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-04-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.
% Licensed under the terms of the LGPL, see the file "license.txt"

test_suite = functiontests(localfunctions);

function testCreateRay2Points(testCase) %#ok<*DEFNU>

p1 = [1 1];
p2 = [2 3];
ray = createRay(p1, p2);

testCase.assertEqual(p1, ray(1,1:2), 'AbsTol', .01);
testCase.assertEqual(p2-p1, ray(1,3:4), 'AbsTol', .01);

function testCreateRay2Arrays(testCase)

p1 = [1 1;1 1];
p2 = [2 3;2 4];
ray = createRay(p1, p2);

testCase.assertEqual(2, size(ray, 1), 'AbsTol', .01);
testCase.assertEqual(p1, ray(:,1:2), 'AbsTol', .01);
testCase.assertEqual(p2-p1, ray(:,3:4), 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% Exact solution of Riemann problem

% Theory in Section 10.2 of:

% 	P. Wesseling: Principles of Computational Fluid Dynamics
% 	Springer, Heidelberg, 2000 ISBN 3-5453-0. XII, 642 pp.
% See http://ta.twi.tudelft.nl/nw/users/wesseling/cfdbook.html

% This program is called by Euler schemes

% Functions called: f

global  PRL  CRL MACHLEFT  gamma  pleft  pright  rholeft  rhoright  uleft...
	uright  tend  lambda		% lambda = dt/dx

gammab = 1/(gamma - 1); gam1 = gamma-1;

% Assumed structure of exact solution
%
%    \         /      |con |       |s|
%     \   f   /       |tact|       |h|
% left \  a  /  state |disc| state |o| right
% state \ n /    2    |cont|   3   |c| state
%   1    \ /          |tinu|       |k|   4
%         |           |ity |       | |

PRL = pright/pleft;
cright = sqrt(gamma*pright/rhoright); cleft = sqrt(gamma*pleft/rholeft);
CRL = cright/cleft;
MACHLEFT = (uleft - uright)/cleft;

p34 = fzero('f',3);		% p34 = p3/p4
p3 = p34*pright; 	alpha = (gamma+1)/(gamma-1);
rho3 = rhoright*(1+alpha*p34)/(alpha+p34);
rho2 = rholeft*(p34*pright/pleft)^(1/gamma);
u2 = uleft-uright+(2/(gamma-1))*cleft*...
	(1-(p34*pright/pleft)^((gamma-1)/(2*gamma)));
c2 = sqrt(gamma*p3/rho2);
spos = 0.5 + ...		% Shock position
	tend*cright*sqrt((gamma-1)/(2*gamma) + (gamma+1)/(2*gamma)*p34)+...
	tend*uright;

conpos = 0.5 + u2*tend + tend*uright;	% Position of contact discontinuity
pos1 = 0.5 + (uleft - cleft)*tend;	% Start of expansion fan
pos2 = 0.5 + (u2+uright-c2)*tend;	% End of expansion fan
xx = 0:0.002:1;
pexact = zeros(size(xx)); uexact= zeros(size(xx)); rhoexact = zeros(size(xx));
machexact = zeros(size(xx));  cexact = zeros(size(xx));
for i = 1:length(xx)
  if xx(i) <= pos1
    pexact(i) = pleft;    rhoexact(i) = rholeft;
    uexact(i) = uleft;    cexact(i)   = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  elseif xx(i) <= pos2
    pexact(i) = pleft*(1+(pos1-xx(i))/(cleft*alpha*tend))^(2*gamma/(gamma-1));
    rhoexact(i) = rholeft*(1+(pos1-xx(i))/(cleft*alpha*tend))^(2/(gamma-1));
    uexact(i) = uleft + (2/(gamma+1))*(xx(i)-pos1)/tend;
    cexact(i) = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  elseif xx(i) <= conpos
    pexact(i) = p3;    	      rhoexact(i) = rho2;
    uexact(i) = u2+uright;    cexact(i)   = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  elseif xx(i) <= spos
    pexact(i) = p3;    rhoexact(i) = rho3;    uexact(i) = u2+uright;
    cexact(i) = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  else
    pexact(i) = pright;    rhoexact(i) = rhoright;
    uexact(i) = uright;    cexact(i)   = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  end
end
entroexact = log(pexact./rhoexact.^gamma);

subplot(2,3,1),     	plot(xx,rhoexact)
subplot(2,3,2),    	plot(xx,uexact)
subplot(2,3,3),    	plot(xx,pexact)
subplot(2,3,4),  	plot(xx,machexact)
subplot(2,3,5),     	plot(xx,entroexact)


text:
function Sen = prepareForActiveInit(Sen,Obs)

if strcmp(Obs.ltype(4:6),'Pnt')
    if Obs.vis
        c = getCell(Obs.exp.e,Sen.imGrid);
        Sen.imGrid.usedCell(c(2),c(1)) = true;
    end
end

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function y = geo_mean( varargin )

%GEO_MEAN   Geometric mean.
%   Y=GEO_MEAN(X), where X is a vector, computes the geometrix mean of X. If any
%   of the elements of X are negative, then Y=-Inf. Otherwise, it is equivalent
%   to Y=PROD(X).^(1/LENGTH(X)). All elements must be real.
%
%   For matrices, GEO_MEAN(X) is a row vector containing the geometric means of
%   the columns. For N-D arrays, GEO_MEAN(X) is an array of the geometric means
%   taken along the first non-singleton dimension of X.
%
%   GEO_MEAN(X,DIM) takes the geometric mean along the dimension DIM of X.
%
%   GEO_MEAN(X,DIM,W), where W is a vector of nonnegative integers, computes a
%   weighted geometric mean Y = PROD(X.^W)^(1/SUM(W)). This is more efficient
%   than replicating the values of X W times. Note that W must be a vector,
%   even if X is a matrix, and its length must be the same as SIZE(X,DIM).
%
%   Disciplined convex programming information:
%       GEO_MEAN is concave  and nondecreasing; therefore, when used in CVX
%       specifications, its argument must be concave.

%
% Check arguments
%

persistent P
if isempty( P ),
    P.map = cvx_remap( { 'real' ; 'concave' ; 'l_convex' ; 'l_concave' } );
    P.map = bsxfun( @and, P.map, ~cvx_remap( 'negative' ) );
    P.funcs = { @geo_mean_1, @geo_mean_2, @geo_mean_2, @geo_mean_2 };
    P.zero = 1;
    P.constant = 1;
    P.reduce = true;
    P.reverse = false;
    P.name = 'geo_mean';
    P.errargs = [];
    P.dimarg = 2;
end
[ sx, x, dim, w ] = cvx_get_dimension( varargin, 2 );
if ~isempty( w ),
    if ~( numel(w)==length(w) && isnumeric(w) && isreal(w) && all(w>=0) ),
        cvx_throw( 'Third argument must be a vector of nonnegative numbers.' );
    elseif ~any( w ),
        cvx_throw( 'The weight vector cannot be all zeros.')
    elseif numel( w ) ~= sx(dim),
        cvx_throw( 'Third argument must be a vector of length %d', sx(dim) );
    else
        w = w(:) / sum(w);
    end
end
y = cvx_reduce_op( P, x, dim, w );

function y = geo_mean_1( x, w )
[ nx, nv ] = size( x );
if isempty( w ),
    w = 1 / nx;
else
    w = repmat( w, [ 1, nv ] );
end
y = prod( x .^ w, 1 );

function y = geo_mean_2( x, w ) %#ok
[ nx, nv ] = size( x );
cvx_begin
    hypograph variable y(1,nv);
    { linearize(x), y } == geo_mean_cone( [nx,nv], 1, w, 'func' ); %#ok
    cvx_setnneg(y);
cvx_end

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

author: yu-jiang

repo: radpbook

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


text:
function y = mat2lpc(x, f)
%MAT2LPC Compresses a matrix using 1-D lossles predictive coding.
%   Y = MAT2LPC(X, F) encodes matrix X using 1-D lossless predictive
%   coding. A linear prediction of X is made based on the
%   coefficients in F. If F is omitted, F = 1 (for previous pixel
%   coding) is assumed. The prediction error is then computed and
%   output as encoded matrix Y.
%
%   See also LPC2MAT.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.4 $  $Date: 2003/06/18 12:38:39 $

error(nargchk(1, 2, nargin));       % Check input arguments
if nargin < 2                       % Set default filter if omitted
    f = 1;
end

x = double(x);                      % Ensure double for computations
[m, n] = size(x);                   % Get dimensions of input matrix
p = zeros(m, n);                    % Init linear prediction to 0
xs = x;     zc = zeros(m, 1);       % Prepare for input shift and pad

for j = 1:length(f)                 % For each filter coefficient ...
    xs = [zc xs(:, 1:end - 1)];     % Shift and zero pad x
    p = p + f(j) * xs;              % Form partial prediction sums
end

y = x - round(p);                   % Compute the prediction error

repo: weiyanmin


text:
function b = multihconj(a, dim)
%MULTIHCONJ  Hermitian conjugating arrays of matrices.
%    B = MULTIHCONJ(A) is equivalent to B = MULTIHCONJ(A, DIM), where
%    DIM = 1.
%
%    B = MULTIHCONJ(A, DIM) is equivalent to
%    B = PERMUTE(A, [1:DIM-1, DIM+1, DIM, DIM+2:NDIMS(A)]), where A is an
%    array containing N P-by-Q matrices along its dimensions DIM and DIM+1,
%    and B is an array containing the Q-by-P Hermitian conjugate (') of
%    those N matrices along the same dimensions. N = NUMEL(A) / (P*Q), i.e.
%    N is equal to the number of elements in A divided by the number of
%    elements in each matrix.
%
%
%    Example:
%       A 5-by-9-by-3-by-2 array may be considered to be a block array
%       containing ten 9-by-3 matrices along dimensions 2 and 3. In this
%       case, its size is so indicated:  5-by-(9-by-3)-by-2 or 5x(9x3)x2.
%       If A is ................ a 5x(9x3)x2 array of 9x3 matrices,
%       C = MULTIHCONJ(A, 2) is a 5x(3x9)x2 array of 3x9 matrices.
%
%    See also MULTITRANSP MULTIHERM.

% This file is part of Manopt: www.manopt.org.
% Original author: Hiroyuki Sato, April 27, 2015.
% Contributors:
% Change log:

    % Setting DIM if not supplied.
    if nargin == 1, dim = 1; end

    % Transposing
    b = multitransp(a, dim);

    %Conjugating
    b = conj(b);

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function o = cat(dr,varargin)
% Concatenate file_array objects.  The result is a non-simple object
% that can no longer be reshaped.
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: cat.m 1143 2008-02-07 19:33:33Z spm $

if dr>32 || dr<0, error('Unknown command option.'); end;
dr   = max(round(dr),1);
d    = ones(nargin-1,16);
tmp  = {};
dpos = 0;
for i=1:nargin-1,
    vi = varargin{i};
    if strcmp(class(vi),'file_array')
        sz                = size(vi);
        d(i,1:length(sz)) = sz;
        svi               = struct(vi);
        svi               = svi(:);
        for j=1:length(svi(:)),
            if length(svi(j).pos)<dr
                svi(j).pos((length(svi(j).pos)+1):dr) = 1;
            end
            svi(j).pos(dr)= svi(j).pos(dr) + dpos;
        end;
        dpos              = dpos + d(i,dr);
        tmp{i}            = svi;
    else
        error(['Conversion to file_array from ' class(vi) ' is not possible.']);
    end;
end;
if any(diff(d(:,[1:(dr-1) (dr+1):end]),1,1))
    error('All matrices on a row in the bracketed expression must have the same number of rows.');
else
    o = vertcat(tmp{:});
    o = class(o,'file_array');
end;

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function r=v_rotpl2ro(u,v,t)
%V_ROTPL2RO find matrix to rotate in the plane containing u and v r=[u,v,t]
% Inputs:
%
%     U(n,1) and V(n,1) define a plane in n-dimensional space
%     T is the rotation angle in radians from U towards V. If T
%       is omitted it will default to the angle between U and V
%
% Outputs:
%
%     R(n,n)   Rotation matrix

%
%      Copyright (C) Mike Brookes 2007-2018
%      Version: $Id: v_rotpl2ro.m 10865 2018-09-21 17:22:45Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

u=u(:);
    n=length(u);
v=v(:);
l=sqrt(u'*u);
if l==0, error('input u is a zero vector'); end
u=u/l;      % normalize
q=v-v'*u*u;        % q is orthogonal to x
l=sqrt(q'*q);
if l==0          % u and v are colinear or v=zero
    [m,i]=max(abs(u));
    q=zeros(n,1);
    q(1+mod(i(1),n))=1;  % choose next available dimension
    q=q-q'*u*u;  % q is orthogonal to x
    l=sqrt(q'*q);
end
q=q/l;          % normalize
if nargin<3
    [s,c]=v_atan2sc(v'*q,v'*u);
    r=eye(n)+(c-1)*(u*u'+q*q')+s*(q*u'-u*q');
else
    r=eye(n)+(cos(t)-1)*(u*u'+q*q')+sin(t)*(q*u'-u*q');
end
if ~nargout && n==3
    v_rotqr2ro(v_rotro2qr(r)); % plot a rotated cube
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function error = verifyquadquadpts(n)
%% VERIFYQUADQUADPTS verify the quadrature formula on [0,1]^2
%
% Author: Huayi Wei < huayiwei1984@gmail.com>
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if (n>10), n=10; end

a = 0;
b = 1;

c = 0;
d = 1;

l = 1;

% get quadrature points
[pts,weight] = quadquadpts(n);
nQuad = size(pts,1);
t1 = 0;
t2 = 0;
for p = 1:nQuad
    % quadrature points in the x-y coordinate
    px = pts(p,1)*a + (1-pts(p,1))*b;
    py = pts(p,2)*c + (1-pts(p,2))*d;

    t1 = t1 + weight(p)*f1(px,py, n);
    t2 = t2 + weight(p)*f2(px,py);
end
t1 = t1*l;
t2 = t2*l;
error(1) = abs(t1 - 1/(4*n*n));
error(2) = abs(t2 - (- cos(1)+1)*(-cos(1)+1));
end

function z = f1(x, y, n)
z = x.^(2*n-1).*y.^(2*n-1);
end

function z = f2(x,y)
z = sin(x)*sin(y);
end

repo: ifem


text:
function d2fdx2 = hessian(f,x)
% HESSIAN Hessian of scalar polynomial SDPVAR object
%
% J = HESSIAN(p)    Hessian w.r.t all variables in p
% J = HESSIAN(p,x)  Hessian w.r.t the SDPVAR variables x
%
% See also INT, JACOBIAN, LINEARIZE, SDISPLAY

if nargin==1
    if isa(f,'sdpvar')
        x = recover(depends(f));
    else
        x = 0;
    end
else
    if length(getvariables(x))<length(x)
        error('x should be a vector of scalar independant variables');
    end
end

if prod(size(f))>1
    error('Hessian only defined for scalars.')
end

if any(ismember([getvariables(f) depends(f)],yalmip('extvariables')))
    error('Hessian is only applicable to polynomial expressions');
end

if isa(f,'double')
    d2fdx2 = zeros(length(x));
    return
end

d2fdx2 = jacobian(jacobian(f,x)',x);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function tl = read_1_line_pinnacle(fid)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

tl = fgetl(fid);
if ischar(tl)
	idx = strfind(tl,'//');
	if ~isempty(idx)
		if idx(1) == 1
			tl = '';
		else
			tl = tl(1:(idx-1));
		end
	end
	tl = ddeblank(tl);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ bool, in,out ] = gsp_check_connectivity( G,param )
%GSP_CHECK_CONNECTIVITY Check if the graph G is aperiodic strongly connected
%   Usage: bool=gsp_check_connectivity( G );
%          bool=gsp_check_connectivity( L );
%          bool=gsp_check_connectivity( W );
%          [bool,in,out]=gsp_check_connectivity( ... );
%
%   Input parameters:
%       G,W,L: Graph, Laplacian matrix or Weight martrix
%       param: Optional parameters
%
%   Output parameters
%       bool: Boolean
%       in  : Nodes without any in connections
%       out : Nodes without any out connections
%
%   Test if each node have at least one in connection and one out
%   connection. If this simple test give good results, the function compute
%   the perron vector of G and test it. It might take some time.
%
%   *param* is an optional structure that contains the following field
%
%   * *param.verbose*: display parameter - 0 no log - 1 display the errors
%

% Date: 6 june 2013
% Author: Nathanael Perraudin

%TODO: Use a clever method

% Handle Input parameters
if nargin<1, error('Not enought inputs parameters!'); end

if nargin<2, param=struct; end

if ~isfield(param, 'verbose'), param.verbose = 0; end

if isstruct(G)
    % If the graph is undirected, use the other function...
    if G.directed == 0
       [bool, in] = gsp_check_connectivity_undirected(G,param);
       out = in;
       return;
    end
    A=G.W;
else
    A=G;
end

if ~gsp_isdirected(A)
     [bool, in] = gsp_check_connectivity_undirected(G,param);
      out = in;
     return;
end

warning('This code is really bad!')

% Number of vertex
N=length(A);

% Remove the diagonal
A=A-diag(diag(A));

% Check the connecivity -- simple
    bool=~boolean((sum(1.-(sum(A,1)>0))+sum(1.-(sum(A,2)>0))));
    in=find(1.-(sum(A,1)>0));
    out=(find(1.-(sum(A,2)>0)))';
if param.verbose
    fprintf('   ---   Test if the graph is strongly connected   ---\n');
end

if bool
% Check the connectivity -- harder

    % Compute the Probablility matrix
    P=A./repmat(sum(A,2),1,N);

    % Compute the perron vector of P
    [phi,max_eig_P] = eigs(P',2);
    % test if max_eig_P==1
    if abs(max_eig_P(1,1)-1)>10e3*eps;
        bool=0;
        if param.verbose
                fprintf('    The maximum eigenvalue of P is not 1. \n');
        end
    else
        if param.verbose
                fprintf('    The maximum eigenvalue of P is 1. \n');
        end
    end
    phi=phi(:,1);
    % Test if the perron vector is positive
    if sum( phi)<0;
        phi=-phi;
    end

    if sum(phi<=10e3*eps)
        bool=0;
        if param.verbose
                fprintf('    Null of negative entry in the perron vector. \n');
        end
    else
        if param.verbose
                fprintf('    Stricly positive perron vector. \n');
        end
    end

    % see code from undirected for comments about this line
    max_eig_P = eigs(speye(size(P))+P,2);
    if abs(max_eig_P(2,2)-2)<10e3*eps;
        bool=0;
        if param.verbose
                fprintf('    Second eigenvalue of P is 1. \n');
        end
    else
        if param.verbose
                fprintf('    Only one unit eigenvalue of P.\n');
        end
    end

elseif param.verbose;
    fprintf('     Not every node has an input connection and an output connection! \n');
    fprintf('         No in connections for nodes: %i\n',in );
    fprintf('         No out connections for nodes: %i\n', out);
end

if param.verbose
    fprintf('   ---   End of tests   ---\n');
    if bool
        fprintf('    The graph is aperiodic strongly connected.\n')
    else
        fprintf('    The graph is not aperiodic stronly connected.\n');
    end
end

end

repo: gspbox


text:
function dap=links2dap(links);
% LINKS2DAP extract OpenDAP Data URLs from OpenSearch links
% Usage: dap=links2dap(links);
%    [links,struc]=opensearch(q);
%    dap=links2dap(links);

% Rich Signell (rsignell@usgs.gov)
k=0;
dap=[];
for i=1:length(links)
    s=char(links{i});
    if ~isempty(regexpi(s,'dodsC')),
        k=k+1;
%        e=regexpi(s,'\.html');
%        dap{k}=s(1:e-1);
        dap{k}=s;
    end
end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
%% with mandrill
x = mandrill;

wavelet = wavelet_factory_2d(size(x));
Sx = scat(x, wavelet);
sx = format_scat(Sx);

%% with 32,32
x = rand(32,32);

filt_opt.J = 3;
filt_opt.L = 6;
scat_opt.oversampling = 2;
[wavelet, filters] = wavelet_factory_2d(size(x), filt_opt, scat_opt);

Sx = scat(x, wavelet);
[s , meta] = format_scat(Sx);

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [text, timeStamp] = ReadPraatTextGrid(transFile)
feature('DefaultCharacterSet', 'UTF8');

trans = my_cat(transFile);

for i=1:length(trans)
    if length(trans{i})>4 && strcmp(trans{i}(1:4), 'item')
        break;
    end
end
i = i+ 6;
idx = regexp(trans{i}, '=');
nInterval = str2num(trans{i}(idx+1:end));

i = i + 1;
for cnt = 1:nInterval
    idx = regexp(trans{i+1}, '=');
    timeStamp(1,cnt) = str2num(trans{i+1}(idx+1:end));

    idx = regexp(trans{i+2}, '=');
    timeStamp(2,cnt) = str2num(trans{i+2}(idx+1:end));

    idx = regexp(trans{i+3}, '"');
    text{cnt} = trans{i+3}(idx(1)+1:idx(end)-1);
    i = i+4;

%      fprintf('%2.2f %2.2f %s\n', timeStamp(1,cnt), timeStamp(2,cnt), text{cnt});
end

cnt = 1;
text2 = {};
for j=1:length(text)
    if length(text{j})>0
        text2{cnt} = text{j};
        timeStamp2(:,cnt) = timeStamp(:,j);
        cnt = cnt + 1;
    end
end

text = text2;
timeStamp = timeStamp2;

author: singaxiong

repo: SignalGraph


text:
function e1_test ( )

%*****************************************************************************80
%
%% E1_TEST tests R4_E1 and R8_E1.
%
%  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, 'E1_TEST:\n' );
  fprintf ( 1, '  Test E1_VALUES, R4_E1, R8_E1.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X           E1(X)\n' );
  fprintf ( 1, '                      R4_E1(X)        Diff\n' );
  fprintf ( 1, '                      R8_E1(X)        Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx1 ] = e1_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_e1 ( single ( x ) );
    fx3 = r8_e1 ( x );

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

  end

  rmpath ( '../test_values' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [a,b,m,c,r,bbx,bby,f,g] = stokes_q1p0(xy,xyp,mv,ev)
%stokes_q1p0  vectorized Q1-P0 matrix generator
%   [A,B,Q,C,G,Bx,By,f,g] = stokes_q1p0(xy,xyp,mv,ev);
%   input
%          xy         Q2 nodal coordinate vector
%          xyp        Q1 nodal coordinate vector
%          mv         Q2 element mapping matrix
%          ev         element mapping matrix
%   output
%          A          Q1 vector diffusion matrix
%          B          Q1-P0 divergence matrix
%          Q          P0 mass matrix
%          C          pressure stabilization matrix
%          G          Q1 vector mass matrix
%          Bx         Q1 x-derivative matrix
%          By         Q1 y-derivative matrix
%          f          velocity rhs vector
%          g          pressure rhs vector
%
%   Natural boundary conditions apply. Dirichlet conditions
%   must be explicitly enforced by calling function flowbc.
%   IFISS function: DJS; 8 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
nngpt=4;
x=xy(:,1); y=xy(:,2);
xp=xyp(:,1); yp=xyp(:,2);
nvtx=length(x); nu=2*nvtx; np=length(xp);
nel=length(ev(:,1)); mp=[1:nel]';
lx=max(x)-min(x); ly=max(y)-min(y);
hx=max(diff(x)); hy=max(diff(y));
fprintf('setting up Q1-P0 matrices...  ')
%
% initialise global matrices
  a = sparse(nu,nu);
  r = sparse(nu,nu);
bbx = sparse(nvtx,nvtx);
bby = sparse(nvtx,nvtx);
 bx = sparse(np,nvtx);
 by = sparse(np,nvtx);
  b = sparse(np,nu);
  m = sparse(np,np);
  f = zeros(nu,1);
  g = zeros(np,1);
%
%
% Gauss point integration rules
if (nngpt==4)        % 2x2 Gauss points
   gpt=1.0e0/sqrt(3.0e0);
   s(1) = -gpt; t(1) = -gpt; wt(1)=1;
   s(2) =  gpt; t(2) = -gpt; wt(2)=1;
   s(3) =  gpt; t(3) =  gpt; wt(3)=1;
   s(4) = -gpt; t(4) =  gpt; wt(4)=1;
elseif (nngpt==1)   % 1x1 Gauss point
   s(1) =    0; t(1) =    0; wt(1)=4;
else
   error('Check Gauss point integration specification')
end
%
% inner loop over elements
for ivtx = 1:4
   xl_v(:,ivtx) = x(ev(:,ivtx));
   yl_v(:,ivtx) = y(ev(:,ivtx));
end
  ae = zeros(nel,4,4);
  re = zeros(nel,4,4);
bbxe = zeros(nel,4,4);
bbye = zeros(nel,4,4);
 bxe = zeros(nel,1,4);
 bye = zeros(nel,1,4);
 mpe = zeros(nel,1,1);
%
% loop over Gauss points
for igpt = 1:nngpt
   sigpt=s(igpt);
   tigpt=t(igpt);
   wght=wt(igpt);
%  evaluate derivatives etc
   [jac,invjac,phi,dphidx,dphidy] = deriv(sigpt,tigpt,xl_v,yl_v);
   for j = 1:4
      for i = 1:4
         ae(:,i,j)  = ae(:,i,j)  + wght*dphidx(:,i).*dphidx(:,j).*invjac(:);
         ae(:,i,j)  = ae(:,i,j)  + wght*dphidy(:,i).*dphidy(:,j).*invjac(:);
         re(:,i,j)  = re(:,i,j)  + wght*phi(:,i).*phi(:,j).*jac(:);
         bbxe(:,i,j) = bbxe(:,i,j) - wght*phi(:,i) .*dphidx(:,j);
         bbye(:,i,j) = bbye(:,i,j) - wght*phi(:,i) .*dphidy(:,j);
      end
      bxe(:,1,j) = bxe(:,1,j) - wght* dphidx(:,j);
      bye(:,1,j) = bye(:,1,j) - wght* dphidy(:,j);
   end
   mpe(:,1,1) = mpe(:,1,1) + wght*jac(:);
% end of Gauss point loop
end
%
% element assembly into global matrices
% component velocity matrices ...
for krow=1:4
   nrow=ev(:,krow);
   for kcol=1:4
      ncol=ev(:,kcol);
      a = a + sparse(nrow,ncol,ae(:,krow,kcol),nu,nu);
      a = a + sparse(nrow+nvtx,ncol+nvtx,ae(:,krow,kcol),nu,nu);
      r = r + sparse(nrow,ncol,re(:,krow,kcol),nu,nu);
      r = r + sparse(nrow+nvtx,ncol+nvtx,re(:,krow,kcol),nu,nu);
      bbx = bbx + sparse(nrow,ncol,bbxe(:,krow,kcol),nvtx,nvtx);
      bby = bby + sparse(nrow,ncol,bbye(:,krow,kcol),nvtx,nvtx);
   end
   kcol=1;
   ncol=mp;
   bx = bx + sparse(ncol,nrow,bxe(:,kcol,krow),np,nvtx);
   by = by + sparse(ncol,nrow,bye(:,kcol,krow),np,nvtx);
end
%
% vector velocity matrices ...
b = [bx,by];
%
% pressure matrices ...
m = sparse(mp,mp,mpe(:,1,1),np,np);
%
% stabilisation matrix
mel=length(mv(:,1));
cm=zeros(mel,4,4); hm=zeros(mel,1);
c=sparse(np,np);
% loop over macroelements
for ivtx = 1:9
   xlm(:,ivtx) = x(mv(:,ivtx));
   ylm(:,ivtx) = y(mv(:,ivtx));
end
elarea=full(diag(m));
hm(:)=sum(reshape(elarea,4,mel))';
hm=hm/4;
%      le=1/he;
%      cm=0.25*hm*[ le41+le12,     -le12,         0,     -le41;
%  	               -le12, le12+le23,     -le23,         0;
%                      0,     -le23, le23+le34,     -le34;
%	 			    -le41,         0,     -le34, le34+le41];
%
xe=xlm(:,9)-xlm(:,5);ye=ylm(:,9)-ylm(:,5);le12=hm;
cm(:,1,1) =cm(:,1,1)+ le12; cm(:,1,2)=-le12;
cm(:,2,2) =cm(:,2,2)+ le12; cm(:,2,1)=-le12;
xe=xlm(:,9)-xlm(:,6);ye=ylm(:,9)-ylm(:,6);le23=hm;
cm(:,2,2) =cm(:,2,2)+ le23; cm(:,2,3)=-le23;
cm(:,3,3) =cm(:,3,3)+ le23; cm(:,3,2)=-le23;
xe=xlm(:,9)-xlm(:,7);ye=ylm(:,9)-ylm(:,7);le34=hm;
cm(:,3,3) =cm(:,3,3)+ le34; cm(:,3,4)=-le34;
cm(:,4,4) =cm(:,4,4)+ le34; cm(:,4,3)=-le34;
xe=xlm(:,9)-xlm(:,8);ye=ylm(:,9)-ylm(:,8);le41=hm;
cm(:,4,4) =cm(:,4,4)+ le41; cm(:,4,1)=-le41;
cm(:,1,1) =cm(:,1,1)+ le41; cm(:,1,4)=-le41;
%
%  macroelement assembly into global matrices
for krow=1:4
   nrow=[0:4:nel-4]+krow;
   for kcol=1:4
      ncol=[0:4:nel-4]+kcol;
      c = c + sparse(nrow,ncol,cm(:,krow,kcol),np,np);
   end
end
%
fprintf('done\n')
return


text:
function binomTable=makeBinomTable(maxDeg)
%%MAKEBINOMTABLE  Create a lookup table of binomial coefficients.
%                 binomTable(n+1,k+1) is the value of binomial(n,k). This
%                 just build's Pascal's triangle. Consider, for example [1].
%                 The table uses about twice the memory it actually needs
%                 just so that it can be easily accessed as a matrix. This
%                 function is easier to use than the output of the pascal
%                 command in Matlab.
%
%INPUTS: maxDeg The maximum degree of the table. This is a value >=0.
%
%OUTPUTS: binomTable A (maxDeg+1)X(maxDeg+1) matrix where
%               binomTable(n+1,k+1) is the value of binomial(n,k). Unused
%               components are set to zero.
%
%This function can be useful when a large number of low-order binomial
%values are needed. The recursion used to fill in the table is from [1].
%
%REFERENCES:
%[1] Stover, Christopher and Weisstein, Eric W. "Pascal's Triangle." From
%    MathWorld--A Wolfram Web Resource.
%    http://mathworld.wolfram.com/PascalsTriangle.html
%
%September 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    binomTable=zeros(maxDeg+1,maxDeg+1);
    binomTable(:,0+1)=1;%Boundary values.
    for n=1:maxDeg
        for k=1:(n-1)
            binomTable(n+1,k+1)=binomTable(n-1+1,k-1+1)+binomTable(n-1+1,k+1);
        end
        binomTable(n+1,n+1)=1;%Boundary value
    end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [time_projection] = assembly_activity(AssemblyTemplates,SpikeCount)

% Activities = assembly_activity(Patterns,Activitymatrix): computes the
% time course of the activity of assembly patterns defined in Patterns in
% the spike matrix Activitymatrix with single bin resolution.
%
% Description of inputs:
%   Activitymatrix: spike matrix. Rows represent neurons, columns represent
%   time bins. Patterns: assembly patterns. Columns denote assembly # and
%   rows neuron #.
%
% Description of output:
%   Activities: Time course of the activity of assemblies. Rows represent
%   assemblies and columns represent time bins.
%
% This framework is described in: Lopes-dos-Santos V, Ribeiro S, Tort ABL
% (2013) Detecting cell assemblies in large neuronal populations, Journal
% of Neuroscience Methods.
%
% Please send bug reports to vitor@neuro.ufrn.br (V�tor)

zSpikeCount = zscore(SpikeCount')';

time_projection=zeros(size(AssemblyTemplates,2),size(zSpikeCount,2));
for assembly_idx = 1:size(AssemblyTemplates,2)

    % computing projector
    ASSEMBLYPROJECTOR=AssemblyTemplates(:,assembly_idx)*AssemblyTemplates(:,assembly_idx)';
    ASSEMBLYPROJECTOR=squeeze(ASSEMBLYPROJECTOR)-diag(diag(squeeze(ASSEMBLYPROJECTOR)));

    % computing activity time course
    for ntime=1:size(zSpikeCount,2)

        time_projection(assembly_idx,ntime)=(zSpikeCount(:,ntime)'*ASSEMBLYPROJECTOR*zSpikeCount(:,ntime));

    end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
% ##2
%==============================================================================
%
%   - data                 PETCT, Omega=(0,140)x(0,151), level=4:7, m=[128,128]
%   - viewer               viewImage2D
%   - interpolation        splineInter
%   - distance             NGF
%   - pre-registration     rigid2D
%   - regularizer          mbCurvature
%   - optimization         Gauss-Newton
% ===============================================================================

close all, help(mfilename);

setup2DPETCTData
imgModel('reset','imgModel','splineInter','regularizer','moments','theta',1e-2);
distance('reset','distance','NGF','edge',25);
trafo('reset','trafo','rigid2D');
regularizer('reset','regularizer','mbCurvature','alpha',10);
[yc,wc,his] = MLIR(ML,...
  'minLevel',4,'maxIterNPIR',25,'parametric',1,'plotMLiter',0);

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


text:
function RDt(option, up)
%RDt detects QRS waves in ECG signals.
%	            RDt(option, up)
%
%	Inputs:
%		option          the option which has led to this function being used
%       up              universal parameters structure
%
%	Outputs:
%       ...
%

fprintf('\n--- Detecting QRS waves ');
log_int_respSig = 1;             % Has value 1 unless this is a final respiratory signal

for subj = up.paramSet.subj_list

    sig_type = option(1:3);
    %% Cycle through each signal of this type
    eval(['sigs = up.paramSet.' sig_type '_sigs;']);
    for sig_no = 1 : length(sigs)
        curr_sig = sigs{sig_no};

        %% Skip if this processing has been done previously
        iden_resp_sig_file_ending
        savepath = [up.paths.data_save_folder, num2str(subj), ending];
        filecontents = whos('-file', savepath);
        var_names = extractfield(filecontents, 'name');
        temp = ~cellfun(@isempty, strfind(var_names, [curr_sig, up.paths.filenames.elim_vhf]));
        rel_var_names = var_names(temp);
        for rel_var_name_no = 1 : length(rel_var_names)
            for current_opt_no = 1 : length(up.al.options.RDt)
                temp = strfind(rel_var_names{rel_var_name_no},'_');
                start_el = temp(1); clear temp
                eval(['save_name = ''' curr_sig, up.paths.filenames.qrss, up.al.options.RDt{current_opt_no}, rel_var_names{rel_var_name_no}(start_el:end) ''';']); clear start_el
                exist_log = check_exists(savepath, save_name);
                if exist_log
                    continue
                end

                %% Load relevant data
                loadpath = [up.paths.data_save_folder, num2str(subj), up.paths.filenames.int_respSigs];
                rel_name = rel_var_names{rel_var_name_no};
                load(loadpath, rel_name);
                eval(['s = ' rel_name ';']);
                eval(['fs = ' rel_name '.fs;']);
                s.v = s.v(:);

                %% Eliminate very low frequencies
                s_filt = elim_sub_cardiac(s, up);
                nan_els = isnan(s_filt.v);
                s_filt.t = s_filt.t(~nan_els);
                s_filt.v = s_filt.v(~nan_els);

                %% Detect QRS Waves
                pk_inds = feval(up.al.options.RDt{current_opt_no}, s_filt, fs, subj, up);

                %% Refine detected QRS spikes
                % The detected peaks are currently somewhere in the middle of the R-area, so need to search for the max ecg value within the tolerance either side of the detected peak:
                max_HR_for_qrs = 300;
                tolerance = round(fs/(max_HR_for_qrs/60));
                if ~isempty(pk_inds)
                    % refine the peaks and store in a new variable (ECG_PKS_inds)
                    ref_pk_inds = nan(length(pk_inds),1);
                    for peak_ind = 1 : length(pk_inds)                             % for each peak
                        % Find tolerance limits, centred on the current peak index:
                        if pk_inds(peak_ind) <= tolerance                          % if the peak is right at the start of the recording, set the lower tolerance limit to the start of the recording
                            lower_lim = 1;
                        else
                            lower_lim = pk_inds(peak_ind)-tolerance;
                        end
                        if (length(s.v)-pk_inds(peak_ind)) <= tolerance % if the peak is right at the end of the recording, set the upper tolerance limit to the end of the recording
                            upper_lim = length(s.v);
                        else
                            upper_lim = pk_inds(peak_ind)+tolerance;
                        end
                        % Find the maximum ecg value within the tolerance limits:
                        [~, max_ind] = max(s.v(lower_lim : upper_lim));
                        % Store the index of this maximum value, referenced to the section_data inds:
                        ref_pk_inds(peak_ind) = lower_lim-1+max_ind;
                        clear max_ind
                    end
                    clear peak_ind
                else
                    % if no peaks were detected then give empty results:
                    ref_pk_inds = [];
                end

                %% Eliminate any peaks which are the same
                ref_pk_inds = unique(ref_pk_inds);

                %% Find troughs
                % search 0.1s either side
                troughs.i = nan(length(ref_pk_inds)-1,1);
                search_min = ref_pk_inds - ceil(0.1*fs);
                search_min(search_min<1) = 1;
                search_max = ref_pk_inds + ceil(0.1*fs);
                search_max(search_max>length(s.v)) = 1;
                for peak_ind = 1 : (length(ref_pk_inds)-1)
                    [~, rel_el] = min(s.v(search_min(peak_ind):search_max(peak_ind)));
                    rel_el = search_min(peak_ind)-1+rel_el;
                    troughs.i(peak_ind) = rel_el;
                end

                %% Save processed data
                eval([save_name '.fs = fs;']);
                eval([save_name '.p.v = s.v(ref_pk_inds);']);
                eval([save_name '.p.t = s.t(ref_pk_inds);']);
                eval([save_name '.tr.v = s.v(troughs.i);']);
                eval([save_name '.tr.t = s.t(troughs.i);']);
                % Identify start and end times of raw data (for resampling)
                eval([save_name '.timings.t_start = s.t(1);']);
                eval([save_name '.timings.t_end = s.t(end);']);
                save_or_append_data

            end

        end

    end

end

end

function pk_inds = GC(s, fs, subj, up)

%% This uses Prof Gari Clifford's "rpeakdetect.m" script
% This script can be downloaded from:
%    http://www.mit.edu/~gari/CODE/ECGtools/ecgBag/rpeakdetect.m
%
% The following is an excerpt from the script:
%
% Written by G. Clifford gari@ieee.org and made available under the
% GNU general public license. If you have not received a copy of this
% license, please download a copy from http://www.gnu.org/
%
% Please distribute (and modify) freely, commenting
% where you have added modifications.
% The author would appreciate correspondence regarding
% corrections, modifications, improvements etc.

% download the script if it isn't in the search path
curr_dir = mfilename('fullpath'); curr_dir = curr_dir(1:end-3);
filepath = [curr_dir, 'rpeakdetect.m'];
if ~exist(filepath, 'file')
    url = 'http://www.mit.edu/~gari/CODE/ECGtools/ecgBag/rpeakdetect.m';
    downloadedfilename = websave(filepath,url);
end

[~, ~, ~, pk_inds, ~, ~]  = rpeakdetect(s.v(:),fs);

end

author: peterhcharlton

repo: RRest


text:
function closescreen()
%CLOSESCREEN Dispose FULLSCREEN() window
%

global frame_java
try frame_java.dispose(); end

author: Sable

repo: mcbench-benchmarks


text:
%IBVS   Implement classical IBVS for point features
%
%  results = ibvs(T)
%  results = ibvs(T, params)
%
%  Simulate IBVS with for a square target comprising 4 points is placed
%  in the world XY plane. The camera/robot is initially at pose T and is
%  driven to the orgin.
%
%  Two windows are shown and animated:
%   1. The camera view, showing the desired view (*) and the
%      current view (o)
%   2. The external view, showing the target points and the camera
%
% The results structure contains time-history information about the image
% plane, camera pose, error, Jacobian condition number, error norm, image
% plane size and desired feature locations.
%
% The params structure can be used to override simulation defaults by
% providing elements, defaults in parentheses:
%
%   target_size    - the side length of the target in world units (0.5)
%   target_center  - center of the target in world coords (0,0,3)
%   niter          - the number of iterations to run the simulation (500)
%   eterm          - a stopping criteria on feature error norm (0)
%   lambda         - gain, can be scalar or diagonal 6x6 matrix (0.01)
%   ci             - camera intrinsic structure (camparam)
%   depth          - depth of points to use for Jacobian, scalar for
%                    all points, of 4-vector.  If null take actual value
%                    from simulation      ([])
%
% SEE ALSO: ibvsplot

% IMPLEMENTATION NOTE
%
% 1.  As per task function notation (Chaumette papers) the error is
%     defined as actual-demand, the reverse of normal control system
%     notation.
% 2.  The gain, lambda, is always positive
% 3.  The negative sign is written into the control law

classdef IBVS_sph < VisualServo

    properties
        pt          % phi, theta
        pt_star

        lambda          % IBVS gain
        eterm

        depth

        proj_lines
    end

    methods

        function ibvs = IBVS_sph(cam, varargin)

            % invoke superclass constructor
            ibvs = ibvs@VisualServo(cam, varargin{:});

            % handle arguments
            opt.eterm = 0.01;
            opt.lambda = 0.04;         % control gain
            opt.depth = [];
            opt.example = false;
            opt.T0 = SE3(0.3, 0.3, -2)*SE3.Rz(0.2);
            opt.Tf = SE3(0, 0, -1.5)*SE3.Rz(1);

            opt = tb_optparse(opt, ibvs.arglist);

            % copy options to IBVS object
            if opt.example
                % run a canned example
                fprintf('---------------------------------------------------\n');
                fprintf('canned example, spherical IBVS with 4 points\n');
                fprintf('---------------------------------------------------\n');
                ibvs.P = mkgrid(2, 1.5, 'pose', SE3(0,0,0.5));
                ibvs.Tf = SE3(0.3, 0.3, -2)*SE3.Rz(0.2);
                ibvs.T0 = SE3.convert(opt.T0);
                %ibvs.T0 = transl(-1,-0.1,-3);%*trotx(0.2);
            else
                ibvs.Tf = SE3(0, 0, -1.5)*SE3.Rz(1);
                ibvs.T0 = SE3.convert(opt.T0);
            end

            ibvs.lambda = opt.lambda;
            ibvs.eterm = opt.eterm;
            ibvs.depth = opt.depth;

            clf
            subplot(121);
            ibvs.camera.plot_create(gca)
            ibvs.camera.plot(ibvs.P, 'ro')

            % this is the 'external' view of the points and the camera
            subplot(122)
            plot_sphere(ibvs.P, 0.08, 'r');
            ibvs.camera.plot_camera('pose', ibvs.T0);
            plotvol([-2 2 -2 2 -3 1])
            view(16, 28);
            grid on
            set(gcf, 'Color', 'w')

            % draw lines from points to centre of sphere
            centre = ibvs.camera.T.t;
            ibvs.proj_lines = gobjects;
            for i=1:numcols(ibvs.P)
                ibvs.proj_lines(i) = plot3([centre(1) ibvs.P(1,i)], [centre(2) ibvs.P(2,i)], [centre(3) ibvs.P(3,i)], 'k')
            end

            set(gcf, 'HandleVisibility', 'Off');

            ibvs.type = 'spherical';
        end

        function init(vs)

            if isempty(vs.Tf)
                vs.Tf = transl(0, 0, -1);
                warning('setting Tf to default');
            end
            %% initialize the vservo variables
            vs.camera.T = vs.T0;    % set camera back to its initial pose
            vs.Tcam = vs.T0;                % initial camera/robot pose

            % final pose is specified in terms of a camera-target pose
            %   convert to image coords
            vs.pt_star = vs.camera.project(vs.P, 'pose', vs.Tf);

            vs.history = [];
        end

        function status = step(vs)
            status = 0;
            Zest = [];

                        % plot the points on mage plane
            vs.camera.clf();
            vs.camera.hold(true);
            vs.camera.plot(vs.P, 'pose', vs.Tcam, 'ro');
            vs.camera.plot(vs.pt_star, 'r*');

            % compute the view
            pt = vs.camera.project(vs.P, 'pose', vs.Tcam);

            % compute image plane error as a column
            e = pt - vs.pt_star;   % feature error
            e(2,:) = angdiff(e(2,:));
            e = e(:);

            % compute the Jacobian
            if isempty(vs.depth)
                % exact depth from simulation (not possible in practice)
                P_C = homtrans(inv(vs.Tcam), vs.P);
                J = vs.camera.visjac_p(pt, P_C(3,:) );
            else
                J = vs.camera.visjac_p(pt, vs.depth );
            end

            % compute the velocity of camera in camera frame
            try
                v = -vs.lambda * pinv(J) * e;
            catch
                status = -1;
                return
            end

            if vs.verbose
                fprintf('v: %.3f %.3f %.3f %.3f %.3f %.3f\n', v);
            end

            % update the camera pose
            Td = trnorm(delta2tr(v));    % differential motion

            vs.Tcam = vs.Tcam .* SE3(Td);       % apply it to current pose

            % update the camera pose
            vs.camera.T = vs.Tcam;

            % draw lines from points to centre of sphere
            centre = vs.Tcam.t;
            for i=1:numcols(vs.P)
                vs.proj_lines(i).XData = [centre(1) vs.P(1,i)];
                vs.proj_lines(i).YData = [centre(2) vs.P(2,i)];
                vs.proj_lines(i).ZData = [centre(3) vs.P(3,i)];
            end

            % update the history variables
            hist.pt = pt(:);
            vel = tr2delta(Td);
            hist.vel = vel;
            hist.e = e;
            hist.en = norm(e);
            hist.jcond = cond(J);
            hist.Tcam = vs.Tcam;

            vs.history = [vs.history hist];

            if norm(e) < vs.eterm
                status = 1;
                return
            end
        end

    end % methods
end % class

author: petercorke

repo: machinevision-toolbox-matlab


text:
function saveTightFigure(h,outfilename)
% SAVETIGHTFIGURE(OUTFILENAME) Saves the current figure without the white
%   space/margin around it to the file OUTFILENAME. Output file type is
%   determined by the extension of OUTFILENAME. All formats that are
%   supported by MATLAB's "saveas" are supported.
%
%   SAVETIGHTFIGURE(H, OUTFILENAME) Saves the figure with handle H.
%
% E Akbas (c) Aug 2010
% * Updated to handle subplots and multiple axes. March 2014.
%

if nargin==1
    hfig = gcf;
    outfilename = h;
else
    hfig = h;
end

%% find all the axes in the figure
hax = findall(hfig, 'type', 'axes');

%% compute the tighest box that includes all axes
tighest_box = [Inf Inf -Inf -Inf]; % left bottom right top
for i=1:length(hax)
    set(hax(i), 'units', 'centimeters');

    p = get(hax(i), 'position');
    ti = get(hax(i), 'tightinset');

    % get position as left, bottom, right, top
    p = [p(1) p(2) p(1)+p(3) p(2)+p(4)] + ti.*[-1 -1 1 1];

    tighest_box(1) = min(tighest_box(1), p(1));
    tighest_box(2) = min(tighest_box(2), p(2));
    tighest_box(3) = max(tighest_box(3), p(3));
    tighest_box(4) = max(tighest_box(4), p(4));
end

%% move all axes to left-bottom
for i=1:length(hax)
    if strcmp(get(hax(i),'tag'),'legend')
        continue
    end
    p = get(hax(i), 'position');
    set(hax(i), 'position', [p(1)-tighest_box(1) p(2)-tighest_box(2) p(3) p(4)]);
end

%% resize figure to fit tightly
set(hfig, 'units', 'centimeters');
p = get(hfig, 'position');

width = tighest_box(3)-tighest_box(1);
height =  tighest_box(4)-tighest_box(2);
set(hfig, 'position', [p(1) p(2) width height]);

%% set papersize
set(hfig,'PaperUnits','centimeters');
set(hfig,'PaperSize', [width height]);
set(hfig,'PaperPositionMode', 'manual');
set(hfig,'PaperPosition',[0 0 width height]);

%% save
saveas(hfig,outfilename);

author: peiyunh

repo: tiny

save_path: github-repos/MATLAB/peiyunh-tiny


text:
function [inputImageFiltered, additionalOutput] = fftImage(inputImage,varargin)
	% Computes FFT on input image.
	% Biafra Ahanonu
	% started: 2013.11.09
	% inputs
		% inputImage - [x y] matrix
	% outputs
		% inputImageFiltered - [x y] matrix
	% example
		% test the lowpass and highpass on an image with a range of options
		% f = fftImage(frame,'runfftTest',1,'bandpassType','lowpass');
		% f = fftImage(frame,'runfftTest',1,'bandpassType','highpass');

	[inputImageFiltered, additionalOutput] = ciapkg.image.fftImage(inputImage,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
%[2006]-"Ant Colony Optimization"

% (9/12/2020)

function ACS = jAntColonySystem(feat,label,opts)
% Parameters
tau   = 1;      % pheromone value
eta   = 1;      % heuristic desirability
alpha = 1;      % control pheromone
beta  = 1;      % control heuristic
rho   = 0.2;    % pheromone trail decay coefficient
phi   = 0.5;    % pheromena coefficient

if isfield(opts,'N'), N = opts.N; end
if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'tau'), tau = opts.tau; end
if isfield(opts,'alpha'), alpha = opts.alpha; end
if isfield(opts,'beta'), beta = opts.beta; end
if isfield(opts,'rho'), rho = opts.rho; end
if isfield(opts,'eta'), eta = opts.eta; end
if isfield(opts,'phi'), phi = opts.phi; end

% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial Tau & Eta
tau = tau * ones(dim,dim);
eta = eta * ones(dim,dim);
% Pre
fitG = inf;
fit  = zeros(1,N);
tau0 = tau;

curve = inf;
t = 1;
% Iterations
while t <= max_Iter
	% Reset ant
	X = zeros(N,dim);
	for i=1:N
    % Set number of features
    num_feat = randi([1,dim]);
    % Ant start with random position
    X(i,1)   = randi([1,dim]);
    k        = [];
    if num_feat > 1
      for d = 2:num_feat
        % Start with previous tour
        k      = [k(1:end), X(i, d-1)];
        % Edge / Probability Selection (4)
        P      = (tau(k(end),:) .^ alpha) .* (eta(k(end),:) .^ beta);
        % Set selected position = 0 probability (4)
        P(k)   = 0;
        % Convert probability (4)
        prob   = P ./ sum(P(:));
        % Roulette Wheel selection
        route  = jRouletteWheelSelection(prob);
        % Store selected position to be next tour
        X(i,d) = route;
      end
    end
  end
  % Binary
  X_bin = zeros(N,dim);
  for i = 1:N
    % Binary form
    ind           = X(i,:);
    ind(ind == 0) = [];
    X_bin(i, ind) = 1;
  end
  % Binary version
  for i = 1:N
    % Fitness
    fit(i) = fun(feat,label,X_bin(i,:),opts);
    % Global update
    if fit(i) < fitG
      Xgb  = X(i,:);
      fitG = fit(i);
    end
  end
  % Tau update
  tour            = Xgb;
  tour(tour == 0) = [];
  tour            = [tour(1:end), tour(1)];
  for d = 1 : length(tour) - 1
    % Feature selected
    x = tour(d);
    y = tour(d + 1);
    % Delta tau
    Dtau = 1 / fitG;
    % Update tau (10)
    tau(x,y) = (1 - phi) * tau(x,y) + phi * Dtau;
  end
  % Evaporate pheromone (9)
  tau = (1 - rho) * tau + rho * tau0;
  % Save
  curve(t) = fitG;
  fprintf('\nIteration %d Best (ACS)= %f',t,curve(t))
  t = t + 1;
end
% Select features based on selected index
Sf = unique(Xgb);
Sf(Sf == 0) = [];
sFeat = feat(:,Sf);
% Store results
ACS.sf = Sf;
ACS.ff = sFeat;
ACS.nf = length(Sf);
ACS.c  = curve;
ACS.f  = feat;
ACS.l  = label;
end

%// Roulette Wheel Selection //
function Index = jRouletteWheelSelection(prob)
% Cummulative summation
C = cumsum(prob);
% Random one value, most probability value [0~1]
P = rand();
% Route wheel
for i = 1:length(C)
	if C(i) > P
    Index = i;
    break;
  end
end
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
% op_freqshift.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% out=op_freqshift(in,f);
%
% DESCRIPTION:
% Apply a frequency shift to the input spectrum by 'f' Hz.
%
% INPUTS:
% in     = input data in matlab structure format.
% f      = frequency shift to apply (in Hz).
%
% OUTPUTS:
% out    = Output following frequency shift.

function out=op_freqshift(in,f);

% if in.dims.coils>0
%     error('ERROR:  Can not operate on data with multilple coils!  ABORTING!!')
% end
% if in.dims.averages>0
%     error('ERROR:  Can not operate on data with multiple averages!  ABORTING!!');
% end
% if in.dims.subSpecs>0
%     error('ERROR:  Can not operate on data with multiple Subspecs!  ABORTING!!');
% end

t=repmat(in.t',[1 in.sz(2:end)]);

fids=in.fids.*exp(-1i*t*f*2*pi);

%re-calculate Specs using fft
specs=fftshift(ifft(fids,[],in.dims.t),in.dims.t);

%plot(in1.ppm,combinedSpecs);

%FILLING IN DATA STRUCTURES
out=in;
out.fids=fids;
out.specs=specs;

%FILLING IN THE FLAGS
out.flags=in.flags;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [PB] = EB2PB(EB)
% Convert computery things from exabytes to petabytes.
% Chad A. Greene 2012
PB = EB*1024 ;


text:
function [params, names] = vargplvmExtractParamNoVardist(model)

% VARGPLVMEXTRACTPARAMNOVARDIST Extract a parameter vector from a variational GP-LVM model.
% ignoring the variational distribution. See vargplvmExtractParam.

if nargout > 1
  returnNames = true;
else
  returnNames = false;
end
params = [];
names = {};

% Inducing inputs
if ~model.fixInducing
    if ~isfield(model, 'learnInducing') || (isfield(model, 'learnInducing') && model.learnInducing)
        params =  [params model.X_u(:)'];
        if returnNames
            for i = 1:size(model.X_u, 1)
                for j = 1:size(model.X_u, 2)
                    X_uNames{i, j} = ['X_u(' num2str(i) ', ' num2str(j) ')'];
                end
            end
            names = {names{:}, X_uNames{:}};
        end
    end
end

% Kernel parameters
if returnNames
  [kernParams, kernParamNames] = kernExtractParam(model.kern);
  for i = 1:length(kernParamNames)
    kernParamNames{i} = ['Kernel, ' kernParamNames{i}];
  end
  names = {names{:}, kernParamNames{:}};
else
  kernParams = kernExtractParam(model.kern);
end
params = [params kernParams];

% beta in the likelihood
if model.optimiseBeta

   if ~isstruct(model.betaTransform)
       fhandle = str2func([model.betaTransform 'Transform']);
       betaParam = fhandle(model.beta, 'xtoa');
   else
      if isfield(model.betaTransform,'transformsettings') && ~isempty(model.betaTransform.transformsettings)
          fhandle = str2func([model.betaTransform.type 'Transform']);
          betaParam = fhandle(model.beta, 'xtoa', model.betaTransform.transformsettings);
      else
          error('vargplvmExtractParam: Invalid transform specified for beta.');
      end
   end

   params = [params betaParam(:)'];

   if returnNames
     for i = 1:length(betaParam)
       betaParamNames{i} = ['Beta ' num2str(i)];
     end
     names = {names{:}, betaParamNames{:}};
   end
end

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
function desc = calHsi(rgb_im, seg, numRegion)

    hsi = rgb2hsi_refine(rgb_im);

    ind={};
    for iReg=1:numRegion
        ind{iReg} = seg(:)==iReg;
    end

    desc = zeros([numRegion 3]);

    for ch=1:3
        for iReg = 1:numRegion
            feature = hsi(:,:,ch);
            feature = feature(ind{iReg});
            desc(iReg, ch) = sum(feature) / sum(ind{iReg});
        end
    end

end

author: kittenish

repo: Image-Shadow-Detection-and-Removal


text:
% OP_GRADSYMV_N_U: assemble the matrix A = [a(i,j)], a(i,j) = (epsilon (gradsym v n)_j, u_i), with n the normal vector.
%
%   mat = op_gradsymv_n_u (spu, spv, msh, epsilon);
%   [rows, cols, values] = op_gradsymv_n_u (spu, spv, msh, epsilon);
%
% INPUT:
%
%   spu:   structure representing the space of trial functions (see sp_scalar/sp_evaluate_col)
%   spv:   structure representing the space of test functions (see sp_scalar/sp_evaluate_col)
%   msh:   structure containing the domain partition and the quadrature rule for the boundary,
%           since it must contain the normal vector (see msh_cartesian/msh_eval_boundary_side)
%   epsilon: coefficient
%
% OUTPUT:
%
%   mat:    assembled matrix
%   rows:   row indices of the nonzero entries
%   cols:   column indices of the nonzero entries
%   values: values of the nonzero entries
%
% Copyright (C) 2014 Adriano Cortes
% Copyright (C) 2014, 2017, 2020 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function varargout = op_gradsymv_n_u (spu, spv, msh, coeff)

  gradv = reshape (spv.shape_function_gradients, spv.ncomp, [], ...
		   msh.nqn, spv.nsh_max, msh.nel);

  ndim = size (gradv, 2);

  shpu = reshape (spu.shape_functions, spu.ncomp, msh.nqn, spu.nsh_max, msh.nel);

  rows = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  cols = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  values = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);

  jacdet_weights = msh.jacdet .* msh.quad_weights .* coeff;

  ncounter = 0;
  for iel = 1:msh.nel
    if (all (msh.jacdet(:,iel)))
      gradv_iel = gradv(:,:,:,:,iel);
      normal_iel = reshape (msh.normal(:,:,iel), [1, ndim, msh.nqn]);
      %Symmetrize gradv
      gradv_iel = 0.5*(gradv_iel + permute (gradv_iel, [2 1 3 4]));

      gradv_n = reshape (sum (bsxfun (@times, gradv_iel, normal_iel), 2), spv.ncomp, msh.nqn, spv.nsh_max, 1);
      shpu_iel = reshape (shpu(:, :, :, iel), spu.ncomp, msh.nqn, 1, spu.nsh_max);

      jacdet_iel = reshape (jacdet_weights(:,iel), [1,msh.nqn,1,1]);

      gradv_n_times_jw = bsxfun (@times, jacdet_iel, gradv_n);
      tmp1 = sum (bsxfun (@times, gradv_n_times_jw, shpu_iel), 1);
      elementary_values = reshape (sum (tmp1, 2), spv.nsh_max, spu.nsh_max);

      [rows_loc, cols_loc] = ndgrid (spv.connectivity(:,iel), spu.connectivity(:,iel));
      indices = rows_loc & cols_loc;
      rows(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = rows_loc(indices);
      cols(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = cols_loc(indices);
      values(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = elementary_values(indices);
      ncounter = ncounter + spu.nsh(iel)*spv.nsh(iel);

    else
      warning ('geopdes:jacdet_zero_at_quad_node', 'op_gradv_n_u: singular map in element number %d', iel)
    end
  end

  if (nargout == 1 || nargout == 0)
    varargout{1} = sparse (rows, cols, values, spv.ndof, spu.ndof);
  elseif (nargout == 3)
    varargout{1} = rows;
    varargout{2} = cols;
    varargout{3} = values;
  else
    error ('op_gradv_n_u: wrong number of output arguments')
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [m3] = cl2m3(cl)
% Convert volume from centiliters to cubic meters.
% Chad Greene 2012
m3 = cl*0.00001;


text:
function [ w, x ] = line_o02 ( )

%*****************************************************************************80
%
%% LINE_O02 returns a 2 point quadrature rule for the unit line.
%
%  Discussion:
%
%    The integration region is:
%
%    - 1.0 <= X <= 1.0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carlos Felippa,
%    A compendium of FEM integration formulas for symbolic work,
%    Engineering Computation,
%    Volume 21, Number 8, 2004, pages 867-890.
%
%  Parameters:
%
%    Output, real W(2), the weights.
%
%    Output, real X(2), the abscissas.
%
  w(1:2) = [ ...
    0.5, ...
    0.5 ];

  x(1:2) = [ ...
    -0.57735026918962576451, ...
     0.57735026918962576451 ];

  return
end


text:
function [d pred f]=astar_search(A,s,h,varargin)
% ASTAR_SEARCH Perform a heuristically guided (A*) search on the graph.
%
% [d pred rank]=astar_search(A,s,h,optionsu) returns the distance map,
% search tree and f-value of each node in an astar_search.
% The search begins at vertex s.  The heuristic h guides the search,
% h(v) should be small close to a goal and large far from a goal.  The
% heuristic h can either be a vector with an entry for each vertex in the
% graph or a function which maps vertices to values.
%
% This method works on non-negatively weighted directed graphs.
% The runtime is O((E+V)log(V)).
%
% ... = astar_search(A,u,...) takes a set of
% key-value pairs or an options structure.  See set_matlab_bgl_options
% for the standard options.
%   options.visitor: a visitor to use with the A* search (see Note)
%   options.inf: the value to use for unreachable vertices
%       [double > 0 | {Inf}]
%   options.target: a special vertex that will stop the search when hit
%       [{'none'} | any vertex number besides the u]; this is ignored if
%       a visitor is set.
%   options.edge_weight: a double array over the edges with an edge
%       weight for each node, see EDGE_INDEX and EXAMPLES/REWEIGHTED_GRAPHS
%       for information on how to use this option correctly
%       [{'matrix'} | length(nnz(A)) double vector]
%
% Note: You can specify a visitor for this algorithm.  The visitor has the
% following optional functions.
%    vis.initialize_vertex(u)
%    vis.discover_vertex(u)
%    vis.examine_vertex(u)
%    vis.examine_edge(ei,u,v)
%    vis.edge_relaxed(ei,u,v)
%    vis.edge_not_relaxed(ei,u,v)
%    vis.black_target(ei,u,v)
%    vis.finish_vertex(u)
% Each visitor parameter should be a function pointer, which returns 0
% if the search should stop.  (If the function does not return anything,
% the algorithm continues.)
%
% Example:
%   load graphs/bgl_cities.mat
%   goal = 11; % Binghamton
%   start = 9; % Buffalo
%   % Use the euclidean distance to the goal as the heuristic
%   h = @(u) norm(xy(u,:) - xy(goal,:));
%   % Setup a routine to stop when we find the goal
%   ev = @(u) (u ~= goal);
%   [d pred f] = astar_search(A, start, h, ...
%       struct('visitor', struct('examine_vertex', ev)));

% David Gleich
% Copyright, Stanford University, 2006-2008

%% History
%  2007-04-20: Added edge weight option
%  2007-07-12: Fixed edge_weight documentation.
%  2008-10-07: Changed options parsing
%%

[trans check full2sparse] = get_matlab_bgl_options(varargin{:});
if full2sparse && ~issparse(A), A = sparse(A); end

options = struct('inf', Inf, 'edge_weight', 'matrix', 'target', 'none');
options = merge_options(options,varargin{:});

edge_weight_opt = 'matrix';

if strcmp(options.edge_weight, 'matrix')
    % do nothing if we are using the matrix weights
else
    edge_weight_opt = options.edge_weight;
end

if strcmp(options.target,'none')
    target = 0; % a flag used to denote "no target" to the mex
elseif isa(options.target, 'double')
    target = options.target;
else
    error('matlab_bgl:invalidParameter', ...
        'options.target is not ''none'' or a vertex number.');
end

if check, check_matlab_bgl(A,struct()); end
if trans, A = A'; end

function hi=vec2func(u)
    hi = h(u);
end

if isa(h,'function_handle')
    hfunc = h;
else
    hfunc = @vec2func;
end

if isfield(options,'visitor')
    [d pred f] = astar_search_mex(A,s,target,hfunc,options.inf,edge_weight_opt,options.visitor);
else
    [d pred f] = astar_search_mex(A,s,target,h,options.inf,edge_weight_opt);
end

% end the main function
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


%==============================================================================
% This code is part of the Finite Element Method app for the Matlab-based toolbox
%  FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR/FAIRFEM
%==============================================================================
%
% classdef TriMesh2 < handle
%
% Finite Element Mesh based on triangular subdivision of rectangular mesh.
%
% Each Cell is divided into two triangles:
%
%  o---------o---------o---------o---------o
%  | \       | \       | \       | \       |
%  |   \     |   \     |   \     |   \     |
%  |     \   |     \   |     \   |     \   |
%  |       \ |       \ |       \ |       \ |
%  o---------o---------o---------o---------o
%  | \       | \       | \       | \       |
%  |   \     |   \     |   \     |   \     |
%  |     \   |     \   |     \   |     \   |
%  |       \ |       \ |       \ |       \ |
%  o---------o---------o---------o---------o
%  | \       | \       | \       | \       |
%  |   \     |   \     |   \     |   \     |
%  |     \   |     \   |     \   |     \   |
%  |       \ |       \ |       \ |       \ |
%  o---------o---------o---------o---------o
%
%  To construct an instance of this class type:
%
%  >> Mesh = TriMesh2(omega,m)
%
% Input:
% 	omega - description of spatial domain
%   m     - number of cells
%
% Properties:
%   xn     - node list
%   tri    - triangle list
%   dim    - space dimension
%   omega  - description of spatial domain
%   m      - number of cells
%   type   - type of partition
%   vol    - volume of triangles
%   nnodes - number of nodes
%   ntri   - number of triangles
%   dx1    - partial derivative operator
%   dx2    - partial derivative operator
%   GRAD   - gradient operator
%   P1     - projection operator for node 1
%   P2     - projection operator for node 2
%   P3     - projection operator for node 3
%   PC     - projection operator for Barycentrum
%   P      - prolongation operator
%   Pt     - restriction operator
%
%  Methods:
%   mfPu   - matrix free prolongation/restriction
%   mfPi   - matrix free edge projector
%   getP   - builds prolongation operator
%   tri2cc - averaging
%
%
% see also
% =========================================================================
classdef TriMesh2 < handle

    properties
        % ===================================
        % node list
        % ===================================
        xn
        % ===================================
        % triangle list
        % ===================================
        tri
        % ===================================
        % space dimension
        % ===================================
        dim  = 2;
        % ===================================
        % description of computational domain
        % ===================================
        omega
        % ===================================
        % number of cells
        % ===================================
        m
        % ===================================
        % type of partition
        % ===================================
        type = 1;
        % ===================================
        % function handle to myself
        % ===================================
        me = @TriMesh2
        % ===================================
        % volume of triangles
        % ===================================
        vol
        % ===================================
        % number of nodes
        % ===================================
        nnodes
        % ===================================
        % number of triangles in mesh
        % ===================================
        ntri
    end

    properties (Access = public, Dependent) % These will be created when first callend and stores persistently
        % ===================================
        % dx1  - partial derivative operator
        % ===================================
        dx1
        % ===================================
        % dx2  - Partial derivative operator
        % ===================================
        dx2
        % ===================================
        % GRAD - Gradient operator
        %
        %         | dx1 |
        %  GRAD = |     |
        %         | dx2 |
        %
        % ===================================
        GRAD
        % ===================================
        % B - Vector gradient operator
        %
        %         | GRAD  0    |
        %  B =    |            |
        %         | 0     GRAD |
        %
        % ===================================
        B
        % ===================================
        % P1 - Projection operator on Node 1
        % ===================================
        P1
        % ===================================
        % P2 - Projection operator on Node 2
        % ===================================
        P2
        % ===================================
        % P3 - Projection operator on Node 3
        % ===================================
        P3
        % ===================================
        % PC - Projection operator on Barycenter
        % ===================================
        PC
        % ===================================
        % P  - Prolongation operator
        % ===================================
        P
        % ===================================
        % Pt - Restriction operator
        % ===================================
        Pt
        % ===================================
        % Boundary indices
        % ===================================
        boundaryIdx
        % ===================================
        % Boundary projector
        % ===================================
        boundaryProj
        mfdx1
        mfdx2
        mfGRAD
    end

    properties (Access = private)
        % These are where the dependent data is actually stored
        dx1_
        dx2_
        GRAD_
        B_
        P1_
        P2_
        P3_
        PC_
        P_
        Pt_
        mfdx1_
        mfdx2_
        mfGRAD_
        boundaryIdx_
        boundaryProj_
    end

    methods
        function this = TriMesh2(omega,m)
            if nargin==0,
                help(mfilename);
                this.runMinimalExample;
                return;
            end
            this.omega = omega;
            this.m     = m;

            this.xn = reshape(getNodalGrid(omega,m),[],2);
            % get indices of bottom left vertices
            nodes = reshape(1:prod(m+1),m+1);
            nodes = reshape(nodes(1:end-1,1:end-1),1,[]);
            % specify triangles
            this.tri   = [nodes; nodes+1; nodes+m(1)+1; nodes+m(1)+2; nodes+m(1)+1; nodes+1];
            this.tri   = reshape(this.tri,3,[])';
            this.nnodes = prod(m+1);
            this.ntri  = size(this.tri,1);
            this.vol = prod((omega(2:2:end)-omega(1:2:end))./m)/2*ones(this.ntri,1);
        end

        function runMinimalExample(~)
            omega = [0 4 2 6]; m = [8 16];
            Mesh  = feval(mfilename,omega,m);
        end

        function x = mfPi(this,x,i)
            % =============================================================
            % function x = mfPi(this,x,i)
            %
            % matrix free edge projector
            % =============================================================
            switch i
                case 1
                    P = this.P1;
                case 2
                    P = this.P2;
                case 3
                    P = this.P3;
                case 'C'
                    P = this.PC;
            end
            if size(x,1) == this.ntri,
                % ajoint
                x = P'*x;
            else
                x = P * x;
            end
        end

        function x = tri2cc(this,x)
            % =============================================================
            % function x = tri2cc(x)
            %
            % averaging or adjoint
            % =============================================================
            if numel(x)==this.ntri,
                x = mean(reshape(x,2,[]),1);
            else
                x = reshape(x,1,[]);
                x = .5*[x;x];
                x = x(:);
            end

        end

        function Pu = mfPuNodal(~,yn,m,flag)
            % =============================================================
            % function Pu = mfPuNodal(~,yn,m,flag)
            %
            % matrix free prolongation/restriction for nodal quantities
            % =============================================================
            d  = numel(yn)/prod(m+1);
            yn = reshape(yn,[m+1, d]);
            switch flag
                case 'Pu' % coarse --> fine
                    yn = reshape(yn,[m+1 2]);
                    Pu = zeros([2*m+1 2]);

                    % include existing nodes
                    Pu(1:2:end,1:2:end,:) = yn;
                    % prolongate in x direction
                    Pu(2:2:end-1,:,:) =  .5*(Pu(1:2:end-2,:,:) + Pu(3:2:end,:,:));
                    % prolongate in y direction
                    Pu(1:2:end,2:2:end-1,:) =  .5*(Pu(1:2:end,1:2:end-2,:) + Pu(1:2:end,3:2:end,:));
                    % prolongate diagonal
                    Pu(2:2:end,2:2:end,:) = .5* (Pu(3:2:end,1:2:end-2,:) + Pu(1:2:end-2,3:2:end,:));

                    Pu = reshape(Pu,[],1);
                case 'PTu'
                    yn = reshape(yn,[m+1 2]);

                    % include parent nodes
                    Pu = yn(1:2:end,1:2:end,:);

                    % distribute in x direction
                    t = .5* yn(2:2:end-1,1:2:end,:);
                    Pu(1:end-1,:,:) = Pu(1:end-1,:,:) + t;
                    Pu(2:end,:,:)   = Pu(2:end,:,:) + t;

                    % distribute in y direction
                    t = .5* yn(1:2:end,2:2:end-1,:);
                    Pu(:,1:end-1,:) = Pu(:,1:end-1,:) + t;
                    Pu(:,2:end,:) = Pu(:,2:end,:) + t;

                    % distribute diagonal
                    t = .5* yn(2:2:end-1,2:2:end-1,:);
                    Pu(2:end,1:end-1,:) = Pu(2:end,1:end-1,:) + t;
                    Pu(1:end-1,2:end,:) = Pu(1:end-1,2:end,:) + t;

                    Pu = reshape(Pu,[],1);
            end
        end

        function P = getPuNodal(~,m)
            % =============================================================
            % function P = getPuNodal(~,m)
            %
            % returns prolongation operator for input of cell-width m
            % =============================================================

            mf = 2*m;
            % indices of fine and coarse grid nodes
            indf = reshape(1:prod(mf+1),mf+1);
            indc = reshape(1:prod(m+1),m+1);

            % allocate space
            I = []; J = []; W = [];

            % include existing nodes
            ii = indf(1:2:end,1:2:end);
            jj = indc(:);
            ww = ones(size(jj));
            I = [I; ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];

            % prolongate in x direction
            ii   = indf(2:2:end-1,1:2:end);
            jj   = [reshape(indc(1:end-1,:),[],1); reshape(indc(2:end,:),[],1)];
            ww = .5*ones(size(jj));
            I = [I; ii(:);ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];

            % prolongate in y direction
            ii   = indf(1:2:end,2:2:end-1);
            jj   = [reshape(indc(:,1:end-1),[],1); reshape(indc(:,2:end),[],1)];
            ww   = .5*ones(size(jj));
            I = [I; ii(:); ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];

            % prolongate diagonally
            ii   = indf(2:2:end-1,2:2:end-1);
            jj   = [reshape(indc(1:end-1,2:end),[],1); reshape(indc(2:end,1:end-1),[],1)];
            ww   = .5*ones(size(jj));
            I = [I; ii(:);ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];

            P = sparse(I,J,W,prod(mf+1),prod(m+1));

        end

        % ========== get methods ========================================
        function dx1 = get.dx1(this)
            if isempty(this.dx1_),
                [this.dx1_, this.dx2_] = getGradientMatrixFEM(this,0);
            end
            dx1 = this.dx1_;
        end

        function dx2 = get.dx2(this)
            if isempty(this.dx2_),
                [this.dx1_, this.dx2_] = getGradientMatrixFEM(this,0);
            end
            dx2 = this.dx2_;
        end

        function GRAD = get.GRAD(this)
            if isempty(this.GRAD_),
                this.GRAD_ = [this.dx1;this.dx2];
            end
            GRAD = this.GRAD_;
        end
        function B = get.B(this)
            if isempty(this.B_),
                this.B_ = blkdiag(this.GRAD,this.GRAD);
            end
            B = this.B_;
        end

        function P1 = get.P1(this)
            if isempty(this.P1_),
                A = speye(this.nnodes);
                this.P1_ = A(this.tri(:,1),:);
            end
            P1 = this.P1_;
        end
        function P2 = get.P2(this)
            if isempty(this.P2_),
                A = speye(this.nnodes);
                this.P2_ = A(this.tri(:,2),:);
            end
            P2 = this.P2_;
        end

        function P3 = get.P3(this)
            if isempty(this.P3_),
                A = speye(this.nnodes);
                this.P3_ = A(this.tri(:,3),:);
            end
            P3 = this.P3_;
        end
        function PC = get.PC(this)
            if isempty(this.PC_),
                this.PC_ = (this.P1+this.P2+this.P3)/3;
            end
            PC = this.PC_;
        end
        function P = get.P(this)
            if isempty(this.P_),
                this.P_ = this.getPuNodal(this.m);
            end
            P = this.P_;
        end
        function Pt = get.Pt(this)
            if isempty(this.Pt_),
                this.Pt_ = this.getPuNodal(this.m/2);
            end
            Pt = this.Pt_;
        end

        function mfdx1 = get.mfdx1(this)
            if isempty(this.mfdx1_),
                [this.mfdx1_, this.mfdx2_] = getGradientMatrixFEM(this,1);
            end
            mfdx1 = this.mfdx1_;
        end

        function mfdx2 = get.mfdx2(this)
            if isempty(this.mfdx2_),
                [this.mfdx1_, this.mfdx2_] = getGradientMatrixFEM(this,1);
            end
            mfdx2 = this.mfdx2_;
        end

        function mfGRAD = get.mfGRAD(this)
            if isempty(this.mfGRAD_),
                this.mfGRAD_ = getGradientMatrixFEM(this,1);
            end
            mfGRAD = this.mfGRAD_;
        end

        function idx = get.boundaryIdx(this)
            if isempty(this.boundaryIdx_),
                id = reshape(1:prod(this.m+1),this.m+1);
                idx = [reshape(id([1,end],:),[],1);reshape(id(2:end-1,[1,end]),[],1)] ;
                this.boundaryIdx_ = idx;
            end
            idx = this.boundaryIdx_;
        end

        function idx = get.boundaryProj(this)
            if isempty(this.boundaryProj_),
                idx = this.boundaryIdx;

                P = speye(prod(this.m+1));
                P = P(idx,:);
                P = kron(speye(2),P);

                this.boundaryProj_ = P;
            end
            idx = this.boundaryProj_;
        end

        % ========== set methods ========================================
        function set.xn(this,xn)
            if numel(xn)~=(2*prod(this.m+1)),
                error('Invalid number of nodes');
            end
            if isempty(this.xn),
                this.xn = xn;
            else

                this.xn = xn;
                % delete operators that are sensitive to nodes
                this.dx1_  = [];
                this.dx2_  = [];
                this.GRAD_ = [];
                this.B_    = [];

                this.vol = volTetraGrid(this,this.xn,'matrixFree',1);
            end
        end
    end
end


text:
function meandpth() % autogenerated function wrapper
    %  meandpth  finds the average depth for a predefined running
    %  window in terms of number of events, and a selected step
    %  and plots the results.
    %                                                     R.Z. 6/94
    %                          Operates on ZG.newcat
    %
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    report_this_filefun(mfilename('fullpath'));

    if ic == 1 || ic == 0
        if isempty(ZG.newcat) ZG.newcat = ZG.primeCatalog; end
        ZG.newcat = ZG.primeCatalog;
        iwln = 100;
        step = 28;
        winlen_days = 2;

        figure_w_normalized_uicontrolunits(...
            'Name','MeanDepth Input Parameters',...
            'visible','off',...
            'NumberTitle','off', ...
            'NextPlot','new', ...
            'Units','Pixel',  'Position',[ZG.welcome_pos 550 200'])
        axis off
        set(gca,'visible','off');

        % creates a dialog box to input some parameters
        %

        freq_field=uicontrol('Style','edit',...
            'Position',[.70 .60 .17 .10],...
            'Units','normalized','String',num2str(iwln),...
            'callback',@callbackfun_001);

        inp2_field=uicontrol('Style','edit',...
            'Position',[.70 .40 .17 .10],...
            'Units','normalized','String',num2str(step),...
            'callback',@callbackfun_002);

        close_button=uicontrol('Style','Pushbutton',...
            'Position', [.60 .05 .15 .15 ],...
            'Units','normalized','callback',@callbackfun_003,'String','Cancel');

        go_button=uicontrol('Style','Pushbutton',...
            'Position',[.25 .05 .15 .15 ],...
            'Units','normalized',...
            'callback',@callbackfun_004,...
            'String','Go');

        txt1 = text(...
            'Position',[0. 0.65 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold' ,...
            'String','Number of events in averaging window:');

        txt2 = text(...
            'Position',[0. 0.40 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold' ,...
            'String','Time-step in days:');

        set(gcf,'visible','on')

    elseif ic == 2

        step = days(step);
        if isempty(ZG.newcat) ZG.newcat = ZG.primeCatalog; end
        len = ZG.newcat.Count;
        xt2  = [ ];
        meand = [ ];
        er = [];
        [t0b, teb] = ZG.newcat.DateRange() ;
        ind = 0;

        me = [];
        b = ZG.primeCatalog;

        sta = ZG.primeCatalog.Date(iwln+1);
        sto = max(ZG.primeCatalog.Date);
        for it=sta+step:step:sto
            ind = ind + 1;
            t = find( abs(b.Date-it)  == min(abs(b.Date-it)));
            if t <= iwln ; t = iwln+1; end
            meand(ind) =  mean(ZG.newcat(t-iwln:t,7)) ;
            er(ind) =  std(ZG.newcat(t-iwln:t,7)) ;
            xt2(ind) = it;        % time is end of window
        end    % for it

        meand = -meand;

        % Find out if figure already exists
        %
        depfg=findobj('Type','Figure','-and','Name','Mean Depth');

        % Set up the Seismicity Map window Enviroment
        %
        if isempty(depfg)

            depfg=figure_w_normalized_uicontrolunits(...
                'Name','Mean Depth',...
                'visible','off',...
                'NumberTitle','off', ...
                'NextPlot','new', ...
                'Units','Pixel',  'Position',[ZG.welcome_pos 550 400']);
            hold on
            axis off

            uicontrol('Style','Pushbutton',...
                'Position',[.9 .80 .10 .05],...
                'Units','normalized',...
                'callback',@(~,~)medispas1('ast'),'String','AS');
            uicontrol('Style','Pushbutton',...
                'Position',[.9 .70 .10 .05],...
                'Units','normalized',...
                'callback',@(~,~)medispas1('lta'),'String','LTA');
            uicontrol('Style','Pushbutton',...
                'Position',[.9 .90 .10 .05],...
                'Units','normalized',...
                'callback',@callbackfun_007,'String','Com');

            new = uicontrol('style','edit','value',winlen_days,...
                'string',num2str(winlen_days), 'background','y',...
                'callback',@callbackfun_008,...
                'units','norm','pos',[.90 .30 .10 .06],'min',0.1,'max',100);

            newlabel = uicontrol('style','text','units','norm','pos',[.85 .30 .05 .06]);
            set(newlabel,'string','winlen_days:','background',color_fbg);

            uicontrol('Units','normal',...
                'Position',[.90 .25 .08 .06],'String','Go',...
                'callback',@callbackfun_009)

        end  % if figure exist

        figure(depfg);
        delete(findobj(depfg,'Type','axes'));
        set(gca,'visible','off');

        %orient tall
        set(gcf,'Units','centimeter','PaperPosition',[1 1 7 8])
        rect = [0.15, 0.15, 0.65, 0.30];
        axes('position',rect)
        p5 = gca;

        % plot errbar
        errorbar(xt2,meand,er)
        hold on
        plot(xt2,meand,'co')
        pl = plot(xt2,meand,'-r')
        hold on
        set(pl,'LineWidth',1.0)
        if ~isempty(ZG.maepi)
            pl =   plot(ZG.maepi.Date,ZG.maepi.Date*0+mean(meand),'xm');
            set(pl,'LineWidth',2.0)
        end

        axis([min(ZG.newcat.Date) max(ZG.newcat.Date+0.1) min(meand*1.1)  max(meand*0.9)])
        v = axis;
        xlabel('Time (years)','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m,'Color','k')
        ylabel('Mean Depth (km)','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m,'Color','k')
        stri = ['Mean Depths and standard deviation ( ' file1 ')'];
        %title(' Mean depths and mean depth error ',...
        %        'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m,'Color','k')
        set(gca,'box','on',...
            'SortMethod','childorder','TickDir','out','FontWeight',...
            'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)

        grid
        hold off

        rect = [0.15,  0.45, 0.65, 0.30];
        axes('position',rect)
        pl =plot(ZG.newcat.Date,-ZG.newcat.Depth,'.b')
        set(pl,'MarkerSize',3')
        set(pl,'LineWidth',1.0)
        hold on
        if ~isempty(ZG.maepi)
            pl =  plot(ZG.maepi.Date,-ZG.maepi.Depth,'xm');
            set(pl,'LineWidth',2.0)
        end
        axis([ v(1) v(2) -max(ZG.newcat.Depth)  -min(ZG.newcat.Depth)])
        ylabel('Depth (km)','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m,'Color','k')
        set(gca,'XTicklabels',[])
        stro = [' ' file1 '; wl = ' num2str(iwln) ' events, inc = ' num2str(step)];
        title(stro,'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m,'Color','k')
        grid

        set(gca,'box','on',...
            'SortMethod','childorder','TickDir','out','FontWeight',...
            'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)

        set(gca,'visible','on');
        set(gcf,'visible','on');

        ic = 1;

    end    % if ic

    function callbackfun_001(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        iwln=str2double(freq_field.String);
        freq_field.String=num2str(iwln);
    end

    function callbackfun_002(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        step=str2double(inp2_field.String);
        inp2_field.String=num2str(step);
    end

    function callbackfun_003(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        close;

    end

    function callbackfun_004(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ic = 2;
        close;
        clear meandpth;
        meandpth;
    end

    function callbackfun_007(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        dispma4();
    end

    function callbackfun_008(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        update_editfield_value(mysrc,myevt);
        winlen_days=mysrc.Value;
        medispas1();
    end

    function callbackfun_009(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        medispas1();
    end

end

author: CelsoReyes