text:
function x = strmv ( uplo, trans, diag, n, a, lda, x, incx )

%*****************************************************************************80
%
%% STRMV computes x: = A*x or x = A'*x for a triangular matrix A.
%
%  Discussion:
%
%    STRMV performs one of the matrix-vector operations
%
%      x := A*x,   or   x := A'*x,
%
%    where x is an n element vector and  A is an n by n unit, or non-unit,
%    upper or lower triangular matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 April 2014
%
%  Author:
%
%    This MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, character UPLO, specifies whether the matrix is an upper or
%    lower triangular matrix as follows:
%    'u' or 'U': A is an upper triangular matrix.
%    'l' or 'L': A is a lower triangular matrix.
%
%    Input, character TRANS, specifies the operation to be performed as
%    follows:
%    'n' or 'N': x := A*x.
%    't' or 'T': x := A'*x.
%    'c' or 'C': x := A'*x.
%
%    Input, character DIAG, specifies whether or not A is unit
%    triangular as follows:
%    'u' or 'U': A is assumed to be unit triangular.
%    'n' or 'N': A is not assumed to be unit triangular.
%
%    Input, integer N, the order of the matrix A.
%    0 <= N.
%
%    Input, real A(LDA,N).
%    Before entry with  UPLO = 'u' or 'U', the leading n by n
%    upper triangular part of the array A must contain the upper
%    triangular matrix and the strictly lower triangular part of
%    A is not referenced.
%    Before entry with UPLO = 'l' or 'L', the leading n by n
%    lower triangular part of the array A must contain the lower
%    triangular matrix and the strictly upper triangular part of
%    A is not referenced.
%    Note that when  DIAG = 'u' or 'U', the diagonal elements of
%    A are not referenced either, but are assumed to be unity.
%
%    Input, integer LDA, the first dimension of A as declared
%    in the calling program. max ( 1, N ) <= LDA.
%
%    Input/output, real X(1+(N-1)*abs( INCX)).
%    Before entry, the incremented array X must contain the n
%    element vector x. On exit, X is overwritten with the
%    tranformed vector x.
%
%    Input, integer INCX, the increment for the elements of
%    X.  INCX must not be zero.
%

%
%  Test the input parameters.
%
  info = 0;
  if  ( ~ lsame ( uplo , 'U' ) && ~ lsame ( uplo , 'L' ) )
    info = 1;
  elseif ( ~ lsame ( trans, 'N' ) && ~ lsame ( trans, 'T' ) && ...
            ~ lsame ( trans, 'C' ) )
    info = 2;
  elseif ( ~ lsame ( diag , 'U' ) && ~ lsame ( diag , 'N' ) )
    info = 3;
  elseif ( n < 0 )
    info = 4;
  elseif ( lda < max ( 1, n ) )
    info = 6;
  elseif ( incx == 0 )
    info = 8;
  end

  if ( info ~= 0 )
    xerbla ( 'DTRMV', info );
    return
  end
%
%  Quick return if possible.
%
  if ( n == 0 )
    return
  end

  nounit = lsame ( diag, 'N' );
%
%  Set up the start point in X if the increment is not unity. This
%  will be  ( N - 1 ) * INCX  too small for descending loops.
%
  if ( incx <= 0 )
    kx = 1 - ( n - 1 ) * incx;
  elseif ( incx ~= 1 )
    kx = 1;
  end
%
%  Start the operations. In this version the elements of A are
%  accessed sequentially with one pass through A.
%
  if ( lsame ( trans, 'N' ) )
%
%  Form x := A*x.
%
    if ( lsame ( uplo, 'U' ) )
      if ( incx == 1 )
        for j = 1 : n
          if ( x(j) ~= 0.0 )
            temp = x(j);
            for i = 1 : j - 1
              x(i) = x(i) + temp * a(i,j);
            end
            if ( nounit )
              x(j) = x(j) * a(j,j);
            end
          end
        end
      else
        jx = kx;
        for j = 1 : n
          if ( x(jx) ~= 0.0 )
            temp = x(jx);
            ix = kx;
            for i = 1 : j - 1
              x(ix) = x(ix) + temp * a(i,j);
              ix = ix + incx;
            end
            if ( nounit )
              x(jx) = x(jx) * a(j,j);
            end
          end
          jx = jx + incx;
        end
      end
    else
      if ( incx == 1 )
        for j = n : -1 : 1
          if ( x(j) ~= 0.0 )
            temp = x(j);
            for i = n : -1 : j + 1
              x(i) = x(i) + temp * a(i,j);
            end
            if ( nounit )
              x(j) = x(j) * a(j,j);
            end
          end
        end
      else
        kx = kx + ( n - 1 ) * incx;
        jx = kx;
        for j = n : -1 : 1
          if ( x(jx) ~= 0.0 )
            temp = x(jx);
            ix = kx;
            for i = n : -1 : j + 1
              x(ix) = x(ix) + temp * a(i,j);
              ix = ix - incx;
            end
            if ( nounit )
              x(jx) = x(jx) * a(j,j);
            end
          end
          jx = jx - incx;
        end
      end
    end
  else
%
%  Form x := A'*x.
%
    if ( lsame ( uplo, 'U' ) )
      if ( incx == 1 )
        for j = n : -1 : 1
          temp = x(j);
          if ( nounit )
            temp = temp * a(j,j);
          end
          for i = j - 1 : -1 : 1
            temp = temp + a(i,j) * x(i);
          end
          x(j) = temp;
        end
      else
        jx = kx + ( n - 1 ) * incx;
        for j = n : -1 : 1
          temp = x(jx);
          ix = jx;
          if ( nounit )
            temp = temp * a(j,j);
          end
          for i = j - 1 : -1 : 1
            ix = ix   - incx;
            temp = temp + a(i,j) * x(ix);
          end
          x(jx) = temp;
          jx = jx - incx;
        end
      end
    else
      if ( incx == 1 )
        for j = 1 : n
          temp = x(j);
          if ( nounit )
            temp = temp * a(j,j);
          end
          for i = j + 1 : n
            temp = temp + a(i,j) * x(i);
          end
          x(j) = temp;
        end
      else
        jx = kx;
        for j = 1 : n
          temp = x(jx);
          ix = jx;
          if ( nounit )
            temp = temp * a(j,j);
          end
          for i = j + 1 : n
            ix = ix + incx;
            temp = temp + a(i,j) * x(ix);
          end
          x(jx) = temp;
          jx = jx + incx;
        end
      end
    end
  end

  return
end


