text:
function xyf_display ( prefix )

%*****************************************************************************80
%
%% XYF_DISPLAY displays points and faces in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Usage:
%
%    xyf_display ( 'prefix' )
%
%    where
%
%    'prefix.xy' contains the point coordinates
%    'prefix.xyf' contains lists of point indices forming faces.
%
%  Parameters:
%
%    Input, string PREFIX, the prefix of the two input files.
%
  fprintf ( 1, '\n' );
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'XYF_DISPLAY\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read an XY file containing coordinates of points in 2D,\n' );
  fprintf ( 1, '  and an XYF file listing point indices forming faces.\n' );
  fprintf ( 1, '  Display the points and faces in a MATLAB graphics window.\n' );
%
%  First argument is the file prefix.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'XYF_DISPLAY:\n' );
    prefix = input ( '  Enter the common prefix of the XY and XYF files (in ''quotes''!).\n' );
  end
%
%  Read the XY data.
%
  xy_filename = strcat ( prefix, '.xy' );

  point_num = xy_header_read ( xy_filename );

  xy_header_print ( point_num );

  xy = xy_data_read ( xy_filename, point_num );
%
%  Read the XYF data.
%
  xyf_filename = strcat ( prefix, '.xyf' );

  [ face_num, face_data_num ] = xyf_header_read ( xyf_filename );

  xyf_header_print ( point_num, face_num, face_data_num );

  [ face_pointer, face_data ] = xyf_data_read ( xyf_filename, ...
    face_num, face_data_num );

  if ( 0 )

    xyf_data_print ( point_num, face_num, face_data_num, face_pointer, ...
      face_data );

  end
%
%  If necessary, move from 0-based indexing.
%
  if ( min ( face_data(1:face_data_num) ) == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Moving FACE_DATA from 0 base to 1 base.\n' );
    face_data(1:face_data_num) = face_data(1:face_data_num) + 1;
  end
%
%  Get the scale.
%
  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:2) = xy_max(1:2) - xy_min(1: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;
%
%  Draw the picture.
%
  face_color = 'g';

  for face = 1 : face_num
    index = face_pointer(face) : face_pointer(face+1) - 1;
    p = face_data(index);
    patch ( xy(1,p), xy(2,p), face_color );
  end

  line_color = 'r';

  for face = 1 : face_num
    index = [face_pointer(face) : face_pointer(face+1) - 1, face_pointer(face)];
    p = face_data(index);
    line ( 'XData', xy(1,p), 'YData', xy(2,p), 'Color', line_color );
  end
%
%  May I take a moment hear to complain that MATLAB's graphic capabilities
%  are awkward and poorly documented?  You can draw a line, if you work hard
%  at it.  But can you draw a point?  Nope!
%
  hold on

  point_size = 40;
  point_color = 'b';
  scatter ( xy(1,:), xy(2,:), point_size, point_color, 'filled' );

  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 ( prefix );
  title ( title_string )

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

  hold off
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'XYF_DISPLAY\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

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

  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 word_num = s_word_count ( s )

%*****************************************************************************80
%
%% S_WORD_COUNT counts the number of "words" in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be examined.
%
%    Output, integer WORD_NUM, the number of "words" in the string.
%    Words are presumed to be separated by one or more blanks.
%
  FALSE = 0;
  TRUE = 1;

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = TRUE;
    elseif ( blank == TRUE )
      word_num = word_num + 1;
      blank = FALSE;
    end

  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 "%s".\n', input_filename );
    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 xy_header_print ( point_num )

%*****************************************************************************80
%
%% XY_HEADER_PRINT prints the header of an XY file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer POINT_NUM, the number of points.
%
  fprintf ( 1, '\n');
  fprintf ( 1, '  Number of points = %d\n', point_num );

  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 "%s".\n', input_filename );
    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
function [ face_pointer, face_data ] = xyf_data_read ( input_filename, ...
  face_num, face_data_num )

%*****************************************************************************80
%
%% XYF_DATA_READ reads the data in an XYF 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';
%    FACE ITEMS, which are indices of points on a face, or -1 to terminate a face.
%
%    The routine ignores comments and blanks and returns
%    the number of face items.
%
%  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 FACE_NUM, the number of faces.
%
%    Input, integer FACE_DATA_NUM, the number of face items.
%
%    Output, integer FACE_POINTER(FACE_NUM+1), pointers to the
%    first face item for each face.
%
%    Output, integer FACE_DATA(FACE_DATA_NUM), the face items.
%
  input_unit = fopen ( input_filename, 'r' );

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

  num = 0;
  face_pointer = zeros ( face_num + 1, 1 );
  face_data = zeros ( face_data_num, 1 );
  face_pointer(1) = 1;
%
%  Use FGETL, not FGETS, so we don't pick up the NEWLINE character.
%
  while ( num < face_num )

    text = fgetl ( input_unit );

    if ( text(1) == '#' || s_len_trim ( text ) == 0 )
      continue
    end

    num = num + 1;

    n = s_word_count ( text );
    face_pointer(num+1) = face_pointer(num) + n;

    ilo = face_pointer(num);
    ihi = face_pointer(num+1) - 1;

    [ temp, count ] = sscanf ( text, '%d' );

    face_data(ilo:ihi) = temp(1:count);

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

  return
end
function xyf_header_print ( point_num, face_num, face_data_num )

%*****************************************************************************80
%
%% XYF_HEADER_PRINT prints the header of an XYF file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer FACE_NUM, the number of faces.
%
%    Input, integer FACE_DATA_NUM, the number of face items.
%
  fprintf ( 1, '\n');
  fprintf ( 1, '  Number of points =     %d\n', point_num );
  fprintf ( 1, '  Number of faces =      %d\n', face_num );
  fprintf ( 1, '  Number of face items = %d\n', face_data_num );

  return
end
function [ face_num, face_data_num ] = xyf_header_read ( input_filename )

%*****************************************************************************80
%
%% XYF_HEADER_READ determines the number of face items in an XYF 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';
%    FACE ITEMS, which are indices of points on a face, or -1 to terminate a face.
%
%    The routine ignores comments and blanks and returns
%    the number of face items.
%
%  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 FACE_NUM, the number of faces.
%
%    Output, integer FACE_DATA_NUM, the number of face items.
%
  face_data_num = 0;
  face_num = 0;
%
%  Open the file.
%
  input_unit = fopen ( input_filename, 'r' );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'XYF_HEADER_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file "%s".\n', input_filename );
    error ( 'XYF_HEADER_READ - Error!' );
  end
%
%  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) == '#' || s_len_trim ( text ) == 0 )
      continue
    end

    n = s_word_count ( text );

    face_data_num = face_data_num + n;

    face_num = face_num + 1;

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

  return
end


text: prior = GaussianStateX([0; 0.04979784; 0;-5], diag([10, 1, 10, 1]));

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function   IM=subplane(Z)
%SUBtract selected PLANE, or remove slope (detilt) from data matrix
%
%Call:
%           IM=subplane(Z)
%Input:
%           Z = (double) data matrix
%Output:
%           IM = Z - plane     (the SUBtracted PLANE is defined by manually selected points)
%
%Vassili Pastushenko	March	2005
%==============================
%figure;
imagesc(Z);
shg
set(gca,'fontsize',15)
%select at least three (x,y) points which define a plane in 3D
title('Click at least three coplanar points')
[x,y]=getpts(gca);

if numel(x)<3,
    title('More points please');
    [xx,yy]=getpts(gca);
    x=[x;xx];y=[y;yy];
end

x=round(x);
y=round(y);
PONT=numel(x);
ZIN=ones(PONT,1);
M=[x y ZIN];
for i=1:PONT
ZIN(i)=Z(y(i),x(i));
end

V=M\ZIN;
[VY,VX]=size(Z);
[X,Y]=meshgrid(1:VX,1:VY);
BAS=V(1)*X+V(2)*Y+V(3);
IM=Z-BAS;
imagesc(IM);
title('Detilted data')
shg

author: Sable

repo: mcbench-benchmarks


text:
function slide21

	figure('units','normalized',...
		'position',[0.1, 0.1, 0.7, 0.7],...
		'menubar','none');

	axes('units','normalized',...
		'position',[0.5, 0.5, 0.4, 0.4])

	x1 = 3;
	x2 = 5;
	y1 = 2;
	y2 = 3;

	patch([x1,x1,x2,x2],[y1,5,5,y1],[1,1,1]); % tank patch
	p1 = patch([x1,x1,x2,x2],[y1,y2,y2,y1],[0,0,1]); % liquid patch

	axis([0 6 0 6])

	t = 1;
	while true

		t = t + 0.05;

		% the upper left corner plots a cosine while the upper right corner
		% plots a sine function
		y = [y1, y1 + 2 + 0.25*cos(t),  y1 + 2 + 0.25*sin(t), y1];

		set(p1,'ydata',y)

		drawnow

	end

end


text:
function results = SOFAgetVersion(flag)
%SOFAgetVersion - Get the version of the SOFA Toolbox
%   Usage: ver = SOFAgetVersion(flag)
%
%   ver = SOFAgetVersion() returns the version of the SOFA Toolbox
%
%   ver = SOFAgetVersion('Toolbox') and SOFAgetVersion('API') do the same (not case-sensitive).
%
%   ver = SOFAgetVersion('SOFA') returns the version of the SOFA supported
%   by the SOFA Toolbox.

% #Author: Piotr Majdak
% #Author: Michael Mihocic: header documentation updated (28.10.2021)
%
% SOFA Toolbox - function SOFAgetVersion
% Copyright (C) Acoustics Research Institute - Austrian Academy of Sciences;
% Licensed under the EUPL, Version 1.2 or - as soon they will be approved by the European Commission - subsequent versions of the EUPL (the "License")
% You may not use this work except in compliance with the License.
% You may obtain a copy of the License at: https://joinup.ec.europa.eu/software/page/eupl
% Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing  permissions and limitations under the License.

if ~exist('flag','var')
	flag='TOOLBOX';
end

switch upper(flag)
	case {'TOOLBOX','API'}
		results = '2.1.5';
	case 'SOFA'
		results = '2.1';
end

author: sofacoustics