text:
function res = normc(A)
	%% ================== File info ==========================
	% Author		: Tiep Vu (http://www.personal.psu.edu/thv102/)
	% Time created	: Tue Jan 26 22:15:40 2016
	% Last modified	: Tue Jan 26 22:15:42 2016
	% Description	: normalize columns of a matrix (each column has Euclidean norm = 1)
	%		This is a built-in function in some new version of MATLAB
	%% ================== end File info ==========================
	B = A.^2;
	C = sqrt(sum(B));
	res = A./repmat(C, size(A,1), 1);
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function [y,zf,u]=randfilt(pb,pa,ny,zi)
%RANDFILT Generate filtered gaussian noise without initial transient
%
% Generates noise by passing gaussian noise through a filter pb/pa
% output signal y is a column vector of length ny
%
% zf and zi are the output and optional input state as defined in filter()
% If zi is not specified, random numbers with the correct covariance will be used.
% output u*u' is the state covariance matrix for filter().

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

% check to see if we must generate zi

if pa(1)~=1
   pb=pb/pa(1); pa=pa/pa(1);
end
if nargin<4 | nargout>2
   lb=length(pb);
   la=length(pa);

   k=max(la,lb)-1;
   l=la-1;
   ii=k+1-l:k;

   % form controllability matrix

   q=zeros(k,k);
   [z,q(:,1)]=filter(pb,pa,1);
   for i=2:k [z,q(:,i)]=filter(pb,pa,0,q(:,i-1)); end

   % we generate m through the step-down procedure
   s=randn(k,1);
   if l
      m=zeros(l,l);
      g=pa;
      for i=1:l
         g=(g(1)*g(1:end-1)-g(end)*g(end:-1:2))/sqrt((g(1)-g(end))*(g(1)+g(end)));
         m(i,i:l)=g;
      end
      s(ii)=triu(toeplitz(pa(1:l)))*(m\s(ii));
      if nargout>2
         u=q;
         u(:,ii)=q(:,ii)*triu(toeplitz(pa(1:l)))/m;
      end
   else
      if nargout>2
         u=q;
      end
   end
   if nargin < 4
      if k
         zi=q*s;
      else
         zi=[];
      end
   end
end
if ny
   [y,zf]=filter(pb,pa,randn(ny,1),zi);
else
   zf=zi;
   y=[];
   end

author: decouples


text:
function nx = setNull(x,e)
% set approx zero to exact zero

if nargin == 1
    e = -10;
end

epsilon = exp(log(10)*e);
nx = x;
nx(abs(nx) < epsilon) = 0;

author: mtex-toolbox

repo: mtex

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


text:
%% DEMOHESSIAN  Short demonstration of Hessians
%

%% Some sample applications of the Hessian toolbox
% Hessians implement second order automatic differentiation in forward mode, which is
% conveniently to implement using the Matlab operator concept.
%

format compact long infsup
setround(0)                           % set rounding to nearest

%% Initialization of Hessians
% The initialization follows the same principles as for gradients, for example

x = hessianinit([ -.3 ; pi ])

%% Access of the Hessian
% Define the function f: R^n->R^n by

f = @(x)( 3*x(1)*x + sin(x).*(sec(x)-sqrt(x)) )

%%
% The number of unknowns is determined by the length of the input vector x.
% For example,

f(1:4)

%%
% The function can be evaluated using the Hessian package as follows:

y = f(hessianinit([1.1 -2.7 pi]'));

%%
% There is direct access of the Hessian of y=f(x) by

y.hx

%%
% However, y.hx contains the Hessians of all three component functions of
% the original function f. To access the Hessians it is recommended to
% use the Hessian of individual components, not components of y.hx:

H3 = y(3).hx

%%
% The matrix H3 is the Hessian of the third component function of f at x.

%% A simple example
% Consider the following function

f = inline(' sin(x-log(x+2))-x*cosh(x)/15 ')

%%
% To plot the function first vectorize f :

F = vectorize(f)

%%
% Plot the function including the x-axis:

x = linspace(-1,3); close; plot( x,F(x), x,0*x )

%% Zeros of a function
% There are two roots near 1.5 and 2.5, and two extrema
% near -0.5 and 2. The roots can be included by verifynlss.
% For this simple function the inclusion is of maximum accuracy.

X1 = verifynlss(f,1.5)
X2 = verifynlss(f,2.5)

%% Extrema of a function
% The extrema can be approximated and included using Hessians.
% The following is one step of a simple Newton iteration starting at x=-0.5 :

x = -0.5;
y = f(hessianinit(x));
x = x - y.hx\y.dx';
y

%% Inclusion of extrema
% Inclusions of the extrema of f are computed by "verifynlss" with
% parameter 'h': This parameter specifies that
% f'(x) = 0 is solved instead of f(x) = 0.

X1 = verifynlss(f,-0.5,'h')
X2 = verifynlss(f,2,'h')

%% Functions in several unknowns
% Function with several unknowns are handled in a similar manner.
% Consider the following test function by N. Gould. It is taken from
% the Coconut collection of test function for global optimization,
% http://www.mat.univie.ac.at/~neum/glopt/coconut/benchmark/Library2.html .

f = inline(' x(3)-1 + sqr(x(1)) + sqr(x(2)) + sqr(x(3)+x(4)) + sqr(sin(x(3))) +  sqr(x(1))*sqr(x(2)) + x(4)-3 + sqr(sin(x(3))) + sqr(x(4)-1) + sqr(sqr(x(2))) + sqr(sqr(x(3)) + sqr(x(4)+x(1))) + sqr(x(1)-4 + sqr(sin(x(4))) + sqr(x(2))*sqr(x(3))) + sqr(sqr(sin(x(4)))) ')

%% Approximation of an extremum
% On the Web-site the global minimum of that function in 4 unknowns is
% given as 5.7444 . We use a couple of a Newton iterations
% starting at x=ones(4,1) to approximate a stationary point:

format short
x = ones(4,1);
for i=1:18
   y = f(hessianinit(x));
   x = x - y.hx\y.dx';
end
y.dx

%%
% Now x is an approximation of a stationary point:
% The gradient of f evaluated at x is not too far from zero.

%% Inclusion of a stationary point
% Using this approximation an inclusion of a stationary point of f is
% computed by (in this case the last parameter 1 is used to see
% intermediate results):

format long
X = verifynlss(f,x,'h',1)

%% Inclusion of a minimum
% The interval vector X includes a root of f', i.e. a stationary point xx of f. To prove
% that f has a minumum at xx we need to prove positive definiteness of
% the Hessian of f evaluated at xx. The interval vector X includes the stationary point
% xx of f. So we compute an inclusion Y of the Hessian at X.
%
% Mathematically,
% for every x in X the following is true: Y.x is an inclusion of f(x),
% Y.dx is an inclusion of f'(x), and Y.hx is an inclusion of the
% Hessian of f at x. Especially, the Hessian of f evaluated at xx is
% included in Y.hx.

Y = f(hessianinit(X));
format _
Y.hx

%%
% This interval matrix contains obviously only diagonally dominant
% matrices, so the stationary point xx of f in X is indeed a (local)
% minimum.
%

%% A model problem in 5000 unknowns I
% The next problem is taken from
%
% http://www.sor.princeton.edu/~rvdb/ampl/nlmodels/cute/bdqrtic.mod
%
%    Source: Problem 61 in
%       A.R. Conn, N.I.M. Gould, M. Lescrenier and Ph.L. Toint,
%       "Performance of a multifrontal scheme for partially separable optimization",
%        Report 88/4, Dept of Mathematics, FUNDP (Namur, B), 1988.
%        Copyright (C) 2001 Princeton University
%        All Rights Reserved
%    see bottom of file test_h.m
%
% The model problem is
%
%     N = length(x);      % model problem: N = 1000, initial approximation x=ones(N,1);
%     I = 1:N-4;
%     y = sum( (-4*x(I)+3.0).^2 ) + sum( ( x(I).^2 + 2*x(I+1).^2 + ...
%               3*x(I+2).^2 + 4*x(I+3).^2 + 5*x(N).^2 ).^2 );
%
% This function is evaluated by
%
%     index = 2;
%     y = test_h(x,index);

%% A model problem in 5000 unknowns II
% The given starting vector is x = ones(5000,1).
% Recall that y = f(hessianinit(x)) has 5000 elements in y.x, 2.5e7 elements
% in y.dx and 1.25e11 elements in y.hx. In full storage this would mean 1 TeraByte
% of storage.
%
% The problem can be solved in the above manner. On my 2.8 GHz
% Laptop this requires about 1 second per Newton iteration, and 5 seconds
% for one function evaluation with interval argument.
%
% The following calculates an inclusion of a stationary point of f by first
% performing a simple Newton iteration followed by a verification step for
% the nonlinear system. The Hessian is treated as full matrix, so the
% computation may take a while.

n = 5000;
index = 2;
tic
X = verifynlss('test_h',ones(n,1),'h',1,index);
tfull = toc
max(relerr(X))

%% A model problem in 5000 unknowns III
% Note the small maximum relative error of the inclusion of the result.
% Verification is faster when solving the nonlinear system treating the
% Hessian as sparse. This is done by the following.

n = 5000;
index = 2;
tic
X = verifynlss('test_h',ones(n,1),'hSparseSPD',1,index);
tsparse = toc
tfull
max(relerr(X))

%%
% Note that verification is faster at the price of a less narrow inclusion
% (for comparison the previous time tfull is displayed).
%

%% Verification of a minimum
% Finally, the Hessian at X is computed which includes the Hessian at the
% stationary point xhat in X.

y = test_h(hessianinit(X),index);
isspd(y.hx)

%%
% The interval Hessian is proved by "isspd" to be symmetric
% positive definite: Mathematically the result "isspd(M)=1" for a symmetric
% (Hermitian) interval matrix M proves that every symmetric (Hermitian) matrix
% A within M is positive definite. In our case especially the Hermitian of f
% at the stationary point xhat. Therefore, xhat is proved to be a local minimum of f.

%% Enjoy INTLAB
% INTLAB was designed and written by S.M. Rump, head of the Institute for Reliable Computing,
% Hamburg University of Technology. Suggestions are always welcome to rump (at) tuhh.de

repo: OpenMAS


text:
function nodesIdxs = nodeBForder(nodeNo,wt)
%NODEBFORDER Nodes in the Breadth-First search order
%  Usage:  nodesIdxs = nodeBForder(nodeNo,wt)
%
%   Input parameters:
%         nodeNo : Id of a node.
%         wt     : Structure containing description of the filter tree.
%
%   Output parameters:
%         nodesIdxs   : Node indexes in the Breadth-First search order.
%
%   `nodeBForder(nodeNo,wt)` For definition of the structure see
%   `wfbinit`. `nodeNo` defaults to the root node if it is empty or equal
%   to 0.
%
%
%   See also: wfbtinit
%

if isempty(nodeNo) || nodeNo==0
   %find root
   nodeNo = find(wt.parents==0);
end

complainif_notposint(nodeNo,'NODEBFORDER');

nodesIdxs = [nodeNo,nodeSubtreeBF(nodeNo,wt)];

function nodesIdxs = nodeSubtreeBF(nodeNo,wt)
%NODESUBTREEBF Node subtree nodes in Breath-First order
%   Usage:  noOut = nodeSubtreeBF(nodeNo,wt);
%
%   Input parameters:
%         nodeNo  : Node index.
%         wt      : Structure containing description of the filter tree.
%
%   Output parameters:
%         noOut   : Nodes in a Breath-First order.
%

% subtreeIdx = [];
%
% children = treeStruct.children{nodeNo}(find(treeStruct.children{nodeNo}~=0));
% subtreeIdx(end+1:end+length(children)) = children;
%
% for ii=1:length(children)
%    tmpSbIdx = nodeSubtreeBF(children(ii),treeStruct);
%    subtreeIdx(end+1:end+length(tmpSbIdx)) = tmpSbIdx;
% end

toGoTrough = [nodeNo];
nodesIdxs = [];
while ~isempty(toGoTrough)
  % chtmp = find(wt.children{toGoTrough(1)}~=0);
   chIdxtmp = wt.children{toGoTrough(1)}(wt.children{toGoTrough(1)}~=0);
   nodesIdxs = [nodesIdxs,chIdxtmp];
   toGoTrough = [toGoTrough(2:end),chIdxtmp];
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
classdef LevelSetFactory < handle

    methods (Access = public, Static)

        function obj = create(d)
            switch d.type
                case 'circle'
                    obj = LevelSetCircle(d);
                case 'circleInclusion'
                    obj = LevelSetWithCircleInclusion(d);
                case 'sphere'
                    obj = LevelSetSphere(d);
                case 'sphereInclusion'
                    obj = LevelSetWithSphereInclusion(d);
                case 'cylinder'
                    obj = LevelSetCylinder(d);
                case 'horizontal'
                    obj = LevelSetHorizontalInclusion(d);
                case {'squareInclusion'}
                    obj = LevelSetSquareInclusion(d);
                case 'smoothSquare'
                    obj = LevelSetSmoothSquareInclusion(d);
                case 'rectangle'
                    obj = LevelSetRectangle(d);
                case 'rectangleInclusion'
                    obj = LevelSetRectangleInclusion(d);
                case 'smoothRectangle'
                    obj = LevelSetSmoothRectangleInclusion(d);
                case 'feasible'
                    obj = LevelSetFeasible(d);
                case 'rand'
                    obj = LevelSetRandom(d);
                case 'holes'
                    obj = LevelSetWithSeveralHoles(d);
                case 'full'
                    obj = LevelSetFull(d);
                case 'horizontalFibers'
                    obj = LevelSetHorizontalFibers(d);
                case 'given'
                    obj = LevelSetGiven(d);
                case 'Vigdergauz'
                    obj = LevelSetVigdergauz(d);
                case 'periodicAndOriented'
                    obj = LevelSetPeriodicAndOriented(d);
                otherwise
                    error('Invalid initial value of design variable.');
            end

        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function varargout = sendemailGUI(varargin)
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @sendemailGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @sendemailGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end

%---------------------------------------------------------------------------------------------------
function sendemailGUI_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = [];

try
        def = varargin{1};
catch
        def = {'','You got an ERPmail!', 1};
end

memail = def{1}';
nm = length(memail);
memmailto = '';
for k=1:nm
        memmailto = sprintf('%s %s', memmailto, memail{k});
end
memmailto  = strtrim(memmailto);
memmailto  = strrep(memmailto, ' ', ', ');
memsubject = def{2};
memerpset  = def{3};
set(handles.edit_to, 'String', memmailto)
set(handles.edit_subject, 'String', memsubject)
set(handles.edit_erpsets, 'String', num2str(memerpset))

%
% Name & version
%
version = geterplabversion;
set(handles.gui_chassis,'Name', ['ERPLAB ' version '   -   Send Email GUI'])
fromemail = getpref('Internet','SMTP_Username');
if isempty(fromemail)
      fromemail = 'email not found';
end
set(handles.text_myemail, 'String', fromemail)

%
% Color GUI
%
handles = painterplab(handles);

%
% Set font size
%
handles = setfonterplab(handles);

% Update handles structure
guidata(hObject, handles);

% help
% helpbutton

% UIWAIT makes sendemailGUI wait for user response (see UIRESUME)
uiwait(handles.gui_chassis);

%---------------------------------------------------------------------------------------------------
function varargout = sendemailGUI_OutputFcn(hObject, eventdata, handles)
% Get default command line output from handles structure
varargout{1} = handles.output;

% The figure can be deleted now
delete(handles.gui_chassis);
pause(0.1)

%---------------------------------------------------------------------------------------------------
function edit_to_Callback(hObject, eventdata, handles)

%---------------------------------------------------------------------------------------------------
function edit_to_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%---------------------------------------------------------------------------------------------------
function edit_subject_Callback(hObject, eventdata, handles)

%---------------------------------------------------------------------------------------------------
function edit_subject_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%---------------------------------------------------------------------------------------------------
function edit_erpsets_Callback(hObject, eventdata, handles)

%---------------------------------------------------------------------------------------------------
function edit_erpsets_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%---------------------------------------------------------------------------------------------------
function pushbutton_cancel_Callback(hObject, eventdata, handles)
handles.output = [];
% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);
%---------------------------------------------------------------------------------------------------
function pushbutton_send_Callback(hObject, eventdata, handles)

mailto  = get(handles.edit_to, 'String');
mailto  = regexp(mailto, ',|\s+', 'split');
mailto(cellfun(@isempty,mailto)) = [];
subject = get(handles.edit_subject, 'String');
msg     = cellstr(char(get(handles.edit_editor, 'String')));
attach  = str2num(get(handles.edit_erpsets, 'String'));

handles.output = {mailto, subject, msg, attach};
%sendmail(mailto, subject, msg);
%pause(0.2)

% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);

%---------------------------------------------------------------------------------------------------
function edit_editor_Callback(hObject, eventdata, handles)

%---------------------------------------------------------------------------------------------------
function edit_editor_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%--------------------------------------------------------------------------
function gui_chassis_CloseRequestFcn(hObject, eventdata, handles)

if isequal(get(handles.gui_chassis, 'waitstatus'), 'waiting')
      % The GUI is still in UIWAIT, us UIRESUME
      uiresume(handles.gui_chassis);
else
      % The GUI is no longer waiting, just close it
      delete(handles.gui_chassis);
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
% ========================================================================
% Fast Multi-Scale Structural Patch Decomposition for Multi-Exposure Image Fusion, TIP,2020
% algorithm Version 1.0
% Copyright(c) 2020, Hui Li, Kede Ma, Yongwei Yong and Lei Zhang
% All Rights Reserved.
% ----------------------------------------------------------------------
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is hereby
% granted, provided that this copyright notice and the original authors'
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining permission of the
% authors. The authors make no representations about the suitability of
% this software for any purpose. It is provided "as is" without express
% or implied warranty.
%----------------------------------------------------------------------
% Please refer to the following paper:
% H. Li et al., "Fast Multi-Scale Structural Patch Decomposition for Multi-Exposure Image Fusion, 2020" In press
% IEEE Transactions on Image Processing
% Please kindly report any suggestions or corrections to xiaohui102788@126.com
%----------------------------------------------------------------------

%% motion detection in dynamic scenes

function imgSeqcor= detection_color(imgSeqColor,r)

[h,w,~,n]=size(imgSeqColor);
N = boxfilter(ones(h, w), r);
i_mean=zeros(h,w,n);

C = 0.03 ^ 2 / 2;
structureThres=0.8;
consistencyThres=0.2;
exposureThres=0.02;
ref=ceil(n/2);

%   ref = selectRef(imgSeqColor);
se = strel('disk',2*r+1);
% se = strel('square',10);

R_img=imgSeqColor(:,:,:,ref);
R_mean=boxfilter(R_img(:,:,1), r)./ N+boxfilter(R_img(:,:,2), r)./ N+boxfilter(R_img(:,:,3), r)./ N;
R_mean= R_mean./3;

R_var= boxfilter(R_img(:,:,1).*R_img(:,:,1), r)./ N+ boxfilter(R_img(:,:,2).*R_img(:,:,2), r)./ N+...
    boxfilter(R_img(:,:,3).*R_img(:,:,3), r)./ N- R_mean.* R_mean*3;

R_var = R_var./3;
R_std=sqrt(max(R_var,0));

i_mean(:,:,ref)=R_mean;
muIdxMap = R_mean < exposureThres | R_mean > 1 - exposureThres;
% R=2*r+1;
% countWindowt = ones(R, R);
% countWindow=countWindowt./sum(countWindowt(:));
i_std=zeros(h,w,n);
sRefMap=zeros(h,w,n);
imgSeqcor=zeros(h,w,3,n);

for i = 1 : n
    if i ~= ref
        img= imgSeqColor(:,:,:,i);

        i_mean(:,:,i)=boxfilter(img(:,:,1), r)./ N+boxfilter(img(:,:,2), r)./ N+boxfilter(img(:,:,3), r)./ N;
        i_mean(:,:,i)= i_mean(:,:,i)./3;

        i_var= boxfilter(img(:,:,1).*img(:,:,1), r)./ N+ boxfilter(img(:,:,2).*img(:,:,2), r)./ N+...
            boxfilter(img(:,:,3).*img(:,:,3), r)./ N- i_mean(:,:,i).* i_mean(:,:,i)*3;

        i_var = i_var./3;
        i_std(:,:,i)=sqrt(max(i_var,0));

        mean_iR=boxfilter(img(:,:,1).*R_img(:,:,1), r)./ N+boxfilter(img(:,:,2).*R_img(:,:,2), r)./ N+...
            boxfilter(img(:,:,3).*R_img(:,:,3), r)./ N;
        mean_iR=mean_iR/3;
        cov_iR = mean_iR -  i_mean(:,:,i).* R_mean;

        %         mean_iR=boxfilter(img(:,:,1).*R_img(:,:,1), r)./ N+boxfilter(img(:,:,2).*R_img(:,:,2), r)./ N+...
        %             boxfilter(img(:,:,3).*R_img(:,:,3), r)./ N-i_mean(:,:,i).* R_mean*3;
        %         cov_iR=mean_iR./3;
        %

        sMap= (cov_iR + C) ./ (R_std.* i_std(:,:,i) + C);
        sRefMapt=(sMap >structureThres);

        %         t= sRefMap(:,:,i) ;
        %         t(muIdxMap) = 1;

        sRefMapt(muIdxMap) = 1;
        %         sRefMap(:,:,i) =sRefMapt;

        sRefMap(:,:,i) = imopen(sRefMapt,se);

        %     sRefMap(:,:,i) =bwareaopen( sRefMap(:,:,i),150);
        %  sRefMap(:,:,i) = imopen(sRefMap(:,:,i),se);

        %   figure,imshow( sRefMap(:,:,i))

        %   sRefMap(:,:,i) = imopen(sRefMap(:,:,i),se);

        %
        cMu  = imhistmatch(i_mean(:,:,i), R_mean,256);
        diff = abs(cMu - R_mean);
        cMap = diff <= consistencyThres;
        sRefMap(:,:,i) =  sRefMap(:,:,i).*cMap;

        %     figure,imshow( sRefMap(:,:,i))

        %     sRefMap(:,:,i) =bwareaopen( sRefMap(:,:,i),150);

        %         figure,imshow( sRefMap(:,:,i))

        %         sRefMap(:,:,i) = imopen(sRefMap(:,:,i),se);

        temp = imhistmatch(imgSeqColor(:,:,:,ref),...
            imgSeqColor(:,:,:,i),256);
        %         temp(temp<0) = 0;
        %         temp(temp>1) = 1;

        %         figure,imshow(temp)

        imgSeqcor(:,:,:,i)=imgSeqColor(:,:,:,i).*repmat(sRefMap(:,:,i),[1 1 3])+...
            temp.*repmat(1-sRefMap(:,:,i),[1 1 3]);
    end
    imgSeqcor(:,:,:,ref)=imgSeqColor(:,:,:,ref);

end

author: thfylsty


text:
function images = publish3View(vw, savePath,plotTitle)
%
% images = publish3View(vw, [savePath], [plotTitle]);
%
% Shows the selected sagittal, coronal, and axial slices
% in a 3View window in a separate figure, with colorbar
% if relevant.
%
% If the optional 'savePath' argument is provided, exports
% the figure with the given image path.
%
% Returns the images in a cell array.
%
% ras, 01/2007.
if notDefined('vw'),      vw = getCurView;      end
if notDefined('savePath'),  savePath = '';          end

% check for 3-view
if ~isfield(vw, 'loc') || isempty(vw.loc)
    error('Not a 3-view or need a location defined in view.loc field.');
end

% get color map
if checkfields(vw, 'ui', 'displayMode')
    mode = sprintf('%sMode', vw.ui.displayMode);
    cmap = vw.ui.(mode).cmap;
end

% ROI prefs
if length(vw.ui.showROIs) <= 1
    switch vw.ui.showROIs
        case -2, roiList = 1:length(vw.ROIs);
        case -1, roiList = vw.selectedROI;
        case 0, roiList = [];
        otherwise, roiList = vw.ui.showROIs;  % manually specify
    end
else
    roiList = vw.ui.showROIs;
end

if roiList==0, roiList = []; end

% get images
[vw axi] = recompute3ViewImage(vw, vw.loc(1), 1);
[vw cor] = recompute3ViewImage(vw, vw.loc(2), 2);
[vw sag] = recompute3ViewImage(vw, vw.loc(3), 3);
axi = ind2rgb(axi, cmap);
cor = ind2rgb(cor, cmap);
sag = ind2rgb(sag, cmap);
images = {sag cor axi};

% open figure
hFig = figure('Color', 'w', 'Position', [0 300 600 400]);
centerfig(hFig, 0);
colormap(cmap);

% sagittal
subplot('Position', [0 0 1/3 1]);
image(sag); axis image; axis off;
drawROIs3View(vw, roiList, gca, 3, vw.loc)
axis([vw.ui.zoom(2,:) vw.ui.zoom(1,:)]);
sagAxesHandle = gca;

% corornal
subplot('Position', [1/3 0 1/3 1]);
image(cor); axis image; axis off;
drawROIs3View(vw, roiList, gca, 2, vw.loc)
axis([vw.ui.zoom(3,:) vw.ui.zoom(1,:)]);
corAxesHandle = gca;

% axial
subplot('Position', [2/3 0 1/3 1]);
image(axi); axis image; axis off;
drawROIs3View(vw, roiList, gca, 1, vw.loc)
axis([vw.ui.zoom(3,:) vw.ui.zoom(2,:)]);
axiAxesHandle = gca;

% Add crosshairs. The function to render crosshairs requires a vw structure
% as input. In order to add cross hairs to the figure to publish, we add
% the 3 axis handles (axial, sagittal, and coronal) to the view struct ui,
% and then call renderCrosshairs. This is a bit of a hack, but it seems
% nicer than writing a new function that renders cross hairs only for
% published figures.
vw.ui.axiAxesHandle = axiAxesHandle;
vw.ui.sagAxesHandle = sagAxesHandle;
vw.ui.corAxesHandle = corAxesHandle;
renderCrosshairs(vw, 1);

if checkfields(vw, 'ui', 'displayMode') && ...
    ismember(vw.ui.displayMode, {'co' 'amp' 'ph' 'map'})

    % put up color bar
    addCbarLegend(vw, hFig, .16);
end

if exist('plotTitle', 'var'), suptitle(plotTitle); end

if ~isempty(savePath)
    saveas(hFig, savePath);
%     close(hFig);
end

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

% /----------------------------------------------------------------/ %
function drawROIs3View(vw, roiList, axs, ori, loc)
% draw the ROIs, parsing the relevant settings, on the axes
% specified by [axs]. This is implemented a little differently
% from drawROIs or drawROIsPerim, in that it doesn't select a given
% orientation each time (which caused some extra time during refreshes).
% roiList is an array into the vw.ROIs field.
ui = vw.ui;

% build prefs for the outline function.
prefs.method = ui.roiDrawMethod;
prefs.axesHandle = axs;

for r = roiList
    R = vw.ROIs(r);
    if ~isempty(R.coords)
        prefs.color = R.color;
        if r==vw.selectedROI, prefs.color=viewGet(vw,'selRoiColor'); end

        switch ori
            case 1, pts = R.coords([2 3], R.coords(1,:)==loc(1));
            case 2, pts = R.coords([1 3], R.coords(2,:)==loc(2));
            case 3, pts = R.coords([1 2], R.coords(3,:)==loc(3));
        end

        if isfield(ui, 'flipLR') && ui.flipLR==1 && ori<3
            % L/R flip affects columns of axi + coronal, but not sag, orientations
            dims = viewGet(vw,'Size');
            pts(2,:) = dims(3) - pts(2,:);
        end

        % restrict to voxels within zoom range

        h = outline(pts, prefs);
        if ishandle(h)
            set(h, 'ButtonDownFcn', sprintf('recenter3View(%s,%i);',vw.name, ori));
        end
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
target_vy = [-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8];
% target_vy = [0];
target_vx = [-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4];
% target_vx = [0];
subfolder_name = 'library5';

N_mid = 11;%floor(nlp.Phase(1).NumNode/2)+1;
N_vx = length(target_vx);
N_vy = length(target_vy);
% gait = param.gait;
t_all = cell(N_vx, N_vy);
q_all = cell(N_vx, N_vy);
dq_all = cell(N_vx, N_vy);
px_all = cell(N_vx, N_vy);
py_all = cell(N_vx, N_vy);
vx_all = cell(N_vx, N_vy);
vy_all = cell(N_vx, N_vy);
for i = 1:N_vx
    vx = target_vx(i);
    for j = 1:N_vy
        vy = target_vy(j);
        data_name = fullfile('local', subfolder_name, sprintf('gait_X%0.1f_Y%.1f.mat', vx, vy));
        param = load(data_name);

        t_all{i,j} = [param.gait(3).tspan];%, param.gait(3).tspan];
        q_all{i,j} = [param.gait(3).states.x];%, param. gait(3).states.x];
        dq_all{i,j} = [param.gait(3).states.dx];%, param. gait(3).states.dx];
        px_all{i,j} = [param.gait(3).states.x(1,:)];%ones(size(t_all{i,j}))*vx;
        py_all{i,j} = [param.gait(3).states.x(2,:)];%ones(size(t_all{i,j}))*vy;
        vx_all{i,j} = [param.gait(3).states.dx(1,:)];%ones(size(t_all{i,j}))*vx;
        vy_all{i,j} = [param.gait(3).states.dx(2,:)];%ones(size(t_all{i,j}))*vy;

    end

end

%%
idx = [3:6 7:9, 12:14];
joint_names = {    'BasePosX'
    'BasePosY'
    'BasePosZ'
    'BaseRotX'
    'BaseRotY'
    'BaseRotZ'
    'qHRight'
    'qARight'
    'qBRight'
    'fourBarARight'
    'fourBarBRight'
    'qHLeft'
    'qALeft'
    'qBLeft'
    'fourBarALeft'
    'fourBarBLeft'};
for k=1:length(idx)
    f = figure(k+100); clf;
    f.Name = joint_names{idx(k)};
    set(f, 'WindowStyle', 'docked');

    ax = axes(f); %#ok<LAXES>
    hold(ax);

    for i = 1:N_vx
        for j = 1:N_vy
            t = t_all{i,j};
            q = q_all{i,j};
            dq = dq_all{i,j};
            py = py_all{i,j};
            px = px_all{i,j};
            vy = vy_all{i,j};
            vx = vx_all{i,j};
            scatter3(ax,t, vx,dq(idx(k),:));
%             scatter3(ax,t,vx,q(idx(k),:));

        end

    end

end
%%
idx = [7:9, 12:14];
for k=1:length(idx)
    f = figure(k+200); clf;
    f.Name = joint_names{idx(k)};
    set(f, 'WindowStyle', 'docked');

    ax = axes(f); %#ok<LAXES>
    hold(ax);

    for i = 1:N_vx
        for j = 1:N_vy
            t = t_all{i,j};
            q = q_all{i,j};
            dq = dq_all{i,j};
            py = py_all{i,j};
            px = px_all{i,j};
            vy = vy_all{i,j};
            vx = vx_all{i,j};
            plot3(ax, vx(N_mid), vy(N_mid), q(idx(k),N_mid),'*','MarkerSize',4);