repo: SOFAtoolbox


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal, Marc Bakry (c) 2017-2019.                 |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%|             marc.bakry@polytechnique.edu                               |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtHmxVibroSlab2d.m                           |
%|    #    |   VERSION    : 0.55                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal & Marc Bakry                  |
%|  ( # )  |   CREATION   : 14.03.2019                                    |
%|  / 0 \  |   LAST MODIF :                                               |
%| ( === ) |   SYNOPSIS   :                                               |
%|  `---'  |                                                              |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Accuracy
tol = 1e-3

% Width of the slab
e = 0.5

% Frequency
f = 900:500:4000

% Incident direction (from bottom)
X0  = [0 1 0];

% Exterior domain (water)
rho0 = 1000;         % density      (kg.m3)
c0   = 1500;         % celerity     (m.s-1)
k0   = 2*pi/c0.*f;   % wave-number  (m-1)
lam0 = c0./f;        % wave-length  (m)

% Interior domain (different from water)
rhoS = 2*rho0;                                 % density                             (kg.m3)
cL   = 2*c0;                                   % celerity of longitudinal waves      (m.s-1)
kL   = 2*pi.*f./cL;                            % wave-number of longitudinal waves   (m-1)
lamL = real(cL)./f;                            % wavelength of longitudinal waves    (m)
cT   = 0;                                      % celerity of transverse waves        (m.s-1)
kT   = 2*pi.*f./cT;                            % wave-number of transverse waves     (m-1)
lamT = real(cT)./f;                            % wavelength of transverse waves      (m)

% Solution (pressure)
sol = zeros(2,length(f));

% Loop for each frequency
for i = 1:length(f)
    % Minimum wavelength
    tmp  = [lam0(i),lamL(i),lamT(i)];
    lmin = min(tmp(tmp>0));

    % Slab mesh
    L    = 60 * lmin;            % 60 wavelength to simulate infinite slab
    nx   = ceil(L/lmin * 6)+1;   % 6 node per wavelength for L
    ny   = ceil(e/lmin * 12)+1;  % 12 node per wavelength for e
    N    = nx * ny;              % Total number of nodes
    mesh = mshSquare(N,[L e])

    % Radiative mesh (fixed number of nodes)
    radiat = mshSquare(1e3,[L L]);

    % Boundary
    bound = swap(mesh.bnd)

    % Measurement points for trans and refl coeff (1 wavelenth from bound)
    Xmes = [0 -e/2-lmin 0 ; 0 e/2+lmin 0];

    % Cut-off function (50% full, 10% decrease)
    cutoff = vibsCutoff(1,L/5,L/10);

    % Green kernel function
    Gxy         = @(X,Y) femGreenKernel(X,Y,'[H0(kr)]',k0(i));
    gradyGxy{1} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]1',k0(i));
    gradyGxy{2} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]2',k0(i));
    gradyGxy{3} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]3',k0(i));

    % Plane wave function
    PW         = @(X) exp(1i*k0(i)*X*X0') .* cutoff(X);
    gradxPW{1} = @(X) 1i*k0(i)*X0(1) .* PW(X);
    gradxPW{2} = @(X) 1i*k0(i)*X0(2) .* PW(X);
    gradxPW{3} = @(X) 1i*k0(i)*X0(3) .* PW(X);

    % Coupling coeff for Brackage-Werner simulation
    beta = 1i*k0(i);

    % Quadrature and finite elements (volumn)
    omega = dom(mesh,3);
    U     = fem(mesh,'P1');

    % Quadrature and finite elements (boundary)
    sigma = dom(bound,3);
    u     = fem(bound,'P1');

    % Left-hand side
    tic
    [A,B,C,D] = vibsHmxBlockOperator(omega,U,sigma,u,cL,cT,rhoS,c0,rho0,f(i),tol);
    toc

    % Add dirichlet condition to x unknows (penalization)
    A(sub2ind(size(A),1:length(U),1:length(U))) = 1e15;

    % Right-hand side
    V    = cell(3,1);
    V{1} = - integral(sigma,ntimes(U,1),PW);
    V{2} = - integral(sigma,ntimes(U,2),PW);
    V{3} = integral(sigma,ntimes(u),gradxPW);

    % Resolution with Schur complement
    Fa     = decomposition(A);
    LHS    = @(V) D*V - C*(Fa \ (B.Ml*(B.Mr*V)) );
    RHS    = V{end} - C*(Fa \ cell2mat(V(1:end-1)) );
    mu     = mgcr(LHS,RHS,[],tol,100);
    lambda = beta*mu;

    % Measure of refexive and transmitted coeff
    tic
    Pmes = 1i/4 .* integral(Xmes,sigma,Gxy,u)*lambda - ...
        1i/4 .* integral(Xmes,sigma,gradyGxy,ntimes(u))*mu;
    Pmes(2) = Pmes(2) + PW(Xmes(2,:));
    toc

    % Save solution
    sol(:,i) = Pmes;
end

% Analytical solution
tic
c     = ones(length(f),1) * [c0 cL c0];
rho   = [rho0 rhoS rho0];
[R,T] = slabVibro(f,rho,c,e);
toc

% Comparison (db)
ref = 20*log10(abs([R ; T]));
sol = 20*log10(abs(sol));

norm(ref-sol)/norm(ref)

% Graphical representation
figure(100)
subplot(1,2,1)
plot(f,ref(1,:),'r',f,sol(1,:),'b+')
grid on
title('Reflexion coeffiscient')
legend({'Analytical','Numerical'})
xlabel('Frequency (Hz)')
ylabel('Amplitude (dB)')

subplot(1,2,2)
plot(f,ref(2,:),'r',f,sol(2,:),'b+')
grid on
title('Transmission coeffiscient')
legend({'Analytical','Numerical'})
xlabel('Frequency (Hz)')
ylabel('Amplitude (dB)')

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function ChannelMat = in_channel_besa_eps(ChannelFile)
% IN_CHANNEL_BESA_EPS:  Read 3D cartesian positions for a set of electrodes from an couple .eps/ela BESA file.
%
% USAGE:  ChannelMat = in_channel_besa_eps(ChannelFile)
%
% INPUTS:
%     - ChannelFile : Full path to the file (either .ela or .eps file)
%
% FORMAT:
%     - .ELA file: contains the labels
%     - .EPS file: contains the positions of the electrodes

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2009-2012

% Check file extension
[fPath, fBase, fExt] = bst_fileparts(ChannelFile);
switch lower(fExt)
    case '.eps'
        epsFile = ChannelFile;
        elaFile = bst_fullfile(fPath, [fBase, '.ela']);
    case '.ela'
        elaFile = ChannelFile;
        epsFile = bst_fullfile(fPath, [fBase, '.eps']);
end

% ELP: Get the positions
if ~isempty(epsFile) && file_exist(epsFile)
    ChannelMat = in_channel_ascii(epsFile, {'-Y','X','Z'}, 0, .01);
    ChannelMat.Comment = 'BESA channels';
else
    disp('BESA> Error: EPS file is missing');
    ChannelMat = [];
end

% ELA: Get the labels/types
if ~isempty(elaFile) && file_exist(elaFile)
    % === READ FILE ===
    % Open file
    fid = fopen(elaFile, 'r');
    if (fid == -1)
        error('Cannot open .ela file.');
    end
    % Read file
    ela = textscan(fid, '%s %s');
    % Close file
    fclose(fid);

    % === FILE TYPE ===
    nChan = length(ela{1});
    % ONE column (channel names)
    if all(cellfun(@isempty, ela{2}))
        chLabels = ela{1};
        % By default, all channels are considered as EEG
        chTypes  = repmat({'EEG'}, [nChan, 1]);
    % TWO columns (channel types+names)
    else
        chTypes  = ela{1};
        chLabels = ela{2};
    end

    % === SET THE LABELS ===
    % Create an empty Channel structure if EPS file is missing
    if isempty(ChannelMat)
        ChannelMat = db_template('channelmat');
        ChannelMat.Comment = 'BESA channels';
        ChannelMat.Channel = repmat(db_template('channeldesc'), [1, nChan]);
        [ChannelMat.Channel.Loc] = deal([0;0;0]);
    end
    % Set the labels
    for i = 1:length(chLabels)
        ChannelMat.Channel(i).Name = chLabels{i};
    end
    % Set the types
    for i = 1:length(chLabels)
        ChannelMat.Channel(i).Type = chTypes{i};
    end
else
    disp('BESA> Error: ELA file is missing');
end

author: brainstorm-tools


text:
% ------------------------------------------------------
% SwarmOps - Heuristic optimization for Matlab
% Copyright (C) 2003-2010 Magnus Erik Hvass Pedersen.
% Please see the file license.txt for license details.
% SwarmOps on the internet: http://www.Hvass-Labs.org/
% ------------------------------------------------------

% Behavioural parameters for Local Unimodal sampling
% from (1). The parameter-array consists of
% the following parameter:
% - gamma
%
% Literature references:
% (1) M.E.H. Pedersen. Tuning & Simplifying Heuristical
%     Optimization. PhD Thesis, University of Southampton,
%     2010.

LUS_DEFAULT = [3.0];

% ------------------------------------------------------


text:
function a = bimarkov_random ( n, key )

%*****************************************************************************80
%
%% BIMARKOV_RANDOM returns the BIMARKOV_RANDOM matrix.
%
%  Discussion:
%
%    This is a random biMarkov or doubly stochastic matrix.
%
%  Example:
%
%    N = 5
%
%    1/5   1/5   1/5   1/5   1/5
%    1/2   1/2    0     0     0
%    1/6   1/6   2/3    0     0
%    1/12  1/12  1/12  3/4    0
%    1/20  1/20  1/20  1/20  4/5
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    0 <= A(I,J) <= 1.0 for every I and J.
%
%    A has constant row sum 1.
%
%    A has constant column sum 1.
%
%    All the eigenvalues of A have modulus 1.
%
%    1 is always an eigenvalue of A, with eigenvector (1,1,...,1).
%
%    The eigenvalue 1 lies on the boundary of all the Gershgorin
%    row or column sum disks.
%
%    Every doubly stochastic matrix is a combination
%      A = w1 * P1 + w2 * P2 + ... + wk * Pk
%    of permutation matrices, with positive weights w that sum to 1.
%    (Birkhoff's theorem, see Horn and Johnson.)
%
%    A is a Markov matrix.
%
%    A is a transition matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Roger Horn, Charles Johnson,
%    Matrix Analysis,
%    Cambridge, 1985.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, integer KEY, a positive value that selects the data.
%
%    Output, real A(N,N), the matrix.
%

%
%  Get a random orthogonal matrix.
%
  a = orth_random ( n, key );
%
%  Square each entry.
%
  a(1:n,1:n) = a(1:n,1:n).^2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% [base,ext] = splitext(name)
function [base,ext] = splitext(name)

k = strfind(name,'.');
if isempty(k),
  base = name;
  ext = '';
  return;
end

k = k(end);
j = strfind(name,'/');
if ~isempty(j) && j(end) > k,
  base = name;
  ext = '';
  return;
end

base = name(1:k-1);
ext = name(k:end);

author: kristinbranson

repo: JAABA


text:
% Initialize path
M_ROOT = fileparts(which(mfilename));
addpath([M_ROOT '/vgg']);
addpath([M_ROOT '/ojw']);
clear M_ROOT

repo: SelfExSR


text:
function [z, mu, sigma] = VBA_zscore (x, flag, dim, omitnan)
% // VBA toolbox //////////////////////////////////////////////////////////
%
% [z, mu, sigma] = VBA_zscore (x, flag, dim, ignore_nans)
% standard score, aka z-score
%
% Compute the stardardized zscore of x, ie z = (x - mean(x)) / std(x)
%
% IN:
%   x: vector or matrix to normalize
%   flag: if true, normalizes x using std(X,1), i.e., by computing the
%     standard deviation(s) using N rather than N-1, where N is the length of
%     the dimension along which VBA_zscore works. VBA_zscore(x,0) is the same as
%     VBA_zscore(x). Default = false;
%   dim: dimension along which the standardization is performed. By
%     default, uses the first non singleton dimension.
%   ignore_nans: if true, use nan-safe normalization (default = false).
%
% OUT: z-score of x
%
% /////////////////////////////////////////////////////////////////////////

if nargin < 2 || isempty (flag)
	flag = 0;
else
    flag = +(flag > 0);
end

if nargin < 3 || isempty (dim)
	dim = find (size(x) > 1, 1);
	if isempty(dim)
		dim = 1;
	end
end

if nargin < 4
    omitnan = false;
end

if omitnan
    mu = VBA_nanmean (x, dim);
    sigma = sqrt(VBA_nanvar (x, dim, flag));
else
    mu = mean (x, dim);
    sigma = std (x, flag, dim);
end
sigma(sigma == 0) = 1;

z = bsxfun(@minus, x, mu);
z = bsxfun(@rdivide, z, sigma);

author: MBB-team

repo: VBA-toolbox


text:
function gX = whitefixedKernGradX(kern, X, X2)

% WHITEFIXEDKERNGRADX Gradient of WHITEFIXED kernel with respect to a point x.
% FORMAT
% DESC computes the gradient of the fixed parameter white noise
% kernel with respect to the input positions.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG x : locations against which gradients are being computed.
% RETURN g : the returned gradients. The gradients are returned in
% a matrix which is numData x numInputs x numData. Where numData is
% the number of data points and numInputs is the number of input
% dimensions in X.
%
% FORMAT
% DESC computes the gradident of the fixed parameter white noise
% kernel with respect to the input positions where both the row
% positions and column positions are provided separately.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG x1 : row locations against which gradients are being computed.
% ARG x2 : column locations against which gradients are being computed.
% RETURN g : the returned gradients. The gradients are returned in
% a matrix which is numData2 x numInputs x numData1. Where numData1 is
% the number of data points in X1, numData2 is the number of data
% points in X2 and numInputs is the number of input
% dimensions in X.
%
% SEEALSO whitefixedKernParamInit, kernGradX, whitefixedKernDiagGradX
%
% COPYRIGHT : Nathaniel J. King, 2006

% KERN

gX = zeros(size(X2, 1), size(X2, 2), size(X, 1));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef ImplicitEquationSolver < handle

    properties (Access = public)
        xsol
    end

    properties (Access = private)
        x0
        functionToSolve
    end

    methods (Access = public)

        function obj =  ImplicitEquationSolver(cParams)
            obj.init(cParams)
        end

        function r = solve(obj)
            F = obj.functionToSolve;
            [rub,rlb] = obj.findRbounds(F);
            r = fzero(F,[rlb,rub]);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.x0 = cParams.x0;
            obj.functionToSolve = cParams.functionToSolve;
        end

        function [rub,rlb] = findRbounds(obj,F)
            r0 = 1;
            F0 = F(r0);
            eps = 10^(-6);
            if F0 >= 0
                r1 = 1 + eps;
                F1 = F(r1);
                while F1 >= 0
                    rnew = obj.newPointBySecant(r0,r1,F0,F1);
                    r0 = r1;
                    F0 = F1;
                    r1 = rnew;
                    F1 = F(r1);
                end
                rub = r1;
                rlb = r0;
            else
                r0 = 1;
                r1 = 1 - eps;
                F1 = F(r1);
                while F1 <= 0
                    rnew = obj.newPointBySecant(r0,r1,F0,F1);
                    r0 = r1;
                    F0 = F1;
                    r1 = max(1e-12,rnew);
                    F1 = F(r1);
                end
                rub = r0;
                rlb = r1;
            end
        end

    end

    methods (Access = private, Static)

        function x2 = newPointBySecant(x0,x1,f0,f1)
            x2 = x1 - (x1-x0)/(f1 - f0)*f1;
        end

    end

end


text:
function old_default = set_matlab_bgl_default(varargin)
% SET_MATLAB_BGL_DEFAULT Sets a default option for the Matlab BGL interface
%
% old_default = set_matlab_bgl_default(options) or
% old_default = set_matlab_bgl_default(...) for key-value pair version
% options.istrans: the input matrices are already transposed [{0} | 1]
% options.nocheck: skip the input checking [{0} | 1]
% options.full2sparse: convert full matrices to sparse [{0} | 1]
%
% to get the current set of default options, call
% options = set_matlab_bgl_default()
%
% These options can make the Matlab BGL interface more efficient by
% eliminating the copying operations that occur between Matlab's structures
% and the BGL structures.  However, they are more difficult to use and are
% disabled by default.
%
% Generally, they are best used when you want to perform a large series of
% computations.
%
% Example:
%   % tranpose the matrix initially...
%   At = A'
%   old_options = set_matlab_bgl_default(struct('istrans',1));
%   % perform a bunch of graph work with At...
%   d1 = dfs(At,1); d2 = dfs(At,2); ...
%   % restore the old options
%   set_matlab_bgl_default(old_options);

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

%% History
%%

persistent default_options;
if ~isa(default_options,'struct')
    % initial default options
    default_options = struct('istrans', 0, 'nocheck', 0, 'full2sparse', 0);
end

if nargin == 0
    old_default = default_options;
else
    old_default = default_options;
    default_options = merge_options(default_options,varargin{:});
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef Entropy < Algorithm

    methods (Access = public)

        function obj = Entropy()
            obj.name = 'Entropy';
            obj.inputPort = DataType.kSignal;
            obj.outputPort = DataType.kFeature;
        end

        function result = compute(~,signal)

            alphabet = unique(signal);
            freq = zeros(size(alphabet));

            for symbol = 1:length(alphabet)
                freq(symbol) = sum(signal == alphabet(symbol));
            end

            P = freq / sum(freq);
            result = -sum(P .* log2(P));
        end

        function metrics = computeMetrics(~,input)
            n = size(input,1);
            flops = n * n;
            memory = n;
            outputSize = Constants.kFeatureBytes;
            metrics = Metric(flops,memory,outputSize);
        end
    end

end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function Problem = ITERmodify_LB_X0_UB(Problem)
X0 = Problem.x0(:);
LB = Problem.lb(:);
UB = Problem.ub(:);
XLabels = Problem.XLabels;

repeat = true;
while repeat
    disp(' ')
    BOUNDS = {'Index', 'Parameter'};
    for ii = 1:length(X0)
        BOUNDS{ii+1,1} = ii;
        BOUNDS{ii+1,2} = XLabels{ii};
    end
    disp(BOUNDS);
    LB_X0_UB = [LB X0 UB];
    openvar('LB_X0_UB');
    disp (' ');
    disp('[LB X0 UB] available for editing in array editor.')
    disp('Close array editor and type ''return'' in command window when finished editing')
    keyboard;
    LB = LB_X0_UB(:,1); UB = LB_X0_UB(:,3); X0 = LB_X0_UB(:,2);
    disp('LB, X0, and UB overwritten by edited values')
    repeat = false;
    if any(UB<LB)
        disp('*******************')
        disp('UB < LB at indices:')
        disp(find(UB<LB))
        beep
        repeat = true;
    elseif any(X0<LB) || any(X0>UB)
        disp('One or more seed values outside bounds. What would you like to do?')
        k = txtmenu([],'Edit in array editor',...
            'Trim seed to bounds','Expand bounds to seed');
        if k
            lowind = X0<LB; hiind = X0>UB;
            if k == 1
                X0(lowind) = LB(lowind);
                X0(hiind) = UB(hiind);
            elseif k ==2
                LB(lowind) = X0(lowind);
                UB(hiind) = X0(hiind);
            end
        else
            repeat = true;
        end
    end
end

Problem.x0 =    X0;
Problem.lb =    LB;
Problem.ub =    UB;


text:
%GENDATL Generation of Lithuanian classes
%
%  A = GENDATL(N,S)
%
% INPUT
%  N  Number of objects per class (optional; default: [50 50])
%  S  Standard deviation for the data generation (optional; default: 1)
%
% OUTPUT
%  A  Dataset
%
% DESCRIPTION
% Generation of Lithuanian classes, a 2-dimensional, 2-class dataset A
% of N objects according to the definition given by Raudys.
% The data is uniformly distributed along two sausages and is superimposed
% by a normal distribution with standard deviation S in all directions.
% Class priors are P(1) = P(2) = 0.5.
%
% SEE ALSO (PRTools Guide)
% DATASETS, PRDATASETS

% Copyright: M. Skurichina, R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: gendatl.m,v 1.3 2007/06/19 11:44:14 duin Exp $

function a = gendatl(N,s)
		if nargin < 1,
		prwarning(3,'Class cardinalities are not specified, assuming [50 50].');
		N = [50 50];
	end
	if nargin < 2,
		prwarning(4,'Standard deviation for the data generation is not specified, assuming 1.');
		s = 1;
	end
	if (length(N) == 1),
		N(2) = N(1);
	end;

	u1	= 2*pi/3*(rand(N(1),1)-0.5*ones(N(1),1));
	u2	= 2*pi/3*(rand(N(2),1)-0.5*ones(N(2),1));
	a 	= [[10*cos(u1) + s*randn(N(1),1) 10*sin(u1) + s*randn(N(1),1)]; ...
				[6.2*cos(u2) + s*randn(N(2),1) 6.2*sin(u2) + s*randn(N(2),1)]];
	lab = genlab(N);
	a 	= prdataset(a,lab,'name','Lithuanian Classes');
	a   = setlablist(a,[1 2]');
	a   = setprior(a,[0.5 0.5]);
return;

author: marianux

repo: ecg-kit

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


text:
function y = USFT_simple(x,shift,boxlen,center,w)
% USFT_simple -- 1d unequispaced Fourier transform
% Usage:
%   y = USFT_simple(x,shift,boxlen,center,w)
% Inputs:
%   x	   vector of length n
%   shift  vector of shifts
%   boxlen half-length of the window associated with shift
%   center  boolean variable
%   w      window of size 2*boxlen
% Outputs:
%   y      vector le length 2n
% Description:
%  Evaluates the FT at the points omega_(j,k) = shift(j) + k,
%                                               -boxlen <= k <boxlen
%
%  For a fixed value of j, and if center = 0,
%
%  y(k) = sum_{-n/2 <= t < n/2} exp(-i 2pi omega_k t/n) x(t), -n/2 <= t < n/2
%
%  If center = 1,
%
%  y(k) = sum_{0 <= t < n} exp(-i 2pi omega_k t/n) x(t), 0 <= t < n
%
%  See Also
%    Evaluate_FT, Adj_USFT_simple
%
% By Emmanuel candes, 2003-2004

  if nargin < 5,
    w = ones(1,2*boxlen);
  end

  if nargin < 4,
    center = 0;
  end

  n  = length(x); n2 = n/2;

  if center == 0,
    t  =  -n2:(n2-1);
  else
    t =  0:(n-1);
  end

  boxcnt = length(shift);
  col    = n2 + [(-boxlen+1):boxlen];

  % Recall that fft_mid0 operates along columns
  X  =  x*ones(1,boxcnt); % Make copies
  if center == 0,
    X  =  fft_mid0(X.* exp(-i*2*pi*t'*shift/n))./sqrt(n);
  else
    X  =  fftshift(fft(X.* exp(-i*2*pi*t'*shift/n)),1)./sqrt(n);
  end
  y    =  (w'*ones(1,boxcnt)).*X(col,:);
  y    =  y(:);

author: thomaskuestner


text:
function [sens] = ft_appendsens(cfg, varargin)

% FT_APPENDSENS concatenates multiple sensor definitions that have been processed
% separately.
%
% Use as
%   combined = ft_appendsens(cfg, sens1, sens2, ...)
%
% A call to FT_APPENDSENS results in the label, pos and ori fields to be
% concatenated, and the tra matrix to be merged. Any duplicate electrodes
% will be removed. The labelold and chanposold fields are kept under the
% condition that they are identical across the inputs.
%
% See also FT_ELECTRODEPLACEMENT, FT_ELECTRODEREALIGN, FT_DATAYPE_SENS,
% FT_APPENDDATA, FT_APPENDTIMELOCK, FT_APPENDFREQ, FT_APPENDSOURCE

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar    varargin
ft_preamble provenance varargin

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input data is valid for this function
% and ensure it is up to the latest standards
for i=1:length(varargin)
  assert(ft_datatype(varargin{i}, 'sens'), 'incorrect input, should be a grad/elec/opto structure');
  varargin{i} = ft_datatype_sens(varargin{i});
end

% do a basic check whether the senstype, units, and coordinate systems match
senstype = cell(1,length(varargin));
for i=1:length(varargin)
  senstype{i} = ft_senstype(varargin{i});
end
typematch = all(strcmp(senstype{1}, senstype));

unitmatch = 1;
if isfield(varargin{1}, 'unit')
  unit = cell(1,length(varargin));
  for i=1:length(varargin)
    unit{i} = varargin{i}.unit;
  end
  unitmatch = all(strcmp(unit{1}, unit));
end

coordsysmatch = 1;
if isfield(varargin{1}, 'coordsys')
  coordsys = cell(1,length(varargin));
  for i=1:length(varargin)
    coordsys{i} = varargin{i}.coordsys;
  end
  coordsysmatch = all(strcmp(coordsys{1}, coordsys));
end

if ~typematch || ~unitmatch || ~coordsysmatch
  ft_error('the senstype, units, or coordinate systems of the inputs do not match');
end

% keep these fields (when present) in the output
sens = keepfields(varargin{1}, {'type', 'unit', 'coordsys'});

% make inventory
haselecpos = 0;
hascoilpos = 0;
hascoilori = 0;
haschanori = 0;
hasoptopos = 0;
hastra     = 0;
haslabelold = 0;
haschanposold = 0;
for i=1:length(varargin)
  % the following fields should be present in any sens structure
  if isfield(varargin{i}, 'label')
    label{i} = varargin{i}.label(:); % ensure column orientation
  end
  if isfield(varargin{i}, 'chanpos')
    chanpos{i} = varargin{i}.chanpos;
  end

  % some the following fields are likely present in a sens structure
  if isfield(varargin{i}, 'elecpos') % EEG
    elecpos{i} = varargin{i}.elecpos;
    haselecpos = 1;
  end
  if isfield(varargin{i}, 'coilpos') % MEG
    coilpos{i} = varargin{i}.coilpos;
    hascoilpos = 1;
  end
  if isfield(varargin{i}, 'coilori') % MEG
    coilori{i} = varargin{i}.coilori;
    hascoilori = 1;
  end
  if isfield(varargin{i}, 'chanori') % MEG
    chanori{i} = varargin{i}.chanori;
    haschanori = 1;
  end
  if isfield(varargin{i}, 'optopos') % NIRS
    optopos{i} = varargin{i}.optopos;
    hasoptopos = 1;
  end
  if isfield(varargin{i}, 'tra') % tra
    tra{i} = varargin{i}.tra;
    hastra = 1;
  end

  % the following fields might be present in a sens structure
  if isfield(varargin{i}, 'labelold')
    labelold{i} = varargin{i}.labelold(:); % ensure column orientation
    haslabelold = 1;
  end
  if isfield(varargin{i}, 'chanposold')
    chanposold{i} = varargin{i}.chanposold;
    haschanposold = 1;
  end
end

% concatenate the main fields and remove duplicates
sens.label = cat(1,label{:});
[dum, labidx] = unique(sens.label, 'stable');
if ~isequal(numel(labidx), numel(sens.label))
  fprintf('removing duplicate labels\n')
  sens.label = sens.label(labidx);
end

sens.chanpos = cat(1,chanpos{:});
[dum, chanidx] = unique(sens.chanpos, 'rows', 'stable');
if ~isequal(numel(chanidx), size(sens.chanpos,1))
  fprintf('removing duplicate channels\n')
  sens.chanpos = sens.chanpos(chanidx,:);
  if ~isequal(labidx, chanidx) % check for matching order
    ft_error('inconsistent order or number of channel labels and positions')
  end
end
if ~isequal(numel(sens.label), size(sens.chanpos,1)) % check for matching number
  ft_error('inconsistent number of channel labels and positions')
end
if haschanori
  sens.chanori = cat(1,chanori{:});
  if ~isequal(numel(chanidx), size(sens.chanpos,1))
    sens.chanori = sens.chanori(chanidx,:); % chanori should match chanpos
  end
end

if hastra && ~any(cellfun(@isempty, tra))
  sens.tra = [];
  for t = 1:numel(tra)
    trarow = [1:size(tra{t},1)]+size(sens.tra,1);
    tracol = [1:size(tra{t},2)]+size(sens.tra,2);
    sens.tra(trarow, tracol) = tra{t};
  end
end

if haselecpos
  sens.elecpos = cat(1,elecpos{:});
  [dum, elecidx, elecidx2] = unique(sens.elecpos, 'rows', 'stable');
  if ~isequal(numel(elecidx), size(sens.elecpos,1))
    fprintf('removing duplicate electrodes\n')
    sens.elecpos = sens.elecpos(elecidx,:);
  end
  if isfield(sens, 'tra')
    % shape duplicates into a single column, if necessary
    for idx = 1:numel(elecidx)
      tmp(:, idx) = sum(sens.tra(chanidx, find(idx==elecidx2)),2); % find and take sum over duplicates
    end
    sens.tra = tmp;
    % check for expected size and non-empty rows
    if ~isequal(size(sens.tra,1), size(sens.chanpos,1)) || ~isequal(size(sens.tra,2), size(sens.elecpos,1)) ...
        || any(any(sens.tra,2)==0) % all channels need to be linked to their origins
      fprintf('removing inconsistent tra matrix\n')
      sens = rmfield(sens, 'tra');
    end
  end
end

if hasoptopos
  sens.optopos = cat(1,optopos{:});
  [dum, optoidx, optoidx2] = unique(sens.optopos, 'rows', 'stable');
  if ~isequal(numel(optoidx), size(sens.optopos,1))
    fprintf('removing duplicate optodes\n')
    sens.optopos = sens.optopos(optoidx,:);
  end
  if isfield(sens, 'tra')
    % shape duplicates into a single column, if necessary
    for idx = 1:numel(optoidx)
      tmp(:, idx) = sum(sens.tra(chanidx, find(idx==optoidx2)),2); % find and take sum over duplicates
    end
    sens.tra = tmp;
    % check for expected size and non-empty rows
    if ~isequal(size(sens.tra,1), size(sens.chanpos,1)) || ~isequal(size(sens.tra,2), size(sens.optopos,1)) ...
        || any(any(sens.tra,2)==0) % all channels need to be linked to their origins
      fprintf('removing inconsistent tra matrix\n')
      sens = rmfield(sens, 'tra');
    end
  end
end

if hascoilpos
  sens.coilpos = cat(1,coilpos{:});
  [dum, coilidx, coilidx2] = unique(sens.coilpos, 'rows', 'stable');
  if ~isequal(numel(coilidx), size(sens.coilpos,1))
    fprintf('removing duplicate coils\n')
    sens.coilpos = sens.coilpos(coilidx,:);
  end
  if isfield(sens, 'tra')
    % shape duplicates into a single column, if necessary
    for idx = 1:numel(coilidx)
      tmp(:, idx) = sum(sens.tra(chanidx, find(idx==coilidx2)),2); % find and take sum over duplicates
    end
    sens.tra = tmp;
    % check for expected size and non-empty rows
    if ~isequal(size(sens.tra,1), size(sens.chanpos,1)) || ~isequal(size(sens.tra,2), size(sens.coilpos,1)) ...
        || any(any(sens.tra,2)==0) % all channels need to be linked to their origins
      fprintf('removing inconsistent tra matrix\n')
      sens = rmfield(sens, 'tra');
    end
  end
end
if hascoilori
  sens.coilori = cat(1,coilori{:});
  if ~isequal(numel(coilidx), size(sens.coilori,1))
    sens.coilori = sens.coilori(coilidx,:); % coilori should match coilpos
  end
end

% keep the following fields only when identical across inputs
if haslabelold && all(isequal(labelold{1}, labelold{:})) % labeloldmatch
  sens.labelold = labelold{1};
end
if haschanposold && all(isequal(chanposold{1}, chanposold{:})) % chanposoldmatch
  sens.chanposold = chanposold{1};
end

% ensure up-to-date and consistent output sensor description
sens = ft_datatype_sens(sens);

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous varargin
ft_postamble provenance sens
ft_postamble history sens
ft_postamble savevar sens

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%
% Plot a motor trajectory on top of the image
%
% Input
%   I [105 x 105] image (double or binary) (large numbers are BLACK)
%   drawings: [nested cell] of strokes in motor space
function plot_image_only(I)
    assert(size(I,1)==105);
    assert(size(I,2)==105);
    newI = 1-I;
    newI(isinf(I)) = 0.5;
    image([1 105],[1 105],repmat(newI,[1 1 3]));
    set(gca,'YDir','reverse','XTick',[],'YTick',[]);
    xlim([1 105]);
    ylim([1 105]);
end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function ecop = ecopula(x)
%ECOPULA Empirical copula based on sample X.
%   ECOP = ECOPULA(X) returns bivariate empirical copula. Extension to
%   n dimensional empirical copula is straightforward.
%
%   Written by Robert Kopocinski, Wroclaw University of Technology,
%   for Master Thesis: "Simulating dependent random variables using copulas.
%   Applications to Finance and Insurance".
%   Date: 2007/05/12
%
%   Reference:
%      [1]  Durrleman, V. and Nikeghbali, A. and Roncalli, T. (2000) Copulas approximation and
%           new families, Groupe de Recherche Operationnelle Credit Lyonnais

[m n] = size(x);

y = sort(x);

for i=1:m
    for j=1:m
        ecop(i,j) = sum( (x(:,1)<=y(i,1)).*(x(:,2)<=y(j,2)) )/m;
    end
end


text:
function sphere_cubed_grid_points_display_test ( )

%*****************************************************************************80
%
%% SPHERE_CUBED_GRID_POINTS_DISPLAY_TEST tests SPHERE_CUBED_GRID_POINTS_DISPLAY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_CUBED_GRID_POINTS_DISPLAY_TEST\n' );
  fprintf ( 1, '  SPHERE_CUBED_GRID_POINTS_DISPLAY_TEST displays points\n' );
  fprintf ( 1, '  on a cubed sphere grid.\n' );

  n = 10;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of divisions on each face = %d\n', n );

  ns = sphere_cubed_grid_point_count ( n );
  fprintf ( 1, '  Total number of points = %d\n', ns );

  xyz = sphere_cubed_grid_points ( n, ns );

  filename = sprintf ( 'sphere_cubed_grid_points_f%d.png', n );

  sphere_cubed_grid_points_display ( ns, xyz, filename );

  return
end


text:
function [wtRes, delRes] = simpleOptKnock(model, targetRxn, deletions, geneDelFlag, minGrowth, doubleDelFlag)
% Simple `OptKnock` is used to check all one gene or reaction deletions for
% growth-coupled metabolite production
%
% USAGE:
%
%    [wtRes, delRes] = simpleOptKnock(model, targetRxn, deletions, geneDelFlag, minGrowth, doubleDelFlag)
%
% INPUTS:
%    model:            COBRA model structure
%    targetRxn:        Target metabolite production reaction
%
% OPTIONAL INPUTS:
%    deletions:        Set of gene or reaction deletions to consider for `KO`
%                      (Default = all reactions)
%    geneDelFlag:      Gene deletion flag (Default = false)
%    minGrowth:        Minimum `KO` growth rate (Default = 0.05)
%    doubleDelFlag:    Double deletions (Default = false)
%
% OUTPUTS:
%    wtRes:            Wild type results
%    delRes:           Deletion strain results
%
% .. Author: - Markus Herrgard 2/14/07
%
% The results structures have fields:
%
%    * `growth` - Growth rate of strain
%    * `minProd` - Minimum prod rate of target metabolite
%    * `maxProd` - Maximum prod rate of target metabolite

if (nargin < 3)
    deletions = model.rxns;
end
if (nargin < 4)
    geneDelFlag = false;
end
if (nargin < 5)
    minGrowth = 0.05;
end
if (nargin < 6)
    doubleDelFlag = false;
end

tol = 1e-7;

%rxnGeneMat is a required field for this function, so if it does not exist,
%build it.
if ~isfield(model,'rxnGeneMat')
    model = buildRxnGeneMat(model);
end

% Number of deletions
nDel = length(deletions);

% Wild type values
solWT = optimizeCbModel(model);
grRounded = floor(solWT.f/tol)*tol;
modelWT = changeRxnBounds(model,model.rxns(model.c==1),grRounded,'l');
modelWT = changeObjective(modelWT,targetRxn,1);
solMax = optimizeCbModel(modelWT);
solMin = optimizeCbModel(modelWT);

wtRes.growth = solWT.f;
wtRes.maxProd = solMax.f;
wtRes.minProd = solMin.f;

if (doubleDelFlag)
    growthRate = sparse(nDel,nDel);
    maxProd = sparse(nDel,nDel);
    minProd = sparse(nDel,nDel);
else
    growthRate = zeros(nDel,1);
    maxProd = zeros(nDel,1);
    minProd = zeros(nDel,1);
end

if (~doubleDelFlag)
    showprogress(0,'Simple OptKnock in progress ...');
end
t0 = cputime;
delCounter = 0;
for i = 1:nDel
    if (~doubleDelFlag)
        if mod(i,10) == 0
            showprogress(i/nDel);
        end
    end
    if (geneDelFlag)
        % Gene deletion
        modelKO = deleteModelGenes(model,deletions{i});
    else
        % Reaction deletion
        modelKO = changeRxnBounds(model,deletions{i},0,'b');
    end
    % Calculate optimal growth rate
    solKO = optimizeCbModel(modelKO);
    %fprintf('Single %s %f\n',deletions{i},solKO.f);
    growthRate(i,1) = solKO.f;
    if (solKO.f > minGrowth && solKO.stat == 1)
        % Max & min production of the metabolite at the optimal growth rate
        grRounded = floor(solKO.f/tol)*tol;
        modelKO = changeRxnBounds(modelKO,modelKO.rxns(modelKO.c==1),grRounded,'l');
        modelKO = changeObjective(modelKO,targetRxn,1);
        solMax = optimizeCbModel(modelKO,'max');
        solMin = optimizeCbModel(modelKO,'min');
        if (~doubleDelFlag)
            maxProd(i,1) = solMax.f;
            minProd(i,1) = solMin.f;
        %fprintf('%f %f\n',solMax.f,solMin.f);
        else
            maxProd(i,i) = solMax.f;
            minProd(i,i) = solMin.f;
            for j = i+1:nDel
                delCounter = delCounter+1;
                if mod(j,50) == 0
                    fComp = delCounter/(nDel*(nDel-1)/2);
                    fprintf('%d\t%f\t%f\n',delCounter,100*fComp,(cputime-t0)/60);
                end
                if (geneDelFlag)
                    % Gene deletion
                    modelKO2 = deleteModelGenes(model,deletions{i});
                    modelKO2 = deleteModelGenes(modelKO2,deletions{j});
                else
                    modelKO2 = changeRxnBounds(model,deletions{i},0,'b');
                    modelKO2 = changeRxnBounds(modelKO2,deletions{j},0,'b');
                end
                % Calculate optimal growth rate
                solKO2 = optimizeCbModel(modelKO2);
                growthRate(i,j) = solKO2.f;
                %fprintf('Double %s %s %f\n',deletions{i},deletions{j},solKO2.f);
                if (solKO2.f > minGrowth && solKO2.stat == 1)
                    grRounded2 = floor(solKO2.f/tol)*tol;
                    modelKO2 = changeRxnBounds(modelKO2,modelKO2.rxns(modelKO2.c==1),grRounded2,'l');
                    modelKO2 = changeObjective(modelKO2,targetRxn,1);
                    solMax2 = optimizeCbModel(modelKO2,'max');
                    solMin2 = optimizeCbModel(modelKO2,'min');
                    if (solMin2.f > 0)
                        fprintf('%s %s %f %f %f\n',deletions{i},deletions{j},solKO2.f,solMax2.f,solMin2.f);
                    end
                    maxProd(i,j) = solMax2.f;
                    minProd(i,j) = solMin2.f;
                end
            end
        end
    end
end

% Store results
delRes.maxProd = maxProd;
delRes.minProd = minProd;
delRes.growth = growthRate;

author: opencobra

repo: cobratoolbox


text:
%EDGSET Set Dirichlet edge constraints.
%
%   [ F, INDROW, AMAT, T ] = EDGSET( PROB, F, AMAT, ISYMM, SET_NULL )
%   Sets Dirichlet edge constraints conditions in the right hand side
%   load vector F and global matrix AMAT with the information in the finite
%   element problem struct PROB.
%
%       Input       Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       prob        struct                 Finite element problem struct
%       f           (neq,1)                Right hand side/load vector
%       amat        (n_a,n_a)              System matrix (sparse or triplet format)
%       isymm       scalar/{0}             Symmetrize BCs if applicable.
%       set_null    scalar/{0}             Set zeros in f vector.
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       f           (neq,1)                Modified right hand side/load vector
%       indrow      (neq,1)                Index to rows (dofs) in rhs which were set
%       amat        (n_a,n_a)              Modified system matrix
%       t           scalar                 Time spent in function

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function out = DN_Unique(x)
% DN_Unique     The proportion of the time series that are unique values
%
%---INPUTS:
%
% x, the input data vector

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

out = length(unique(x))/length(x);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function posteriori = slposterioritrue(condprops, nums, priori, op)
%SLPOSTERIORITRUE Computes the posteriori that samples belong to true class
%
% $ Syntax $
%   - posteriori = slposterioritrue(condprops, nums, priori)
%   - posteriori = slposterioritrue(condprops, nums, priori, 'log')
%
% $ Arguments $
%   - condprods:      the conditional probabilities of classes: C x n
%   - nums:           the number of samples belong to the classes: 1 x C
%   - priori:         the prior probabilities of classes: 1 x C
%   - posteriori:     the resulting posterior probabilities: 1 x n
%
% $ Description $
%   - posteriori = slposteriori(condprops, nums, priori) Computes the
%     posterior probability that the samples belong to the true class
%     according to the given conditional probabilities of all samples
%     to all classes and the priori of the classes.
%     In the input argument, each column in condprops represent the
%     conditional probabilities of that sample belong to all the
%     C classes. The samples from the same underlying classes should be
%     put together and the ownership is given by nums.
%     The priori can be given by an 1 x C row vector or [], which
%     means that all classes have equal prior.
%
%   - posteriori = slposterioritrue(condprops, nums, priori, 'log') means
%     that the input condprops are given by its logarithm.
%
% $ History $
%   - Created by Dahua Lin, on Sep 16th, 2006
%

%% parse and verify input arguments

if nargin < 2
    raise_lackinput('slposterioritrue', 2);
end

if ~isnumeric(condprops) || ndims(condprops) ~= 2
    error('sltoolbox:invalidarg', ...
        'The condprops should be a 2D numeric matrix');
end
[C, n] = size(condprops);

if ~isequal(size(nums), [1, C])
    error('sltoolbox:sizmismatch', ...
        'The nums should be an 1 x C row vector');
end
if nargin < 3 || isempty(priori)
    priori = [];
else
    if ~isequal(size(priori), [1, C])
        error('sltoolbox:sizmismatch', ...
            'The priori should be a an 1 x C row vector');
    end
end

if nargin >= 4 && strcmpi(op, 'log')
    is_log = true;
else
    is_log = false;
end

%% compute

if is_log
    if isempty(priori)
        P = sladdvec(condprops, -max(condprops, [], 1), 2);
    else
        P = sladdvec(condprops, log(priori)', 1);
        P = sladdvec(P, -max(condprops, [], 1), 2);
    end
    P = exp(P);
else
    if isempty(priori)
        P = condprops;
    else
        P = slmulvec(condprops, priori', 1);
    end
end

tP = sum(P, 1);
[sp, ep] = slnums2bounds(nums);
posteriori = zeros(1, n);

for k = 1 : C
    sk = sp(k);
    ek = ep(k);
    posteriori(sk:ek) = P(k, sk:ek);
end
posteriori = posteriori ./ tP;

author: lmthang


text:
function b = hdf5isvalid(fname)

val = H5F.is_hdf5(fname);
if val > 0
    b = true;
else
    b = false;
    fprintf('File %s is not an HDF5 file.\n', fname);
end

author: BUNPC


text:
function op = tfocs_prox( f, prox_f, VECTOR_SCALAR )
% OP = TFOCS_PROX( F, PROX_F )
%   combines F and PROX_F into the appropriate TFOCS-compatible object.
%
%   F is any function (with known proximity operator),
%   and PROX_F is its proximity operator, defined as:
%
%   PROX_F( Y, t ) = argmin_X  F(X) + 1/(2*t)*|| X - Y ||^2
%
%   To use this, please see the file PROX_L1 as an example
%
%   The basic layout of a file like PROX_L1 is as follows:
%   ( for the function F(X) = q*||X||_1 )
%
%       function op = prox_l1(q)
%       op = tfocs_prox( @f, @prox_f )
%
%         function v = f(x)
%           ... this function calculates the function f ...
%         end
%         function v = prox_f(y,t)
%           ... this function calculates the prox-function to f ...
%         end
%       end
%
%   Note: in the above template, the "end" statements are very important.
%
% OP = TFOCS_PROX( F, PROX_F, 'vector')
%   will signal the routine that it is OK to allow "vector" stepsizes
%   (this corresponds to solving
%       PROX_F( Y, T ) = argmin_X F(X) + 1/2( X-Y )'*diag(1./T)*( X-Y )
%    where T is a vector of the same size as X ).
% ... = TFOCS_PROX( F, PROX_F, 'scalar' )
%   will throw an error if a vector stepsize is attempted. This is the default.
%
%   Also, users may wish to test their smooth function
%   with the script TEST_NONSMOOTH
%
%   See also prox_l1, test_nonsmooth

if nargin < 3 || isempty(VECTOR_SCALAR)
    VECTOR_SCALAR = 'scalar';
end

if strcmpi(VECTOR_SCALAR,'scalar')
    op = @fcn_impl;
elseif strcmpi(VECTOR_SCALAR,'vector')
    op = @fcn_impl_vector;
else
    error('bad option for VECTOR_SCALAR parameter');
end

function [ v, x ] = fcn_impl(x, t )
    if nargin < 1,
        error( 'Not enough arguments.' );
    end
    if nargin == 2,
        if numel(t) ~= 1, error('The stepsize must be a scalar'); end
        x  = prox_f(x,t);
    elseif nargout == 2,
        error( 'This function is not differentiable.' );
    end
    v = f(x);
end

function [ v, x ] = fcn_impl_vector(x, t )
    if nargin < 1,
        error( 'Not enough arguments.' );
    end
    if nargin == 2,
        x  = prox_f(x,t);
    elseif nargout == 2,
        error( 'This function is not differentiable.' );
    end
    v = f(x);
end

end

% TFOCS v1.2 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2012 California Institute of Technology and CVX Research.
% See the file TFOCS/license.{txt,pdf} for full license information.

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function varargout = hat(varargin)
% VL_HAT  Hat operator
%   H = VL_HAT(OM) returns the skew symmetric matrix by taking the "hat"
%   of the 3D vector OM.
%
%   See also: VL_IHAT(), VL_HELP().
[varargout{1:nargout}] = vl_hat(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
% Attenuation Compensation Using Time Reversal Example
%
% This example demonstrates how the acoustic attenuation present in
% photoacoustic forward problem can be compensated for using time reversal
% image reconstruction. It builds on the 2D Time Reversal Reconstruction
% For A Circular Sensor Example.
%
% For a more detailed discussion of this example and the underlying
% techniques, see B. E. Treeby, E. Z. Zhang, and B. T. Cox, "Photoacoustic
% tomography in absorbing acoustic media using time reversal," Inverse
% Problems, vol. 26, no. 11, p. 115003, 2010.
%
% author: Bradley Treeby
% date: 6th September 2010
% last update: 25th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% =========================================================================
% FORWARD SIMULATION
% =========================================================================

% define the size of the simulation grid and the PML
simulation_size = 512;              % [grid points]
PML_size = 20;                      % [grid points]
x = 52e-3;                          % [m]
y = x;                              % [m]

% reduce the number of grid points in Nx and Ny by the size of the PML so
% that using 'PMLInside' set to false will still give the correct
% simulation size
Nx = simulation_size - 2*PML_size;  % [grid points]
Ny = Nx;                            % [grid points]
dx = x/Nx;                          % [m]
dy = dx;                            % [m]

% create the computational grid
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium_non_absorbing.sound_speed = 1510;	% [m/s]

% create a duplicate of the propagation medium structure and append the
% absorption properties
medium = medium_non_absorbing;
medium.alpha_power = 1.5;
medium.alpha_coeff = 3;              % [dB/(MHz^y cm)]

% store maximum supported frequency
f_max = kgrid.k_max*medium.sound_speed/(2*pi);

% load the shepp logan phantom (note if the simulation or PML sizes are
% changed, the loaded data will need to be resized)
load EXAMPLE_shepp_logan

% smooth the phantom and assign it to the initial pressure
shepp_logan = smooth(kgrid, shepp_logan, true);
source.p0 = shepp_logan;

% define a circular Cartesian sensor mask
sensor_radius = 25e-3;              % [m]
sensor_points = 200;
cart_sensor_mask = makeCartCircle(sensor_radius, sensor_points);
sensor.mask = cart_sensor_mask;

% create the time array used for the simulation, with t_max defined using
% Huygens' principle to avoid artifact trapping in the reconstruction
t_max = 2*sensor_radius/medium.sound_speed;
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, [], t_max);

% set the input options, switching off the smoothing (the input has already
% been smoothed), setting the PML to be outside the defined grid, casting
% to 'single' to speed up the example, and switching off visualisation
input_args = {'Smooth', false, 'PMLInside', false, 'PMLSize', PML_size, 'DataCast', 'single', 'PlotSim', false};

% run the forward simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% add noise to the recorded sensor data
signal_to_noise_ratio = 40;         % [dB]
sensor_data = addNoise(sensor_data, signal_to_noise_ratio, 'peak');

% =========================================================================
% IMAGE RECONSTRUCTION WITHOUT ABSORPTION COMPENSATION
% =========================================================================

% create a second computation grid for the reconstruction to avoid the
% inverse crime
PML_size = 25;                      % [grid points]
Nx = simulation_size - 2*PML_size;  % [grid points]
Ny = Nx;                            % [grid points]
dx = x/Nx;                          % [m]
dy = dx;                            % [m]
kgrid_recon = makeGrid(Nx, dx, Ny, dy);

% attach the original time array
kgrid_recon.t_array = kgrid.t_array;

% remove the initial pressure field from the source structure
source = rmfield(source, 'p0');

% create a continuous binary sensor mask with the same radius as the
% Cartesian sensor mask used in the forward simulation
pixel_radius = round(sensor_radius/kgrid_recon.dx);
binary_sensor_mask = makeCircle(kgrid_recon.Nx, kgrid_recon.Ny, floor(kgrid_recon.Nx/2) + 1, floor(kgrid_recon.Ny/2) + 1, pixel_radius);

% assign the sensor mask to the sensor structure
sensor.mask = binary_sensor_mask;

% interpolate the simulated sensor data onto the continuous binary sensor
% mask to remove any gaps and assign to the time reversal field
sensor.time_reversal_boundary_data = interpCartData(kgrid_recon, sensor_data, cart_sensor_mask, binary_sensor_mask, 'linear');

% re-assign the input options (the PML_size has changed)
input_args = {'Smooth', false, 'PMLInside', false, 'PMLSize', PML_size, 'DataCast', 'single', 'PlotSim', false};

% run the time-reversal reconstruction using the non absorbing medium
p0_recon = kspaceFirstOrder2D(kgrid_recon, medium_non_absorbing, source, sensor, input_args{:});

% =========================================================================
% CHOOSING THE CUTOFF FREQUENCY
% =========================================================================

% get the average frequency spectrum of the simulated sensor data
num_signals = length(sensor_data(:,1));
[as_f, as] = spect(sensor_data(1, :), 1/dt);
for index = 2:num_signals
    [as_f, sp] = spect(sensor_data(index, :), 1/dt);
    as = as + sp;
end
as = as/num_signals;

% compute the relative power spectrum
ps = log10(as.^2);
offset = max(ps(:));
ps = ps - offset;

% get the frequency spectrum of a single measurement
as_sing = spect(sensor_data(1,:), 1/dt);
ps_sing = log10(as_sing.^2) - offset;

% scale the frequency variable
[f_sc, scale, prefix] = scaleSI(as_f(end));

% define the cutoff frequency for the filter
f_cutoff = 3e6;

% =========================================================================
% IMAGE RECONSTRUCTION WITH ABSORPTION COMPENSATION
% =========================================================================

% create the filter to regularise the absorption parameters
medium.alpha_filter = getAlphaFilter(kgrid_recon, medium, f_cutoff);

% reverse the sign of the absorption proportionality coefficient
medium.alpha_sign = [-1, 1];        % [absorption, dispersion];

% run the time-reversal reconstruction
p0_recon_compensated = kspaceFirstOrder2D(kgrid_recon, medium, source, sensor, input_args{:});

% =========================================================================
% VISUALISATION
% =========================================================================

% plot the amplitude spectrum
figure;
plot(as_f*scale, ps_sing, 'r-');
ylabel('Power Spectrum [dB]');
xlabel(['Frequency [' prefix 'Hz]']);
hold on;
plot(as_f*scale, ps, 'k-');
ylim = [-8 0];
plot([f_cutoff*scale f_cutoff*scale], [ylim(1) ylim(2)], 'k--');
plot([f_max*scale f_max*scale], [ylim(1) ylim(2)], 'k--');
set(gca, 'XLim', [0 15], 'Ylim', ylim);

% plot the filter, reducing the number of plotting points by a factor of 10
figure;
mesh(medium.alpha_filter(1:10:end, 1:10:end), 'EdgeColor', 'black');
axis tight;
view([-29, 46]);

% set the pressure outside the sensor mask to be zero
binary_sensor_map = makeDisc(kgrid_recon.Nx, kgrid_recon.Ny, kgrid_recon.Nx/2, kgrid_recon.Ny/2, pixel_radius - 50);
p0_recon(binary_sensor_map ~= 1) = 0;
p0_recon_compensated(binary_sensor_map ~= 1) = 0;

% shrink the data sets
phantom_padding = 50;
shepp_logan = shepp_logan(1+ phantom_padding:end-phantom_padding, 1+ phantom_padding:end-phantom_padding);
p0_recon = p0_recon(1+ phantom_padding:end-phantom_padding, 1+ phantom_padding:end-phantom_padding);
p0_recon_compensated = p0_recon_compensated(1+ phantom_padding:end-phantom_padding, 1+ phantom_padding:end-phantom_padding);

% plot the initial pressure distribution
figure;
subplot(2, 2, 1), imagesc(kgrid_recon.y_vec(1+ phantom_padding:end-phantom_padding)*1e3, kgrid_recon.x_vec(1+ phantom_padding:end-phantom_padding)*1e3, shepp_logan, [0 1]);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
colormap(flipud(gray));

% plot the reconstruction without attenuation compensation
subplot(2, 2, 2), imagesc(kgrid_recon.y_vec(1+ phantom_padding:end-phantom_padding)*1e3, kgrid_recon.x_vec(1+ phantom_padding:end-phantom_padding)*1e3, p0_recon, [0 1]);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
colormap(flipud(gray));

% plot the reconstruction with attenuation compensation
subplot(2, 2, 3), imagesc(kgrid_recon.y_vec(1+ phantom_padding:end-phantom_padding)*1e3, kgrid_recon.x_vec(1+ phantom_padding:end-phantom_padding)*1e3, p0_recon_compensated, [0 1]);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
colormap(flipud(gray));

% plot a reconstruction profile
figure;
plot(kgrid.y_vec(1+ phantom_padding:end-phantom_padding)*1e3, shepp_logan(end/2,:), 'k');
hold on;
plot(kgrid_recon.y_vec(1+ phantom_padding:end-phantom_padding)*1e3, p0_recon(end/2,:), 'r-');
plot(kgrid_recon.y_vec(1+ phantom_padding:end-phantom_padding)*1e3, p0_recon_compensated(end/2, :), 'b-');
set(gca, 'XLim', [-20, 20], 'YLim', [-0.1 1.1]);
ylabel('Pressure Magnitude [au]');
xlabel('y-position [mm]');
legend('Original', 'No Compensation', 'With Attenuation Compensation', 'Location', 'Best');

repo: Aero-matlab


text:
function r = i4int_to_r8int ( imin, imax, i, rmin, rmax )

%*****************************************************************************80
%
%% I4INT_TO_R8INT maps an integer interval to an R8 interval.
%
%  Formula:
%
%    R := RMIN + ( RMAX - RMIN ) * ( I - IMIN ) / ( IMAX - IMIN )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer IMIN, IMAX, the range.
%
%    Input, integer I, the integer to be converted.
%
%    Input, real RMIN, RMAX, the range.
%
%    Output, real R, the corresponding value in [RMIN,RMAX].
%
  if ( imax == imin )

    r = 0.5 * ( rmin + rmax );

  else

    r = (  ( imax - i        ) * rmin   ...
        +  (        i - imin ) * rmax ) ...
        /  ( imax     - imin );

  end

  return
end

author: johannesgerer


text:
function [f,h] = kafbox_profiler_plotconvergence(algorithms,...
    mse_curves,resinds)

figure; hold all
set(gcf,'Position',[200, 200, 500 300])

titles = cell(length(algorithms),1);

for i=1:size(resinds,1)
    algo = algorithms{resinds(i,1)};

    ls = '-'; % line style
    if isfield(algo.figstyle,'line')
        ls = algo.figstyle.line;
    end

    lw = 1; % line width
    if isfield(algo.figstyle,'linewidth')
        lw = algo.figstyle.linewidth;
    end

    curve = mse_curves{resinds(i,1)}{resinds(i,2)};
    xs= ~isnan(curve);
    inds = 1:length(mse_curves{resinds(i,1)}{resinds(i,2)});
    plot(inds(xs),10*log10(curve(xs)),'color',algo.figstyle.color,...
        'LineWidth',lw,'LineStyle',ls)

    titles{i} = algo.name;
end

h = legend(titles);
grid on; box on
xlabel('iteration');
ylabel('(N)MSE');

% Expand figure axes
f = gcf;
style = hgexport('factorystyle'); style.Bounds = 'tight';
hgexport(f,'-clipboard',style,'applystyle', true);
drawnow;

repo: kafbox


text:
function m = size(t,idx)
%SIZE Size of a ttensor.
%
%   D = SIZE(T) returns the size of the tensor.
%
%   I = size(T,DIM) returns the size of the dimension specified by
%   the scalar DIM.
%
%   See also TTENSOR, TTENSOR/NDIMS.
%
%MATLAB Tensor Toolbox.
%Copyright 2010, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by Brett Bader and Tamara Kolda.
% http://csmr.ca.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2010) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in tensor_toolbox/LICENSE.txt
% $Id: size.m,v 1.8 2010/03/19 23:46:31 tgkolda Exp $

if ndims(t) == 0
    m = [];
end

if exist('idx','var')
    m = size(t.u{idx}, 1);
else
    for i = 1 : ndims(t)
	m(i) = size(t.u{i}, 1);
    end
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [ a, b ] = p08_lim ( dim_num )

%*****************************************************************************80
%
%% P08_LIM returns the integration limits for problem 08.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Output, real A(DIM_NUM), B(DIM_NUM), the lower and upper
%    limits of integration.
%
  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  return
end

repo: jburkardt-m


text:
function Tracks = Pops2Tracks(WolfPops, UAV)
%POPS2TRACKS 将狼群转换为航迹

SearchAgents = size(WolfPops.Pos, 1);          % 种群数量
UAVnum = UAV.num;                                      % 无人机个数
dim = UAV.PointDim;                                      % 仿真维度
v = WolfPops.Pos(:, end-UAVnum+1:end);   % 协同无人机速度
P = WolfPops.Pos(:, 1:end-UAVnum);            % 协同无人机航迹 xy

Tracks = cell(SearchAgents, 1);
for agent = 1 : SearchAgents
    a.V = v(agent, :)';
    P_a = P(agent, :);
    a.P = cell(UAVnum, 1);
    for i =1:UAVnum
        PointNum = UAV.PointNum(i);
        % 三维仿真这里得改 PointNum*dim
        P_ai = P_a(1 : PointNum*dim);
        P_ai = reshape(P_ai, dim, PointNum);
        P_a = P_a(PointNum*dim+1 : end);
        a.P(i) =  {P_ai};
    end
    Tracks(agent) = {a};
end

end

repo: Grey-Wolf-Optimizer-for-Path-Planning


text:
function varargout = spy(X)
%SPY (overloaded)

 if isa(X,'blkvar')
    X = sdpvar(X);
 end

Z = reshape(sum(abs(X.basis),2),X.dim(1),X.dim(2));
Z = Z~=0;
if nargout==0
    spy(Z)
else
    varargout{1}=Z;
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
classdef ApoapsisAltitudeTermCondition < AbstractEventTerminationCondition
    %PeriapsisAltitudeTermCondition Summary of this class goes here
    %   Detailed explanation goes here

    properties
        apoalt(1,1) double = 0; %km
        bodyInfo KSPTOT_BodyInfo
    end

    methods
        function obj = ApoapsisAltitudeTermCondition(apoalt)
            obj.apoalt = apoalt;
        end

        function evtTermCondFcnHndl = getEventTermCondFuncHandle(obj)
            evtTermCondFcnHndl = @(t,y) obj.eventTermCond(t,y);
        end

        function initTermCondition(obj, initialStateLogEntry)
            obj.bodyInfo = initialStateLogEntry.centralBody;
        end

        function name = getName(obj)
            name = sprintf('Apoapsis Altitude (%.3f km)', obj.apoalt);
        end

        function tf = shouldBeReinitOnRestart(obj)
            tf = true;
        end

        function params = getTermCondUiStruct(obj)
            params = struct();

            params.paramName = 'Apo. Altitude';
            params.paramUnit = 'km';
            params.useParam = 'on';
            params.useStages = 'off';
            params.useTanks = 'off';
            params.useEngines = 'off';
            params.useStopwatches = 'off';

            params.value = obj.apoalt;
            params.refStage = LaunchVehicleStage.empty(1,0);
            params.refTank = LaunchVehicleEngine.empty(1,0);
            params.refEngine = LaunchVehicleEngine.empty(1,0);
            params.refStopwatch = LaunchVehicleStopwatch.empty(1,0);
        end

        function optVar = getNewOptVar(obj)
            optVar = ApoapsisAltitudeOptimizationVariable(obj);
        end

        function optVar = getExistingOptVar(obj)
            optVar = obj.optVar;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end
    end

    methods(Static)
        function termCond = getTermCondForParams(paramValue, stage, tank, engine, stopwatch)
            termCond = ApoapsisAltitudeTermCondition(paramValue);
        end
    end

    methods(Access=private)
        function [value,isterminal,direction] = eventTermCond(obj, t,y)
            rVect = y(1:3);
            vVect = y(4:6);
            cartElem = CartesianElementSet(t, rVect(:), vVect(:), obj.bodyInfo.getBodyCenteredInertialFrame());
            kepElem = cartElem.convertToFrame(obj.frame).convertToKeplerianElementSet();

            value =  kepElem.getAltitudeApoapsis() - obj.apoalt;
            isterminal = 1;
            direction = 0;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function X = randnd(beta,varargin)
%% X = randnd(beta,varargin)
% This function is a "decorator" for the Matlab in-build 'randn' function.
% It takes an additional first argument (beta) which shapes the spectral
% characteristic of the data (in all dimensions) as f^beta. The output data
% is scaled to zero mean and unit standard deviation.
%
% For help on 'varargin' go to Matlab help for the 'randn' function.
%
% EXAMPLE 1 - Input robustness
% A = randnd(-1,[2^5+1 1]); % Handle odd dimensionality
% A = randnd((rand-0.5)*4,[134 12 1 1 1 26 1 9]); % Handle large number of dimentions
% A = randnd(-1,[1 0 3]); % Handle 0's
% A = randnd(2,[1 1 1 256]); % Handle leading singleton dimensions
%
% EXAMPLE 2 - Pink map
% A = randnd(-1,2^8); % Generate 256 x 256 grid with random 1/f (pink) noise
% contourf(A,'edgecolor','none'); axis square; % Display it as a contour plot
%
% EXAMPLE 3 - Brownian cube
% A = randnd(-2,[2^4 2^4 2^4]); % Generate 16 x 16 x 16 cube of f^-2 (Brownian) noise
% [X,Y,Z] = meshgrid(1:2^4); % Generate 3d meshgrid
% scatter3(X(:),Y(:),Z(:),100,A(:),'filled');
% axis vis3d; axis off;
%
% EXAMPLE 4 - Single blue noise
% p = single(magic(10)); % Single precision matrix of interest
% A = randnd(1,size(p),'like',p); % Generate f (blue) noise 'like' p
% isa(A,'single') % X is also a single precision random number
%
% EXAMPLE 5 - Proof that beta scales power spectral density correctly
% beta = -2;
% A = randnd(beta,[1 2^11]);
% Afft = abs(fft(A.^2)); % Get the magnitude of power (i.e. signal^2) spectrum
% Afft = Afft(2:length(Afft)/2); % Get the spectrum between DC and Nyquist
% k = 1:length(Afft); % Create wavenumber (frequency) vector
% P = polyfit(log(k),log(Afft),1); % Fit a straight line to the loglog plot
% loglog(k,Afft,'.',k,exp(P(2))*k.^P(1),'-'); grid on; % Plot on loglog axis
% xlim([min(k) max(k)]); title(['Beta = ',num2str(beta),'; Slope = ',num2str(P(1))]);
%
%
% Based on similar functions by Jon Yearsley and Hristo Zhivomirov
% Written by Marcin Konowalczyk
% Timmel Group @ Oxford University

%% Parse the input
narginchk(0,Inf); nargoutchk(0,1);

if nargin < 2 || isempty(beta); beta = 0; end % Default to white noise
assert(isnumeric(beta) && isequal(size(beta),[1 1]),'''beta'' must be a number');
assert(-2 <= beta && beta <= 2,'''beta'' out of range'); % Put on reasonable bounds

%% Generate N-dimensional white noise with 'randn'
X = randn(varargin{:});
if isempty(X); return; end; % Usually happens when size vector contains zeros

% Squeeze prevents an error if X has more than one leading singleton dimension
% This is a slight deviation from the pure functionality of 'randn'
X = squeeze(X);

% Return if white noise is requested
if beta == 0; return; end;

%% Generate corresponding N-dimensional matrix of multipliers
N = size(X);
% Create matrix of multipliers (M) of X in the frequency domain
M = [];
for j = 1:length(N)
    n = N(j);

    if (rem(n,2)~=0) % if n is odd
        % Nyquist frequency bin does not show up in odd-numbered fft
        k = ifftshift(-(n-1)/2:(n-1)/2);
    else
        k = ifftshift(-n/2:n/2-1);
    end

    % Spectral multipliers
    m = (k.^2)';

    if isempty(M);
        M = m;
    else
        % Create the permutation vector
        M_perm = circshift(1:length(size(M))+1,1,2);
        % Permute a singleton dimension to the beginning of M
        M = permute(M,M_perm);
        % Add m along the first dimension of M
        M = bsxfun(@plus,M,m);
    end
end
clear j M_perm k n m
% Reverse M to match X (since new dimensions were being added form the left)
M = permute(M,length(size(M)):-1:1);
assert(isequal(size(M),size(X)),'Bad programming error'); % This should never occur

% Shape the amplitude multipliers by (beta/2) which corresponds to shaping
% the power by beta
M = M.^(beta/2);

% Set the DC component to zero
M(1,1) = 0;

%% Multiply X by M in frequency domain
Xstd = std(X(:));
Xmean = mean(X(:));
X = real(ifftn(fftn(X).*M));

% Force zero mean unity standard deviation
X = X - mean(X(:));
X = X./std(X(:));

% Restore the standard deviation and mean from before the spectral shaping.
% This ensures the random sample from randn is truly random. After all, if
% the mean was always exactly zero it would not be all that random.
X = X + Xmean;
X = X.*Xstd;
end

repo: DBCNN


text:
% this script tests GLCM features between CERR and pyradiomics on a wavelet filtered image.
%
% RKP, 03/22/2018

%% Load image
glcmParamFileName = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing','tests_for_cerr','test_glcm_radiomics_extraction_settings.json');
cerrFileName = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing','data_for_cerr_tests','CERR_plans','head_neck_ex1_20may03.mat.bz2');

planC = loadPlanC(cerrFileName,tempdir);
indexS = planC{end};

paramS = getRadiomicsParamTemplate(glcmParamFileName);
strNum = getMatchingIndex(paramS.structuresC{1},{planC{indexS.structures}.structureName});
scanNum = getStructureAssociatedScan(strNum,planC);

scanType = 'wavelet';
dirString = 'HHH';

%% Calculate features using CERR
harFeat3DdirS = calcGlobalRadiomicsFeatures...
            (scanNum, strNum, paramS, planC);
harlCombS = harFeat3DdirS.Wavelets_Coif1__HHH.glcmFeatS.AvgS;
cerrGlcmV = [harlCombS.autoCorr, harlCombS.jointAvg, harlCombS.clustPromin, harlCombS.clustShade, harlCombS.clustTendency, ...
harlCombS.contrast, harlCombS.corr, harlCombS.diffAvg, harlCombS.diffEntropy, harlCombS.diffVar, harlCombS.dissimilarity, ...
harlCombS.energy, harlCombS.jointEntropy, harlCombS.invDiff, harlCombS.invDiffMom, harlCombS.firstInfCorr, ...
harlCombS.secondInfCorr, harlCombS.invDiffMomNorm, harlCombS.invDiffNorm, harlCombS.invVar, ...
harlCombS.sumAvg, harlCombS.sumEntropy, harlCombS.sumVar];

% %% Calculate features using pyradiomics
% % image and mask for a structure
% testM = single(planC{indexS.scan}(scanNum).scanArray) - ...
%     single(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset);
% mask3M = zeros(size(testM),'logical');
% [rasterSegments, planC, isError] = getRasterSegments(strNum,planC);
% [maskBoundBox3M, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);
% mask3M(:,:,uniqueSlices) = maskBoundBox3M;
%
% dx = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
% dy = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
% dz = mode(diff([planC{indexS.scan}(scanNum).scanInfo(:).zValue]));
% pixelSize = [dx dy dz]*10;
%
% teststruct = PyradWrapper(testM, mask3M, pixelSize, scanType, dirString);
%
% %teststruct = PyradWrapper(testM, mask3M, scanType, dirString);
% pyradGlcmNamC = {'Autocorrelation', 'JointAverage', 'ClusterProminence', 'ClusterShade',  'ClusterTendency', ...
%     'Contrast', 'Correlation', 'DifferenceAverage', 'DifferenceEntropy', 'DifferenceVariance', 'Dissimilarity', ...
%     'JointEnergy', 'JointEntropy','Id','Idm', 'Imc1' , ...
%     'Imc2', 'Idmn','Idn','InverseVariance', 'sumAverage', 'SumEntropy', 'sumVariance'};
%
% pyradGlcmNamC = strcat(['wavelet','_', dirString,'_glcm_'],pyradGlcmNamC);
% pyRadGlcmV = [];
% for i = 1:length(pyradGlcmNamC)
%     if isfield(teststruct,pyradGlcmNamC{i})
%         pyRadGlcmV(i) = teststruct.(pyradGlcmNamC{i});
%     else
%         pyRadGlcmV(i) = NaN;
%     end
% end
%
% %% Compare
% glcmDiffV = (cerrGlcmV - pyRadGlcmV) ./ cerrGlcmV * 100

%% Compare using previously calculated values of pyradiomics glcm
saved_pyRadGlcmV = [929.516284841440,30.4950632056943,5316.83199378276,-0.303806680516285,20.9971321089971,22.7275310589439,-0.0380659828547803,2.78099962843313,3.04970358686672,14.7880425011987,NaN,0.0446949242970363,6.28336538646016,0.479852104826365,0.424210035556516,-0.0921069324860780,0.664903434331990,0.993049866792659,0.955511453133062,0.367836063144979,NaN,3.75942560575858,NaN];
glcmDiffV = (cerrGlcmV - saved_pyRadGlcmV) ./ cerrGlcmV * 100

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function centroid = findBallFcn(greenBall1, thresh, imageType, axH)

%% Find Green Object
% This script reads in an image file and then attempts to find a green
% object in the image. It is designed to find one green ball and highlight
% that ball on the original image

% Copyright 2013 The MathWorks, Inc.

error(nargchk(3, 4, nargin, 'struct'));

if isempty(greenBall1)
    return;
end
if ischar(greenBall1)
    greenBall1 = imread(greenBall1);
end
%% Step 1: Read image into MATLAB
% First we read the specified image from the file and bring it into MATLAB
% as a variable. We also display the image to ensure it is correct.
% greenBall1 = imread('greenBall3.jpg');
% imtool(greenBall1);

%% Step 2: Identify Unique Characteristics of Object of Interest

%%
% Extract each color
% Next we using indexing to extract three 2D matrices from the 3D image
% data corresponding to the red, green, and blue components of the image.
r = greenBall1(:, :, 1);
g = greenBall1(:, :, 2);
b = greenBall1(:, :, 3);

%%
% View different color planes
% figure
% subplot(2,2,1),imagesc(r)
% subplot(2,2,2),imagesc(g)
% subplot(2,2,3),imagesc(b)

%%
% Calculate Green
% Then we perform an arithmetic operation on the matrices as a whole to try
% to create one matrix that represents an intensity of green.
justGreen = g - r/2 - b/2;
% colorPlanesPlot(r,g,b,justGreen);

%%
% close all

%% Step 3: Isolate Object of Interest

%%
% Threshold the image
% Now we can set a threshold to separate the parts of the image that we
% consider to be green from the rest.
% bw = justGreen > 50;
% imagesc(bw);
% colormap(gray);
if nargin == 4
    bw = justGreen > thresh;
else
    bw = justGreen > 80;
end
%%
% Remove small unwanted objects
% We can use special functions provided by the Image Processing toolbox to
% quickly perform common image processing tasks. Here we are using
% BWAREAOPEN to remove groups of pixels less than 30.
ball1 = bwareaopen(bw, 30);
% imagesc(ball1);

%% Step 4: Find center of green object
% Now we are using REGIONPROPS to extract the centroid of the group of
% pixels representing the ball.
% s  = regionprops(ball1, {'centroid','area'});
% imshow(greenBall1);
% if isempty(s)
%   title('No ball found!');
% else
%   [~, id] = max([s.Area]);
%   hold on, plot(s(id).Centroid(1),s(id).Centroid(2),'wp','MarkerSize',20,'MarkerFaceColor','r'), hold off
%   title(['Center location is (',num2str(s(id).Centroid(1),4),', ',num2str(s(id).Centroid(2),4),')'])
% end
s  = regionprops(ball1, {'centroid','area'});
if isempty(s)
    centroid = [];
else
    [maxArea, id] = max([s.Area]); %#ok<ASGLU>
    centroid = s(id).Centroid;
end
switch imageType
   case 'video'
      imshow(greenBall1, 'Parent', axH);
      if ~isempty(centroid)
         line(centroid(1), centroid(2), 'Parent', axH, 'Color', 'w', 'Marker', 'p', 'MarkerSize', 20, 'MarkerFaceColor', 'r')
      end
   case 'bw'
      imshow(ball1, 'Parent', axH);
end
%% Step 5: Verify estimated location
% Finally we will plot the center on the original image to clearly evaluate
% how well we have found the center.
% imagesc(greenBall1);
% hold on, plot(s(id).Centroid(1),s(id).Centroid(2),'wp','MarkerSize',20,'MarkerFaceColor','r'), hold off

author: Sable

repo: mcbench-benchmarks


text:
function month_name = month_to_month_name_common ( m )

%*****************************************************************************80
%
%% MONTH_TO_MONTH_NAME_COMMON returns the name of a Common month.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the month index.
%
%    Output, string MONTH_NAME, the month name.
%
  name = [ ...
    'January  '; 'February '; 'March    '; 'April    '; ...
    'May      '; 'June     '; 'July     '; 'August   '; ...
    'September'; 'October  '; 'November '; 'December ' ];

  if ( m < 1 || 12 < m )

    month_name = '?????????';

  else

    month_name = name(m,:);

  end

  return
end