%             scatter3(ax,t,vx,q(idx(k),:));
            [S, V] = meshgrid(-0.5:0.01:0.5, -1:0.02:1);
            L = zeros(size(S));
            for ii = 1:size(S, 1)
                for jj = 1:size(S, 2)
                    L(ii, jj) = P(k,:)*[1; S(ii, jj); V(ii, jj)];
                    %L(i, j) = l(n);
                end
            end
            surface(ax, S, V, L);
        end

    end

end

%%
idx = [7:9, 12:14];
X = zeros(2, N_vx*N_vy);
Y = zeros(length(idx), N_vx*N_vy);
dY = zeros(length(idx), N_vx*N_vy);
ii = 1;
for i = 1:N_vx
    for j = 1:N_vy
        t = t_all{i,j};
        q = q_all{i,j};
        dq = dq_all{i,j};

        %         X(:,ii) = [q(1:2,N_mid);dq(1:2,N_mid)];
        X(:,ii) = dq(1:2,N_mid);
        Y(:,ii) = q(idx,N_mid);
        dY(:,ii) = dq(idx,N_mid);
        ii = ii+1;
    end
end
% x = [ones(1,45); X];
x = [ones(1,81); X];
P = Y/x;
dP = dY/x;

%%
y1 = dY(4,:);
x1 = X(1,:);
x2 = X(2,:);
%%
0.0008997

0.06608

%% Save Training Data
save('Ayonga3DinsertionFunction', 'P','dP')

%%
% P = zeros(length(idx), 3);
% dP = zeros(length(idx), 3);
%
% for i=1:length(idx)
%
%     P(:,i) = Y(i,:)/X;
%     dP(:,i) = dY(i,:)/X;
% end

%%
f = figure(2);clf;
f.Name = 'Sideway Velocity Feature';
set(f, 'WindowStyle', 'docked');
ax = axes(f);
hold(ax);
for i = 1:N_vx
    t = t_all{i,5};
    dq = dq_all{i,5};
    plot(ax, t, dq(1,:));
    plot(ax, t(N_mid), dq(1,N_mid),'*','MarkerSize',4);
end
grid on
%%
f = figure(3);clf;
f.Name = 'Forward Velocity Feature';
set(f, 'WindowStyle', 'docked');
ax = axes(f);
hold(ax);
for j = 1:N_vy
    t = t_all{3,j};
    dq = dq_all{3,j};
    plot(ax, t, dq(2,:));
    plot(ax, t(N_mid), dq(2,N_mid),'*','MarkerSize',4);
end
grid on
%%
f = figure(4);clf;
f.Name = 'Sideway Position Feature';
set(f, 'WindowStyle', 'docked');
ax = axes(f);
hold(ax);
for i = 1:N_vx
    t = t_all{i,5};
    q = q_all{i,5};
    plot(ax, t, q(1,:));
    plot(ax, t(N_mid), q(1,N_mid),'*','MarkerSize',4);
end
grid on
%%
f = figure(5);clf;
f.Name = 'Forward Position Feature';
set(f, 'WindowStyle', 'docked');
ax = axes(f);
hold(ax);
for j = 1:N_vy
    t = t_all{3,j};
    q = q_all{3,j};
    plot(ax, t, q(2,:));
    plot(ax, t(N_mid), q(2,N_mid),'*','MarkerSize',4);
end
grid on

%%
f = figure(6);clf;
f.Name = 'Sideway Position Feature';
set(f, 'WindowStyle', 'docked');
ax = axes(f);
hold(ax);
for i = 1:N_vx
    q = q_all{i,5};
    dq = dq_all{i,5};
    plot(ax, q(1,:), dq(1,:));
    plot(ax, q(1,N_mid), dq(1,N_mid),'*','MarkerSize',4);
end
grid on
%%
f = figure(7);clf;
f.Name = 'Forward Position Feature';
set(f, 'WindowStyle', 'docked');
ax = axes(f);
hold(ax);
for j = 1:N_vy
    q = q_all{3,j};
    dq = dq_all{3,j};
    plot(ax, q(2,:), dq(2,:));
    plot(ax, q(2,N_mid), dq(2,N_mid),'*','MarkerSize',4);
end
grid on

author: ayonga

repo: frost-dev

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


text:
function x=complexBarycentricCoords2Pt(c,v)
%%COMPLEXBARYCENTRICCOORDS2PT Given a set of complex barycentric
%   coordinates and a st of vertices, get the real 2D point that is
%   specified by the coordinates.
%
%INPUTS: c A numVertXnumPts set of complex coordinate vectors to convert.
%        v This is either a 2XN set of real vertices, or a length N vctor
%          of complex scalar vertices (x coordinate is real and y
%          coordinate is imaginery). These define the polygon over which
%          the vertices are specified.
%
%OUTPUTS: x The 2XnumPts set of the points given by the complex barycentric
%           coordinates.
%
%Complex barycentric coordinates are discussed in [1].
%
%REFERENCES:
%[1] O. Weber, M. Ben-Chen, C. Gotsman, and K. Hormann, "A complex view of
%    barycentric mappings," Computer Graphics Forum, vol. 30, no. 5, pp.
%    1533-1542, Aug. 2011.
%
%October 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numPts=size(c,2);

if(isvector(v))
    zj=v;
else
    zj=v(1,:).'+1j*v(2,:).';
end

x=zeros(2,numPts);
for curPt=1:numPts
    xComplex=sum(c(:,curPt).*zj);
    x(:,curPt)=[real(xComplex);imag(xComplex)];
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 ph = atlasUndoPhase(atlasParams, ph)
%
%      ph = atlasUndoPhase(atlasParams, ph)
%
% Purpose:
%   Convert the atlas phase data by shifting their phase and scaling the
%   phase range.  The amount of shifting and scaling is stored in the
%   phaseShift/phaseScale fields.  I am not sure why this code is used at
%   all.  I think we should be using shiftPhase instead and this is left
%   over from version 1.0.  I am trying to get rid of it -- BW
%
% HISTORY:
% 2003.08.29 RFD (bob@white.stanford.edu) got tired of duplicating this
% code, so here's a simple function.

warning('atlasUndoPhase is obsolete.')

% Checking.
if(length(atlasParams) ~= length(ph) | length(atlasParams(1).phaseShift) ~= size(ph{1},3))
    warning('atlasParams and ph should have the same number of scans and slices!');
end

%  Create a scaled and shifted version of the atlas parameters.
for(scan=1:length(atlasParams))
    for(slice=1:length(atlasParams(scan).phaseShift))
        ph{scan}(:,:,slice) = mod((squeeze(ph{scan}(:,:,slice)) ...
            - atlasParams(scan).phaseShift(slice)) / atlasParams(scan).phaseScale(slice), 2*pi);
    end
end

return;


text:
function layer = genNetworkLSTM(para)
if isfield(para, 'LastActivation4MSE')==0
    para.LastActivation4MSE = 'linear';
end

layer{1}.name = 'Input';        % this is an input layer
layer{end}.inputIdx = 1;    % specifies the index of GCC in Visible_tr
layer{end}.dim = [1 1]*double(para.inputDim);             % [input dim; output dim];

if isfield(para, 'ProjectionSize') &&  ~isempty(para.ProjectionSize) && para.ProjectionSize>0
    layer{end+1}.name = 'Affine';
    layer{end}.prev = -1;
    layer{end}.W = [];
    layer{end}.b = [];
    layer{end}.dim = [para.ProjectionSize layer{end-1}.dim(1)];
    layer{end}.update = 1;
end

for i=1:length(para.hiddenLayerSizeLSTM)
    layer{end+1}.name = 'LSTM';
    layer{end}.prev = -1;   % this is the index of GCC features in Visible_tr; It is the offset to be added to current layer index
    layer{end}.W = []; % to be initialized randomly or by pretraining
    layer{end}.b = [];
    layer{end}.usePastState = para.usePastState(i);
    layer{end}.dim = [para.hiddenLayerSizeLSTM(i) layer{end-1}.dim(1)];
    layer{end}.update = 1;

    if isfield(para, 'useAffineBtwLSTM') && para.useAffineBtwLSTM && i<length(para.hiddenLayerSizeLSTM)
        layer{end+1}.name = 'Affine';
        layer{end}.prev = -1;
        layer{end}.W = [];
        layer{end}.b = [];
        layer{end}.dim = [1 1] * layer{end-1}.dim(1);
        layer{end}.update = 1;
    end
end

layer2 = genNetworkFeedForward_v2(layer{end}.dim(1), para.hiddenLayerSizeFF, para.outputDim, para.costFn, para.LastActivation4MSE);

layer = [layer layer2(2:end)];

if isfield(para, 'labelDelay')
    layer{end}.labelDelay = para.labelDelay;
end
if isfield(para, 'costFrameSelection');
    layer{end}.costFrameSelection = para.costFrameSelection;
end

layer = FinishLayer(layer);
end

author: singaxiong

repo: SignalGraph


text:
function test_ft_plot_crosshair

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_plot_crosshair

close all

result = {};
figure; ft_plot_axes([], 'unit', 'mm'); ft_plot_crosshair([100 100 100], 'color', 'm'); view([1 0.5 0.5]); result{end+1} = getframe;
figure; ft_plot_axes([], 'unit', 'mm'); ft_plot_crosshair([100 100 100], 'color', 'y'); view([1 0.5 0.5]); result{end+1} = getframe;
figure; ft_plot_axes([], 'unit', 'mm'); ft_plot_crosshair([100 100 100], 'color', [0 0 0]); view([1 0.5 0.5]); result{end+1} = getframe;
figure; ft_plot_axes([], 'unit', 'mm'); ft_plot_crosshair([100 100 100], 'color', 'powderblue'); view([1 0.5 0.5]); result{end+1} = getframe;
figure; ft_plot_axes([], 'unit', 'mm'); ft_plot_crosshair([100 100 100], 'color', 'limegreen'); view([1 0.5 0.5]); result{end+1} = getframe;

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequal(result{i}, result{j}));
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function disp(F)
%DISP   Display a DISKFUN to the command line.
%
% See also DISKFUN/DISPLAY.

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

loose = strcmp( get(0, 'FormatSpacing'), 'loose' );

% Get display style and remove trivial empty DISKFUN case.
if ( isempty(F) )
    fprintf('    empty diskfun\n')
    if ( loose )
        fprintf('\n');
    end
    return
end

% Get information that we want to display:
len = length(F);                          % Numerical rank
vscl = vscale(F);                         % vertical scale

% Display the information:
disp('     diskfun object ')
fprintf('       domain        rank    vertical scale\n');
fprintf('      unit disk   %6i          %3.2g\n', len, vscl);

if ( loose )
    fprintf('\n');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Fig. 5.8   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami

clf
n9=1;
d9=[1 8 32 0];
rlocus(n9,d9)
title('Fig.5.08 The root locus for L(s)=1/s(s^2+8s+32)')
axis([-14 2 -6 6])
z=0:.1:.9;
 wn= 1:6;
 sgrid(z, wn)
hold on
x=[-2.67 2 -2.67 2];
y=[0 4.67*sqrt(3) 0 -4.67*sqrt(3)];
plot(x,y)
r=roots([1 0 32]);
 plot(r,'*')
 hold off

author: Sable

repo: mcbench-benchmarks


text:
function [cc] = cdtbal2(pp,ee,tt)
%CDTBAL2 compute the modified circumballs associated with a
%constrained 2-simplex Delaunay triangulation in R^2.
%   [CC] = CDTBAL2(PP,EE,TT) returns the smallest enclosing
%   balls associated with the triangles in [PP,TT], such th-
%   at CC = [XC,YC,RC.^2]. Such balls never lie outside the
%   boundaries of the associated CDT. See TRICON2 for info-
%   mation regarding the edge array EE.

%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 01/10/2017

%---------------------------------------------- basic checks
    if (~isnumeric(pp) || ~isnumeric(ee) || ...
        ~isnumeric(tt) )
        error('cdtbal2:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

%---------------------------------------------- basic checks
    if (ndims(pp) ~= +2 || ndims(ee) ~= +2 || ...
        ndims(tt) ~= +2 )
        error('cdtbal2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end
    if (size(pp,2)~= +2 || size(ee,2) < +5 || ...
        size(tt,2) < +6 )
        error('cdtbal2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

%----------------------------------------- calc. circumballs
    cc = tribal2(pp,tt);

%------------------------ replace with face-balls if smaller
    cc = minfac2(cc,pp,ee,tt,1,2,3) ;
    cc = minfac2(cc,pp,ee,tt,2,3,1) ;
    cc = minfac2(cc,pp,ee,tt,3,1,2) ;

end

function [cc] = minfac2(cc,pp,ee,tt,ni,nj,nk)
%MINFAC2 modify the set of circumballs to constrain centres
%to the boundaries of the CDT.
%   [CM] = MINFAC2(CC,PP,EE,TT,NI,NJ,NK) returns the set of
%   modified circmballs CM, where any ball CC lying outside
%   the boundaries of the CDT [PP,EE,TT] is replaced by the
%   edge-centred diametric ball. [NI,NJ] are the local inde-
%   xes associated with an edge to test. NK is the local in-
%   dex of the opposite vertex.

%------------------------------------------------ outer edge
    EF = ee(tt(:,ni+3),5) > +0 ;

%------------------------------------------------ edge balls
    bc = (pp(tt(EF,ni),:)+pp(tt(EF,nj),:))*.50;

%------------------------------------------------ edge radii
    br = sum((bc(:,1:2)-pp(tt(EF,ni),:)).^2,2)...
       + sum((bc(:,1:2)-pp(tt(EF,nj),:)).^2,2);
    br = br * +0.5 ;

%------------------------------------------- enclosing radii
    ll = sum((bc(:,1:2)-pp(tt(EF,nk),:)).^2,2);

%------------------------------------------- replace if min.
    bi = br >= ll ...
       & br <= cc(EF,3) ;
    ei = find(EF) ;
    ti = ei  (bi) ;

%------------------------------------------- replace is min.
    cc(ti,1:2) = bc(bi,:) ;
    cc(ti,  3) = br(bi,:) ;

end

author: CHLNDDEV


text:
% Mat file utility functions
%
% Author: Jonathan Karr
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 11/23/2010
classdef MatFileUtil
    methods (Static)
        function diff(file1, file2)
            struct1 = load(file1);
            struct2 = load(file2);

            edu.stanford.covert.util.StructUtil.diff(struct1, struct2);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function comet(f, g, ignored)
%COMET   Two-dimensional comet plot.
%   COMET(F) displays a comet graph of the CHEBFUN F, and COMET(F, G) displays a
%   comet graph of the CHEBFUN F versus the CHEBFUN G.
%
%   A comet graph is an animated graph in which a thick dot (the comet head)
%   traces the data points on the screen.
%
% See also COMET3.

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

if ( nargin == 1 ) % COMET(F)
    data = plotData(f);
else               % COMET(F, G)
    if ( nargin > 2 )
        warning('CHEBFUN:CHEBFUN:comet:nargin', ...
            'Third input to @CHEBFUN/COMET() is ignored.');
    end
    data = plotData(f, g);
end

comet(data.xLine, data.yLine, 0);

end


text:
function laplacian_test06 ( )

%*****************************************************************************80
%
%% LAPLACIAN_TEST06 tests L1DD_LU and similar routines.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAPLACIAN_TEST06\n' );
  fprintf ( 1, '  Compute LU factors for the Laplacian:\n' );
  fprintf ( 1, '  L1DD_LU for Dirichlet/Dirichlet BC;\n' );
  fprintf ( 1, '  L1DN_LU for Dirichlet/Neumann BC;\n' );
  fprintf ( 1, '  L1ND_LU for Neumann/Dirichlet BC;\n' );
  fprintf ( 1, '  L1NN_LU for Neumann/Neumann BC;\n' );
  fprintf ( 1, '  L1PP_LU for Periodic BC;\n' );

  n = 5;

  for test = 1 : 2

    if ( test == 1 )
      h = 1.0;
    else
      h = 1.0 / ( n + 1 );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using spacing H = %g\n', h );

    a = l1dd ( n, h );
    [ l, u ] = l1dd_lu ( n, h );
    r8mat_print ( n, n, l, '  L1DD L factor:' );
    r8mat_print ( n, n, u, '  L1DD U factor:' );
    err = lu_error ( n, a, l, u );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  L1DD LU error = %g\n', err );

    a = l1dn ( n, h );
    [ l, u ] = l1dn_lu ( n, h );
    r8mat_print ( n, n, l, '  L1DN L factor:' );
    r8mat_print ( n, n, u, '  L1DN U factor:' );
    err = lu_error ( n, a, l, u );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  L1DN LU error = %g\n', err );

    a = l1nd ( n, h );
    [ l, u ] = l1nd_lu ( n, h );
    r8mat_print ( n, n, l, '  L1ND L factor:' );
    r8mat_print ( n, n, u, '  L1ND U factor:' );
    err = lu_error ( n, a, l, u );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  L1ND LU error = %g\n', err );

    a = l1nn ( n, h );
    [ l, u ] = l1nn_lu ( n, h );
    r8mat_print ( n, n, l, '  L1NN L factor:' );
    r8mat_print ( n, n, u, '  L1NN U factor:' );
    err = lu_error ( n, a, l, u );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  L1NN LU error = %g\n', err );

    a = l1pp ( n, h );
    [ l, u ] = l1pp_lu ( n, h );
    r8mat_print ( n, n, l, '  L1PP L factor:' );
    r8mat_print ( n, n, u, '  L1PP U factor:' );
    err = lu_error ( n, a, l, u );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  L1PP LU error = %g\n', err );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function T56=T56(h)

T56=270.5;


text:
function E = compose(disc, y)
%FEVAL   Evaluation functional for CHEBCOLLOC.
%   FEVAL(DISC, LOC, DIRN) returns a functional that evaluates the Chebyshev
%   polynomial represented by a COLLOC discretization at the given point LOC as
%   approached from the direction DIRN (either +1 or -1).
%
% See also CHEBCOLLOC.FEVAL

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

n = disc.dimension;

% Find the collocation points and create an empty functional.
[x, ~, v] = functionPoints(disc);
offset = cumsum([0 ; n(:)]);
N = offset(end);
E = zeros(1, N);

% Evaluate the given input at x:
y = y(x);
% Find which point is in which subinterval
intnum = chebfun.whichInterval(disc.domain, y);
intnum = max(intnum, 1); intnum = min(intnum,numel(n));
% Loop over each interval:
for i = unique(intnum)'
    j = i == intnum;
    active = offset(i) + (1:n(i));
    E(j,active) = barymat(y(j), x(active), v(active));
end

end


text:
function xy_display ( xy_filename )

%*****************************************************************************80
%
%% XY_DISPLAY plots a set of points in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Usage:
%
%    xy_display ( 'xy_filename' )
%
%  Parameters:
%
%    Input, string XY_FILENAME, the name of the file
%    containing the coordinates of the points.
%
  fprintf ( 1, '\n' );
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'XY_DISPLAY\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read a file of POINT_NUM points in 2 dimensions;\n' );
  fprintf ( 1, '  Display the points in a MATLAB graphics window.\n' );
%
%  First argument is the point file.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'XY_DISPLAY:\n' );
    xy_filename = input ( '  Enter the name of the point file (in ''quotes''!).\n' );
  end
%
%  Read the data.
%
  point_num = xy_header_read (  xy_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".\n', xy_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of points POINT_NUM  = %d\n', point_num );

  xy = xy_data_read ( xy_filename, point_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', xy_filename );

  r8mat_transpose_print_some ( 2, point_num, xy, 1, 1, ...
    2, 5, '  First 5 nodes:' );

  scatter ( xy(1,:), xy(2,:), 'filled' );

  xy_min(1) = min ( xy(1,:) );
  xy_max(1) = max ( xy(1,:) );

  xy_min(2) = min ( xy(2,:) );
  xy_max(2) = max ( xy(2,:) );

  xy_range(1) = xy_max(1) - xy_min(1);
  xy_range(2) = xy_max(2) - xy_min(2);

  margin = 0.025 * max ( xy_range(1), xy_range(2) );

  x_min = xy_min(1) - margin;
  x_max = xy_max(1) + margin;
  y_min = xy_min(2) - margin;
  y_max = xy_max(2) + margin;

  xlabel ( '--X axis--' )
  ylabel ( '--Y axis--' )
%
%  The TITLE function will interpret underscores in the title.
%  We need to unescape such escape sequences!
%
  title_string = s_escape_tex ( xy_filename );
  title ( title_string )

  axis ( [ x_min, x_max, y_min, y_max ] );
  axis equal

  fprintf ( 1, '\n' );
  fprintf ( 1, 'XY_DISPLAY\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, an optional title.
%
  incx = 5;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d       ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%12f', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end
function s2 = s_escape_tex ( s1 )

%*****************************************************************************80
%
%% S_ESCAPE_TEX de-escapes TeX escape sequences.
%
%  Discussion:
%
%    In particular, every occurrence of the characters '\', '_',
%    '^', '{' and '}' will be replaced by '\\', '\_', '\^',
%    '\{' and '\}'.  A TeX interpreter, on seeing these character
%    strings, is then likely to return the original characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, the string to be de-escaped.
%
%    Output, string S2, a copy of the string, modified to avoid TeX escapes.
%
  s1_length = length ( s1 );

  s1_pos = 0;
  s2_pos = 0;
  s2 = [];

  while ( s1_pos < s1_length )

    s1_pos = s1_pos + 1;

    if ( s1(s1_pos) == '\' || ...
         s1(s1_pos) == '_' || ...
         s1(s1_pos) == '^' || ...
         s1(s1_pos) == '{' || ...
         s1(s1_pos) == '}' )
      s2_pos = s2_pos + 1;
      s2 = strcat ( s2, '\' );
    end

    s2_pos = s2_pos + 1;
    s2 = strcat ( s2, s1(s1_pos) );

  end

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
%% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LEN, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
function xy = xy_data_read ( input_filename, point_num )

%*****************************************************************************80
%
%% XY_DATA_READ reads data from an XY file.
%
%  Discussion:
%
%    The number of points in the file can be determined by calling
%    XY_HEADER_READ first.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Output, real XY(2,POINT_NUM), the point coordinates.
%
  input_unit = fopen ( input_filename, 'r' );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'XY_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file.\n' );
    error ( 'XY_DATA_READ - Error!' );
  end

  i = 0;
  xy = zeros ( 2, point_num );
%
%  Use FGETL, not FGETS, so we don't pick up the NEWLINE character.
%
  while ( i < point_num )

    text = fgetl ( input_unit );

    if ( text(1) == '#' )
      continue;
    end

    if ( s_len_trim ( text ) == 0 )
      continue;
    end

    temp = sscanf ( text, '%f %f' );

    i = i + 1;
    xy(1,i) = temp(1);
    xy(2,i) = temp(2);

  end
%
%  Close the file.
%
  fclose ( input_unit );

  return
end
function point_num = xy_header_read ( input_filename )

%*****************************************************************************80
%
%% XY_HEADER_READ determines the number of pairs of data in an XY file.
%
%  Discussion:
%
%    This routine assumes that the file contains exactly three kinds of
%    records:
%
%    COMMENTS which begin with a '#' character in column 1;
%    BLANKS which contain nothing but 'whitespace';
%    XY coordinates, which each contain one pair of real values.
%
%    The routine ignores comments and blank lines and returns
%    the number of lines containing XY coordinates.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer POINT_NUM, the number of points in the file.
%
  point_num = 0;
%
%  Open the file.
%
  input_unit = fopen ( input_filename, 'r' );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'XY_HEADER_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file.\n' );
    error ( 'XY_HEADER_READ - Error!' );
  end
%
%  Read every line in the file.
%
%  Use FGETL, not FGETS, so we don't pick up the NEWLINE character.
%
  while ( 1 )

    text = fgetl ( input_unit );

    if ( text == -1 )
      break;
    end

    if ( text(1) == '#' )
      continue;
    end

    if ( s_len_trim ( text ) == 0 )
      continue;
    end

    sscanf ( text, '%f %f' );

    point_num = point_num + 1;

  end
%
%  Close the file.
%
  fclose ( input_unit );

  return
end


text:
function [infstack] = Intf_filt(infstack,SHP,phi_PL,Coh_cal,reference_ind,BroNumthre,Cohthre)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   This file is part of TomoSAR.
%
%   TomoSAR is distributed in the hope that it will be useful,
%   but without warranty of any kind; without even the implied
%   warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%   See the Apache License for more details.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Author : Dinh Ho Tong Minh (INRAE) and Yen Nhi Ngo, Jan. 2022
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if not(exist('BroNumthre', 'var'))
     BroNumthre=5;
end

if not(exist('Cohthre', 'var'))
     Cohthre=0.65;
end

[~,~,n_interf] = size(infstack);

phi_PL(:,:,reference_ind) = [];
mask_coh = Coh_cal > Cohthre;
mask_PS = SHP.BroNum>BroNumthre;  %PS keep

mask = and(mask_PS,mask_coh);
mask = repmat(mask,[1,1,n_interf]);

infstack(mask) = abs(infstack(mask)).*exp(1i*phi_PL(mask));

author: DinhHoTongMinh

repo: TomoSAR


text:
%% This function is used to get the endeffector force of the KUKA iiwa 7 R 800.
% function [ f ] = getEEF_Force( t )
% This function is used to get the endeffector force of the KUKA iiwa 7 R 800.
% t: is the TCP/IP connection
% f: is 1x3 cell array,  X,Y and Z components of the force acting on the endeffector
% Copy right, Mohammad SAFEEA, 3rd of May 2017

repo: KST-Kuka-Sunrise-Toolbox


text:
function [xi,w]=thirdOrderTetrahedronCubPoints()
%%THIRDORDERTETRAHEDRONCUBPOINTS Obtain third-order cubature points
%   for integration over a tetrahedron in 3D. The points and weights are
%   for the tetrahedron with vertices (1,0,0), (0,1,0), (0,0,1), and
%   (0,0,0), but can be transformed to any tetrahedron using
%   transformSimplexTetrahedronPts.
%
%INPUTS: None
%
%OUTPUTS: xi This is a 3XnumCubPoints set of points for the standard
%            tetrahedron.
%          w A 1XnumCubPoints set of cubature weights. This sums to the
%            volume of the standard tetrahedron (1/6).
%
%This function implements the points given in [1] (8 points).
%
%EXAMPLE:
%Given the vertices of the simplex, we compare a third-order moment
%computed using these cubature points to one computed using
%monomialIntSimplex. The results are the same within typical finite
%precision limits.
% [xi,w]=thirdOrderTetrahedronCubPoints();
% alpha=[1;1;1];
% theMoment=findMomentFromSamp(alpha,xi,w);
% intVal=monomialIntSimplex(alpha);
% RelErr=(theMoment-intVal)/intVal
%
%REFERENCES:
%[1] F. D. Witherden and P. E. Vincent, "On the identification of symmetric
%    quadrature rules for finite element methods," Computer and Mathematics
%    with Applications, vol. 69, no. 10, pp. 1232-1241, May 2015.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

M=[-0.34367339496723662642072827083693243093,  -0.34367339496723662642072827083693243093,   -0.9689798150982901207378151874892027072,   0.18162379004944980942342872025562069427;
   -0.34367339496723662642072827083693243093,   -0.9689798150982901207378151874892027072,  -0.34367339496723662642072827083693243093,   0.18162379004944980942342872025562069427;
    -0.9689798150982901207378151874892027072,  -0.34367339496723662642072827083693243093,  -0.34367339496723662642072827083693243093,   0.18162379004944980942342872025562069427;
   -0.34367339496723662642072827083693243093,  -0.34367339496723662642072827083693243093,  -0.34367339496723662642072827083693243093,   0.18162379004944980942342872025562069427;
   -0.78390550020314279176487322158837338344,  -0.78390550020314279176487322158837338344,   0.35171650060942837529461966476512015033,   0.15170954328388352390990461307771263906;
   -0.78390550020314279176487322158837338344,   0.35171650060942837529461966476512015033,  -0.78390550020314279176487322158837338344,   0.15170954328388352390990461307771263906;
    0.35171650060942837529461966476512015033,  -0.78390550020314279176487322158837338344,  -0.78390550020314279176487322158837338344,   0.15170954328388352390990461307771263906;
   -0.78390550020314279176487322158837338344,  -0.78390550020314279176487322158837338344,  -0.78390550020314279176487322158837338344,   0.15170954328388352390990461307771263906];
w=M(:,4);
xi=M(:,1:3)';
%Transform the points to the standard tetrahedron.
v1=[-1, 1,-1,-1;
    -1,-1,-1, 1;
    -1,-1, 1, -1];
v2=[1,0,0,0;
    0,1,0,0;
    0,0,1,0];
[A,d]=affineTransBetweenTetrahedra(v1,v2);
xi=bsxfun(@plus,A*xi,d);
w=w/8;
end

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


text:
function obj = trim_mask(obj)
% Exclude empty voxels from mask information in obj.volInfo structure, and re-make obj.volInfo
%
% ..
%    Tor Wager, 2013
% ..

obj = replace_empty(obj);
whomit = all(obj.dat == 0 | isnan(obj.dat), 2);

obj.dat(whomit, :) = [];
if ~isempty(obj.removed_voxels) && length(obj.removed_voxels) == obj.volInfo.n_inmask
    obj.removed_voxels(whomit) = [];
end
obj.volInfo.image_indx(obj.volInfo.wh_inmask(whomit)) = 0;
obj.volInfo.wh_inmask(whomit) = [];
obj.volInfo.xyzlist(whomit, :) = [];
obj.volInfo.cluster(whomit) = [];
obj.volInfo.n_inmask = length(obj.volInfo.wh_inmask);

if isa(obj, 'fmri_data')
    obj.mask.dat(whomit, :) = [];
    if ~isempty(obj.mask.removed_voxels) && length(obj.mask.removed_voxels) == obj.mask.volInfo.n_inmask
        obj.mask.removed_voxels(whomit) = [];
    end
    obj.mask.volInfo.image_indx(obj.mask.volInfo.wh_inmask(whomit)) = 0;
    obj.mask.volInfo.wh_inmask(whomit) = [];
    obj.mask.volInfo.xyzlist(whomit, :) = [];
    obj.mask.volInfo.cluster(whomit) = [];
    obj.mask.volInfo.n_inmask = length(obj.mask.volInfo.wh_inmask);
end

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function strat = beatTheBookie_timeseries(dat_dir, files, bet, marg, nValidOdds, hours, discard_games_randomly, prop2discard)

% loop through hours
% apply formula
% keep all bookies and time when the criteria is met
% place bet only for 1 or the three options for all the bookies
% available. Discard the other options

money = 0;
accuracy = [];
m = 2; % counter for the money
acc = 1; % counter for accuracy
st2 = 1;
st = 1;
z = 0;

for fi = 1 : length(files)

    if mod(fi,1000) == 0
        fprintf('Game #%d \n', fi);
    end

    fid = fopen([dat_dir files(fi).name], 'r');
    C = textscan(fid, repmat('%f ' , [1,72*3]), 'delimiter', ',');
    fclose(fid);

    data = cell2mat(C);

    % sanity check: games with no bets
    if sum(isnan(data(:))) == (32 * 216)
        fprintf('Game number %d, %s, has only nans \n', fi, files(fi).name);
        z = z + 1;
        continue;
    end

    % get result of the game
    name = strrep(files(fi).name, '.txt', '');
    C = strsplit(name, '_');
    score_home = str2double(C{end-1});
    score_away = str2double(C{end});

    if score_home > score_away
        result = 1;
    elseif score_home == score_away
        result = 2;
    else
        result = 3;
    end

    possible_earn = 0;

    for h = 1 : length(hours{1})

        avg_home = nanmean(data(:,hours{1}(h)));[max_home, bk_home_id] = nanmax(data(:, hours{1}(h)));
        avg_draw = nanmean(data(:,hours{2}(h)));[max_draw, bk_draw_id] = nanmax(data(:, hours{2}(h)));
        avg_away = nanmean(data(:,hours{3}(h)));[max_away, bk_away_id] =  nanmax(data(:, hours{3}(h)));

        averages = [avg_home avg_draw avg_away];
        maximums = [max_home max_draw max_away];
        bookies = [bk_home_id bk_draw_id bk_away_id];

        % check that there is at least one valid bet. If the bets are still not
        % online for this time, continue to the next game
        if isnan(avg_home) || isnan(avg_draw) || isnan(avg_away)
            continue
        end

        if isnan(max_home) && isnan(max_draw) && isnan(max_away)
            continue
        end

        % Apply formula and decide whether to bet
        earn_margin(1) = ((1 ./ avg_home - marg) * max_home - 1) * (sum(~isnan(data(:,hours{1}(h)))) > nValidOdds);
        earn_margin(2) = ((1 ./ avg_draw - marg) * max_draw - 1) * (sum(~isnan(data(:,hours{2}(h)))) > nValidOdds) ;
        earn_margin(3) = ((1 ./ avg_away - marg) * max_away - 1) * (sum(~isnan(data(:,hours{3}(h)))) > nValidOdds);

        if sum(earn_margin > 0) >= 1

            [~, id] = max(earn_margin);
            possible_earn = bet  * (maximums(id) - 1);
            break;
        end

    end

    % If match did not reach the criteria for betting, move to the next
    % match
    %if sum(bets_id(:)) == 0
    if possible_earn == 0
        continue
    end

    % Artibtrarily skip some games to simulate a real situation of betting
    % where we miss some of the games
    if discard_games_randomly
        num = rand(1);
        if num < prop2discard
            continue
        end
    end

    % Decide where to place the bet: home, draw or away
    if possible_earn > 0
        aux_possible_earn = possible_earn;
        bet_result = id;

        % There is an error in one of the odds: it was set to 126. I fix this here
        % I set a maximum earning of 15 for the bets
        if aux_possible_earn > 750
            continue
        end

        % calculate loss / earning
        if isequal(bet_result, result)
            money(m) = money(m-1) + aux_possible_earn;
            m = m + 1;
            accuracy(acc) = 1;
            ids(acc) = id;
            max_odds(acc) = maximums(id);
            mean_odds(acc) = averages(id);
            acc = acc + 1;
        else
            money(m) = money(m-1) - bet;
            m = m + 1;
            accuracy(acc) = 0;
            ids(acc) = id;
            max_odds(acc) = maximums(id);
            mean_odds(acc) = averages(id);
            acc = acc + 1;
        end

    end

    clear data bets_id possible_earn
end

strat.money = money;
strat.name = 'BeatTheBookies';
strat.accuracy = accuracy;
strat.ids = ids;
strat.max_odds = max_odds;
strat.mean_odds = mean_odds;

end

repo: BeatTheBookie


text:
function [projOrder, d, dRow, dVar] = getProjOrder(L)
%GETPROJORDER   Get projection order of a LINOP.
%   Each boundary and continuity constraint in a LINOP forces a reduction in the
%   total number of rows in the discrete operator, so that the composite is
%   square. The reduction is found by down-projection of the result of applying
%   the operator.
%
%   GETPROJORDER(DISC) returns a matrix of dimensions by which each column of
%   the system should be down-projected in order to end with a square system.

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

d = L.diffOrder;
dRow = max(d, [], 2);
dVar = max(d, [], 1);
projOrder = max(dVar, 0);

end


text:
%  Figure 4.16      Feedback Control of Dynamic Systems, 6e
%                   Franklin, Powell, Emami
%% script to generate Figure 4.16(a), 4.16 (b)
% case I, P control by reaction curve data.
% First, input the plant model and compute the reaction curve.
sysp=tf(1,[600 70 1])
set(sysp,'td', 5)

% From the curve, kp = 6.92 for proportional control
sysdP = tf(6.92,1)
% To get the closed loop response, we must use a pade approximation to the
% delay
syspade=pade(sysp,3);
sysforward=syspade*sysdP;
sysback=tf(1,1)
syscl=feedback(sysforward,sysback)
figure(1)
clf
step(syscl)
title('Figure 4.16a Closed loop response using Reaction Curve data')
gtext('P')
hold on
%case II, PI control
%kp = 6.22, TI=43.3
sysdPI=6.22*tf([43.3 1],[43.3 0])
sysforward=syspade*sysdPI;
syscl=feedback(sysforward,sysback)
step(syscl)
gtext('PI')
nicegrid
%
sysdP = tf(6.92/3,1)
% To get the closed loop response, we must use a pade approximation to the
% delay
syspade=pade(sysp,3);
sysforward=syspade*sysdP;
sysback=tf(1,1)
syscl=feedback(sysforward,sysback)
figure(2)
clf
step(syscl)
title('Figure 4.16b Closed loop response using Reaction Curve data')
gtext('P')
hold on
%case II, PI control
%kp = 6.22, TI=43.3
sysdPI=(6.22/3)*tf([43.3 1],[43.3 0])
sysforward=syspade*sysdPI;
syscl=feedback(sysforward,sysback)
step(syscl)
gtext('PI')
nicegrid;


text:
function lipidMask = op_CSICreateLipidMask(MRSIStruct)
    checkArguments(MRSIStruct);
    dimensionIndexes = checkDimensions(MRSIStruct);

end

function checkArguments(MRSIStruct)
    checkSpatialFT(MRSIStuct);
    checkSpectralFT(MRSIStuct);
end

function [dimensionIndexs] = checkDimensions(MRSIStruct)
    dataSize = getSize(MRSIStruct);
    dims = getDimension(MRSIStruct);

end

author: CIC-methods

repo: FID-A

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


text:
function v = load_sequence(inpat, from, to)

% load first frame to get dimensions, channels
im = single(imread(sprintf(inpat,from)));
sz = [size(im,1) size(im,2)];
ch =  size(im,3);
nf = to - from + 1;

% allocate video volume
v = single(zeros(sz(1), sz(2), ch, nf));

% load sequence
disp(sprintf('Reading %s', inpat))
for i = [from:to],

	filename = sprintf(inpat,i);
%	disp(sprintf('Reading %s', filename))
	im = single(imread(filename));
	v(:,:,:,i - from + 1) = im;

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function SF = SF_evaluation(MF)
RF = diff(MF,1,1);
RF1 = sqrt(mean(mean(RF.^2)));
CF = diff(MF,1,2);
CF1 = sqrt(mean(mean(CF.^2)));
SF = sqrt(RF1^2+CF1^2);
end

author: Linfeng-Tang

repo: Image-Fusion


text:
function [iclust, lambda, lambda0] = getIclust(stat, cl)

iclust = zeros(cl.Ly, cl.Lx);
lambda = zeros(cl.Ly, cl.Lx);
lambda0 = zeros(cl.Ly, cl.Lx);

for j = 1:numel(stat)
   ipix = stat(j).ipix(:);

   inew = stat(j).lam(:)>lambda(ipix) + 1e-6;
%    lambda(ipix(inew)) = stat(j).lam(inew);
%    lambda0(ipix(inew)) = stat(j).lambda(inew);

   lambda(ipix(inew)) = stat(j).lam(inew);
   lambda0(ipix(inew)) = stat(j).lam(inew);

   iclust(ipix(inew)) = j;
end

author: cortex-lab


text:
function [dv, deltaV1, deltaV2, deltaV1R, deltaV2R, xfrOrbit, deltaV1NTW, deltaV2NTW] = twoBurnOrbitChangeObjFunc(x, iniOrbit, finOrbit, gmuXfr)
%twoBurnOrbitChangeObjFunc Summary of this function goes here
%   Detailed explanation goes here
    burn1TA = AngleZero2Pi(x(1));
    burn2TA = AngleZero2Pi(x(2));
    xfrArcTOF = x(3);

%     if(iniOrbit(2)<1E-5 && finOrbit(2)<1E-5)
%         if(finOrbit(2)<1E-5)
%             finOrbit(2) = finOrbit(2) + 1E-5;
%         else
%             iniOrbit(2) = iniOrbit(2) + 1E-5;
%         end
%     end

%     if(iniOrbit(3)<1E-3 && finOrbit(3)<1E-3)
%         finOrbit(3) = finOrbit(3) + 1E-3;
%     end

    if(length(iniOrbit) >= 8)
        gmuXfr1 = iniOrbit(8);
    else
        gmuXfr1 = gmuXfr;
    end

    if(length(finOrbit) >= 8)
        gmuXfr2 = finOrbit(8);
    else
        gmuXfr2 = gmuXfr;
    end

    [rVect1,vVect1Pre]=getStatefromKepler(iniOrbit(1), iniOrbit(2), iniOrbit(3), iniOrbit(4), iniOrbit(5), burn1TA, gmuXfr1);
    [rVect2,vVect2Post]=getStatefromKepler(finOrbit(1), finOrbit(2), finOrbit(3), finOrbit(4), finOrbit(5), burn2TA, gmuXfr2);

    xfrArcTOFDays = xfrArcTOF/(86400);
    [vVect1PostSW,vVect2PreSW]=orbit.lambert(rVect1', rVect2', xfrArcTOFDays, 0, gmuXfr);
    [vVect1PostLW,vVect2PreLW]=orbit.lambert(rVect1', rVect2', -xfrArcTOFDays, 0, gmuXfr);

    deltaV1SW = vVect1PostSW' - vVect1Pre;
    deltaV1LW = vVect1PostLW' - vVect1Pre;
    deltaV2SW = vVect2Post - vVect2PreSW';
    deltaV2LW = vVect2Post - vVect2PreLW';

    dvSW = norm(deltaV1SW) + norm(deltaV2SW);
    dvLW = norm(deltaV1LW) + norm(deltaV2LW);
    if(dvSW < dvLW)
        dv = dvSW;
        deltaV1 = deltaV1SW;
        deltaV2 = deltaV2SW;
        [sma, ecc, inc, longAscNode, ArgPeri, TA1] = getKeplerFromState(rVect1,vVect1PostSW,gmuXfr);
        [~, ~, ~, ~, ~, TA2] = getKeplerFromState(rVect2,vVect2PreSW,gmuXfr);
        xfrOrbit = [sma, ecc, inc, longAscNode, ArgPeri, TA1, TA2];
        if(ecc<1.0)
            xfrOrbit(6) = AngleZero2Pi(TA1);
            xfrOrbit(7) = AngleZero2Pi(TA2);
        end

        deltaV1NTW = getNTWdvVect(deltaV1, rVect1, vVect1Pre);
        deltaV2NTW = getNTWdvVect(deltaV2, rVect2, vVect2PreSW);
    else
        dv = dvLW;
        deltaV1 = deltaV1LW;
        deltaV2 = deltaV2LW;
        [sma, ecc, inc, longAscNode, ArgPeri, TA1] = getKeplerFromState(rVect1,vVect1PostLW,gmuXfr);
        [~, ~, ~, ~, ~, TA2] = getKeplerFromState(rVect2,vVect2PreLW,gmuXfr);
        xfrOrbit = [sma, ecc, inc, longAscNode, ArgPeri, TA1, TA2];
        if(ecc<1.0)
            xfrOrbit(6) = AngleZero2Pi(TA1);
            xfrOrbit(7) = AngleZero2Pi(TA2);
        end

        deltaV1NTW = getNTWdvVect(deltaV1, rVect1, vVect1Pre);
        deltaV2NTW = getNTWdvVect(deltaV2, rVect2, vVect2PreLW);
    end
    deltaV1R = rVect1;
    deltaV2R = rVect2;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [fet clu spktimes wav]= ConvertKlusta2Matlab(shank,basepath,...
    basename,wvformExtract,saveFiles,numCores,keepMUA)
% USAGE
%
% [fet clu spktimes wav]= ConvertKlusta2Matlab(shank,basepath,basename,wvformExtract,saveFiles,numCores)
%
% INPUTS
%
% shank - the shank number (as a number, not text) equalling the name of
%         the folder under basepath with the data of interst.  Default = 1.
% basepath - directory path to the main recording folder with .dat and .xml
%            as well as shank folders made by makeProbeMapKlusta2.m (default is
%            current directory matlab is pointed to)
% basename - shared file name of .dat and .xml (default is last part of
%            current directory path, ie most immediate folder name)
% wvformExtract - binary (0 or 1) to choose whether to extract raw waveforms
% saveFiles - binary (0 or 1) to choose whether to save clu, res, fet, wav(.spk)
%             files or just return them as outputs
% numCores  - double between 1 and inf that sets the number of cores for a
%             parloop to run (set 2-12 on SSD or RAID0 systems, set to 1
%             for data stored on a single hard drive
% keepMUA   - binary (0 or 1) to choose whether to keep MUA as a separate
%             cluster (#1)
%
% OUTPUTS
%
% fet - features of waveforms
% clu - cluster ID's for waveforms
% spktimes(res) - spike times of waveforms(in sample #'s, not seconds)
% wav- waveform shapes taken from raw .dat file
%
%Converts .kwik/kwx files from Klusta into klusters-compatible
% fet,res,clu,spk files.  Works on a single shank of a recording, assumes a
% 16bit .dat and an .xml file is present in "basepath" (home folder) and
% that they are named basename.dat and basename.xml.  Also assumes that
% subdirectories in that basepath are made for each shank with names
% specified by numbers (ie 1,2,3,4..8).  In each shank folder should be
% .kwik and .kwx files made by klusta with names as follows:
% basename_sh[shankumber].kwik/kwx.
%
%
% Brendon Watson 2016
% edited by David Tingley 1/2017
% edited by Aza and Antonio 7/2017

% p = inputParser;
% addRequired(p,'shank',@isnumeric)
% addRequired(p,'basepath',@isstr)
% addRequired(p,'basename',@isstr)
%
% addParameter(p,'saveFiles',[],@isnumeric)
% addParameter(p,'numCores',[1],@isnumeric)
% p.parse(varargin{:})
%
% shank = p.Results.shank;
% basepath = p.Results.basepath;
% saveFiles = p.Results.saveFiles;
% numCores = p.Results.numCores;
% basename = p.Results.basename;

%
% Handle inputs
if ~exist('shank','var')
    shank = 1;
end
if ~exist('basepath','var');
    basepath = cd;
end
if ~exist('saveFiles','var');
    saveFiles = 0;
end
if ~exist('numCores','var');
    numCores = 1;
end
if ~exist('basename','var');
    [~,basename] = fileparts(pwd);
elseif strcmp(basename, 'lfpfile')
    d = dir('*lfp');
    basename = d.name(1:end-4);
end

if saveFiles % make file names
    cluname = fullfile(basepath,[basename '.clu.' num2str(shank)]);
    resname = fullfile(basepath,[basename '.res.' num2str(shank)]);
    fetname = fullfile(basepath,[basename '.fet.' num2str(shank)]);
    spkname = fullfile(basepath,[basename '.spk.' num2str(shank)]);
end

% Start grabbing data
datpath = fullfile(basepath,[basename '.dat']);
tkwik = fullfile(basepath,num2str(shank),[basename '_sh' num2str(shank) '.kwik']);
tkwx = fullfile(basepath,num2str(shank),[basename '_sh' num2str(shank) '.kwx']);
kwikinfo = h5info(tkwik,['/channel_groups/' num2str(shank) '/clusters/main']);
clu = h5read(tkwik,['/channel_groups/' num2str(shank) '/spikes/clusters/main']);
cluster_names = unique(clu);

try % different klusta versions have different names for this field
    totalch = h5readatt(tkwik,'/application_data/spikedetekt','nchannels');
catch
    totalch = h5readatt(tkwik,'/application_data/spikedetekt','n_channels');
end
try % some files don't have these variables saved in them?
    sbefore = h5readatt(tkwik,'/application_data/spikedetekt','extract_s_before');
    safter = h5readatt(tkwik,'/application_data/spikedetekt','extract_s_after');
catch
    wvform_nsamples = h5readatt(tkwik,'/application_data/spikedetekt','waveforms_nsamples');
    sbefore = wvform_nsamples/2;
    safter = wvform_nsamples/2;
end
channellist = h5readatt(tkwik,['/channel_groups/' num2str(shank)],'channel_order')+1;

%% Getting spiketimes
spktimes = h5read(tkwik,['/channel_groups/' num2str(shank) '/spikes/time_samples']);

%% spike extraction from dat
wav = [];
tsampsperwave = (sbefore+safter);
ngroupchans = length(channellist);
valsperwave = tsampsperwave * ngroupchans;
if wvformExtract & ~exist(spkname)  % only create .spk file if it doesn't already exist
    disp('loading from .dat file')
    dat=memmapfile(datpath,'Format','int16');
    wvforms_all = zeros(length(spktimes),valsperwave,'int16');

    if numCores>1 % if using multiple cores, start parfor loop
        delete(gcp('nocreate')); parpool(numCores)
        parfor j=1:length(spktimes)
            try
                byteList = [];
                for i = 1:length(channellist)
                    ind1 = (double(spktimes(j)).*totalch); % correction for problems with high values
                    ind11 = int64(ind1) - int64(sbefore.*totalch) + 1;
                    ind2 = (double(spktimes(j)).*totalch);
                    ind22 = int64(ind2) + int64(safter.*totalch);
                    ind = ind11:ind22;

                    byteList = [byteList; ind(double(channellist(i)):totalch:end)];
                end
                wvforms = dat.data(byteList);
            catch
                disp(['Error extracting spike at sample ' int2str(double(spktimes(j))) '. Saving as zeros']);
                disp(['Time range of that spike was: ' num2str(double(spktimes(j))-sbefore) ' to ' num2str(double(spktimes(j))+safter) ' samples'])
                wvforms = zeros(valsperwave,1,'int16');
            end
            if rem(j,50000) == 0
                disp([num2str(j) ' out of ' num2str(length(spktimes)) ' done'])
            end
            wvforms_all(j,:)=wvforms(:);
        end
    else % if we are using only one core, use a regular for loop
        for j=1:length(spktimes)
            try
                byteList = [];
                for i = 1:length(channellist)
                    ind1 = (double(spktimes(j)).*totalch); % correction for problems with high values
                    ind11 = int64(ind1) - int64(sbefore.*totalch) + 1;
                    ind2 = (double(spktimes(j)).*totalch);
                    ind22 = int64(ind2) + int64(safter.*totalch);
                    ind = ind11:ind22;

                    byteList = [byteList; ind(double(channellist(i)):totalch:end)];
                end
                wvforms = dat.data(byteList);
            catch
                disp(['Error extracting spike at sample ' int2str(double(spktimes(j))) '. Saving as zeros']);
                disp(['Time range of that spike was: ' num2str(double(spktimes(j))-sbefore) ' to ' num2str(double(spktimes(j))+safter) ' samples'])
                wvforms = zeros(valsperwave,1,'int16');
            end
            if rem(j,50000) == 0
                disp([num2str(j) ' out of ' num2str(length(spktimes)) ' done'])
            end
            wvforms_all(j,:)=wvforms(:);
        end

    end

    clear dat
    wvforms_all = reshape(wvforms_all',size(wvforms_all,1)*size(wvforms_all,2),1);
    wav = reshape(wvforms_all,ngroupchans,tsampsperwave,[]);
    wav = permute(wav,[3 1 2]);
elseif wvformExtract & exist(spkname) % load from .spk if it exists
    disp('loading from .spk file')
    fid = fopen(spkname,'r');
    wav = fread(fid,[1 inf],'int16');
    wav = reshape(wav,ngroupchans,tsampsperwave,[]);
    wav = permute(wav,[3 1 2]);
end
%% Spike features

fets = h5read(tkwx,['/channel_groups/' num2str(shank) '/features_masks']);
fets = double(squeeze(fets(1,:,:)));
fet = fets';
fetMultiplier = double(intmax('int32'))/max(abs(fets(:))); % masked klustakwik has small floats that need to be expanded before rounding below
fet = fet .* fetMultiplier;

%% writing to clu, res, fet, spk
if saveFiles

    %clu
    for i=1:length(kwikinfo.Groups)  % sometimes .kwik has more groups than clu (empty?)..
        group(i) = h5readatt(tkwik,kwikinfo.Groups(i).Name,'cluster_group');
        temp = strsplit(kwikinfo.Groups(i).Name,'/');
        clusterID(i) = str2num(temp{length(temp)}); % fix for if channel ordering is mixed up (we don't know that the first clusters are noise/MUA)
    end

    cluster_names = unique(clu);
    for i=1:length(cluster_names)
        [a b] =ismember(cluster_names(i),double(clusterID));
        if isempty(b)
            warning('could not find cluster name?')
        end
        if keepMUA
            if group(b) == 0 || group(b) == 3
                clu(clu == clusterID(b)) = 0; % noise
            elseif group(b) == 1
                clu(clu == clusterID(b)) = 1; % MUA
            end

        else
            if group(b)~=2
                clu(clu == clusterID(b)) = 0;  % re-name unsorted and noise as MUA/Noise cluster for FMATToolbox
            end
        end
    end

    clu = cat(1,length(unique(clu)),double(clu));
    fid=fopen(cluname,'w');
    fprintf(fid,'%.0f\n',clu);
    fclose(fid);
    clear fid

    fid=fopen(resname,'w');
    fprintf(fid,'%.0f\n',spktimes);
    fclose(fid);
    clear fid

    SaveFetIn(fetname,fet);
    %spk

    if wvformExtract & ~exist(spkname)
        fid=fopen(spkname,'w');
        fwrite(fid,wvforms_all,'int16');
        fclose(fid);
        clear fid
    end
    disp(['Shank ' num2str(shank) ' done'])
end

function SaveFetIn(FileName, Fet, BufSize)
if nargin<3 | isempty(BufSize)
    BufSize = inf;
end

nFeatures = size(Fet, 2);
formatstring = '%d';
for ii=2:nFeatures
    formatstring = [formatstring,'\t%d'];
end
formatstring = [formatstring,'\n'];

outputfile = fopen(FileName,'w');
fprintf(outputfile, '%d\n', nFeatures);

if isinf(BufSize)
    fprintf(outputfile,formatstring,round(Fet'));
else
    nBuf = floor(size(Fet,1)/BufSize)+1;

    for i=1:nBuf
        BufInd = [(i-1)*nBuf+1:min(i*nBuf,size(Fet,1))];
        fprintf(outputfile,formatstring,round(Fet(BufInd,:)'));
    end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function w = loadsac(w, filename)
   %LOADSAC  creates a waveform from a SAC file
   %
   % old usage, use waveform with a datasource instead
   %
   % waveform = loadsac(waveform, filename);
   % To call, the first argument must be a waveform.  However this waveform
   % will never be accessed-- it will be overwritten by the SAC information.
   % Therefore, you can use "waveform" as your first argument, so that a
   % generic waveform object will be created for you
   %
   % If the filename is a cell of several SAC file names, such as:
   % fn = {'BOB.SHZ.AV','TOM.EHZ.AV','HARRY.EHZ.AV'}
   % then, w = loadsac(waveform,fn) will return a matrix of waveforms, the
   % same size and shape  as fn.   fn(N,M) will be the SAC file name for the
   % waveform denoted by w(N,M).  In the example above, loadsac will return a
   % 1x3 waveform object, with all three files.
   %
   % Programming suggestion?select files using a dialog box, by using the
   % following code:
   %   % display the dialog box used for selecting files, and
   %   % allow for multiple selections
   %   [f,d] = uigetfile({'*.*','all Files (*.*)'}, 'MultiSelect', 'on');
   %
   %   % f will have one or more file names, so concatenate them
   %   % with the directory
   %   filename = strcat(d,f);
   %
   %   %call loadsac to load all these files into waveform objects
   %   w = loadsac(waveform, filename);
   %
   % The SAC header is parsed out, with the following equivelent fields used
   % by waveform:
   % _waveform_ 	_sac_
   % STATION		KSTNM
   % CHANNEL		KCMPNM
   % FREQUENCY		1 / DELTA (ODELTA is not used)
   % START			NZYEAR, NZJDAY, NZHOUR, NZMIN, NZSEC, NZMSEC, B
   % UNITS			IDEP is parsed.  IUNK becomes "Counts"
   % All header fields with values - that is, header fields that do not
   % contain the value  -12345 (either numeric or text) - are put into
   % user-defined fields within the waveform.
   % These fields can then be accessed via get/set.  User-defined fields
   % does not include those that are parsed into waveform's predefined  fields
   %
   % For more information about the SAC header, and the interpretation of the
   % fields, check out Lawrence Livermore's web site.  Search for "seismic
   % sac" and you'll find it.
   %
   % -------------------------------------------------------------------------
   % The functionality of loadsac has been redistributed to the datasource
   % class.  This will allow greater flexibility in traversing directory/file
   % structures and loading arbitrary time-ranges from the SAC files.
   % -------------------------------------------------------------------------
   %
   %  Here are examples of loading sac files using the datasource class
   %  sacfile of interest: "2008/myfile.MSOM.035.sac"
   %
   % % here is a datasource that finds files such as:"2008/myfile.MSOM.035.sac"
   % ds = datasource('sac',...        % type of data
   %    '%04d/myfile.%s.%03d.sac',... % filename described in sprintf terms
   %    'year','station','jday');     % used to determine filename at runtime
   %
   % %  and here is a datasource tied to a specific file
   % ds2 = datasource('sac','specificpath/specificfilename');
   %
   % %     Here's a scnlobject filled with wildcards,
   % %     that will grab any station or channel
   % anyscnl = scnlobject('*','*','*','*');
   %
   % %     and here's one that is specific to our station of interest
   % scnl = scnlobject('MSOM','BHZ','AV'); %% specific station
   %
   % % load any station for a timespan of interest
   % w = waveform(ds,anyscnl,'2/4/2008 04:40','2/4/2008 04:53');
   %    or
   % % load our specific sacfile, making sure to grab the whole thing
   % w = waveform(ds2,anyscnl,'1/1/2000','12/31/2009');

   % VERSION: 1.1 of waveform objects
   % AUTHOR: Celso Reyes
   % LASTUPDATE: 3/14/2009

   w = load_sac(filename);
end

author: geoscience-community-codes

repo: GISMO