text:
function fem2d_bvp_serene_test01 ( )

%*****************************************************************************80
%
%% FEM2D_BVP_SERENE_TEST01 carries out test case #1.
%
%  Discussion:
%
%    Use A1, C1, F1, EXACT1, EXACT_UX1, EXACT_UY1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 June 2014
%
%  Author:
%
%    John Burkardt
%
  nx = 5;
  ny = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_BVP_SERENE_TEST01\n' );
  fprintf ( 1, '  Solve - del ( A del U ) + C U = F \n' );
  fprintf ( 1, '  on the unit square with zero boundary conditions.\n' );
  fprintf ( 1, '  A1(X,Y) = 1.0\n' );
  fprintf ( 1, '  C1(X,Y) = 0.0\n' );
  fprintf ( 1, '  F1(X,Y) = 2*X*(1-X)+2*Y*(1-Y).\n' );
  fprintf ( 1, '  U1(X,Y) = X * ( 1 - X ) * Y * ( 1 - Y )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The grid uses %d by %d nodes.\n', nx, ny );
  node_num = fem2d_bvp_serene_node_num ( nx, ny );
  fprintf ( 1, '  The number of nodes is %d\n', node_num );
%
%  Geometry definitions.
%
  x = linspace ( 0.0, 1.0, nx );
  y = linspace ( 0.0, 1.0, ny );

  show11 = 0;
  u = fem2d_bvp_serene ( nx, ny, @a1, @c1, @f1, x, y, show11 );

  if ( nx * ny <= 25 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '     I     J    X         Y         U         Uexact    Error\n' );
    fprintf ( 1, '\n' );

    k = 0;

    for j = 1 : ny

      if ( mod ( j, 2 ) == 1 )
        inc = 1;
      else
        inc = 2;
      end

      for i = 1 : inc : nx
        k = k + 1;
        uexact = exact1 ( x(i), y(j) );
        fprintf ( 1, '  %4d  %4d  %8f  %8f  %8f  %8f  %8e\n', ...
          i, j, x(i), y(j), u(k), uexact, abs ( u(k) - uexact ) );
      end
    end

  end

  e1 = fem2d_l1_error_serene ( nx, ny, x, y, u, @exact1 );
  e2 = fem2d_l2_error_serene ( nx, ny, x, y, u, @exact1 );
  h1s = fem2d_h1s_error_serene ( nx, ny, x, y, u, @exact_ux1, @exact_uy1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  l1 error   = %g\n', e1 );
  fprintf ( 1, '  L2 error   = %g\n', e2 );
  fprintf ( 1, '  H1S error  = %g\n', h1s );

  return
end
function value = a1 ( x, y )

%*****************************************************************************80
%
%% A1 evaluates A function #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the evaluation point.
%
%    Output, real VALUE, the value of A(X).
%
  value = 1.0;

  return
end
function value = c1 ( x, y )

%*****************************************************************************80
%
%% C1 evaluates C function #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the evaluation point.
%
%    Output, real VALUE, the value of C(X).
%
  value = 0.0;

  return
end
function value = exact1 ( x, y )

%*****************************************************************************80
%
%% EXACT1 evaluates exact solution #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the evaluation point.
%
%    Output, real VALUE, the value of the solution.
%
  value = x .* ( 1.0 - x ) .* y .* ( 1.0 - y );

  return
end
function value = exact_ux1 ( x, y )

%*****************************************************************************80
%
%% EXACT_UX1 evaluates the derivative dUdX of exact solution #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the evaluation point.
%
%    Output, real VALUE, the value of dUdX.
%
  value = ( 1.0 - 2.0 * x ) .* ( y - y .* y );

  return
end
function value = exact_uy1 ( x, y )

%*****************************************************************************80
%
%% EXACT_UY1 evaluates the derivative dUdY of exact solution #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the evaluation point.
%
%    Output, real VALUE, the value of dUdX.
%
  value = ( x - x .* x ) .* ( 1.0 - 2.0 * y );

  return
end
function value = f1 ( x, y )

%*****************************************************************************80
%
%% F1 evaluates right hand side function #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the evaluation point.
%
%    Output, real VALUE, the value of the right hand side.
%
  value = 2.0 * x .* ( 1.0 - x ) ...
        + 2.0 * y .* ( 1.0 - y );

  return
end

author: johannesgerer


text:
function img = discreteHalfPlane(varargin)
%DISCRETEHALFPLANE Discretize a half plane
%
%   A Halfplane is the set of point delimited by a straight line.
%   Only the points located 'on the left' of the line belong to the
%   halfplane.
%
%   IMG = discreteHalfPlane(DIM, LINE)
%   DIM is the size of image, with the format [x0 dx x1;y0 dy y1]
%   LINE is a 1x4 array of the form [x0 y0 dx dy], x0 and y0 being
%   coordinate of a point belonging to the boundary line, and dx and dy
%   being direction vectors of the boundary line of the halfplane.
%
%   IMG = discreteHalfPlane(DIM, POINT, DIRECTION)
%   POINT is a point belonging to the boundary line
%   DIRECTION is a 2x1 vector containing direction vector of the boundary
%   line.
%
%   IMG = discreteHalfPlane(LX, LY, ...);
%   Specifes the pixels coordinates with the two row vectors LX and LY.
%
%   Example
%   img = discreteHalfPlane([1 1 100;1 1 100], [50 50], 30, 10);
%
%   See Also
%     imShapes, discreteDisc, discreteSquare
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2006-10-12
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%   HISTORY
%   04/01/2007: concatenate transforms before applying them
%   04/03/2009: use meshgrid
%   29/04/2009: update transforms
%   29/05/2009: use more possibilities for specifying grid
%   22/01/2010: fix auto center with odd image size

% compute coordinate of image voxels
[lx, ly, varargin] = parseGridArgs(varargin{:});
[x, y]   = meshgrid(lx, ly);

% default parameters
center = [lx(ceil(end/2)) ly(ceil(end/2))];
theta   = 0;

% process input parameters
if length(varargin)==1
    % first argument contains all parameters
    var = varargin{1};
    center = var(1:2);
    theta = atan2(var(4), var(3));

elseif ~isempty(varargin)
    % parameters are given as separate arguments
    center = varargin{1};
    var = varargin{2};
    theta = atan2(var(2), var(1));
end

% transforms voxels according to square orientation
tra     = createTranslation(-center);
rot     = createRotation(-theta);
[x, y]  = transformPoint(x, y, rot*tra); %#ok<ASGLU>

% create image : simple threshold over 1 dimension
img = y > -1e-14;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [L,S] = PCP(M,lam,tol)
%
% [L,S] = PCP(M,lam,opts)

% This code solves the following model
%
% min_A { lam*||S(:)||_1 + ||L||_* }
% s.t. M = S+L

% where M is the data matrix, which will be decomposed into
% S sparse matrix S and low-rank matrix L.
%
% lam -- S small positive parameter

%% parameter setting
beta = .25/mean(abs(M(:))); %
maxit = 1000;

%% initialization
[m,n] = size(M);
S = zeros(m,n);
L = zeros(m,n);
Lambda = zeros(m,n); % the dual variable

% main
for iter = 1:maxit

    nrmLS = norm([S,L],'fro');
    % dS, dL record the change of S and L, only used for stopping criterion

    %% S - subproblem
    % S = argmin_A  lam*||S||_1 - <Lambda, S+L-M> + (beta/2) * ||S+L-M||.^2
    % Define element wise softshinkage operator as
    %     softshrink(z; gamma) = sign(z).* max(abs(z)-gamma, 0);
    % S has closed form solution: S=softshrink(Lambda/beta + M - L; lam/beta)
    % (see my slide page 42 Equation (66).
    X = Lambda / beta + M;
    Y = X - L;
    dS = S;
    S = sign(Y) .* max(abs(Y) - lam/beta, 0); % softshinkage operator
    dS = S - dS;

    %% L - subproblem
    % L = argmin_B ||L||_* -<Lambda, S+L-M> + + (beta/2) * ||S+L-M||.^2
    % L has closed form solution (singular value thresholding)
    % see my slide page 42, Equation (65).
    Y = X - S;
    dL = L;

    %[U,D,V] = svd(Y,'econ'); % use 'econ' is more efficient especially when Y is large
    [U,D,V] = svdecon(Y); % fastest

    VT=V';
    D = diag(D);
    ind = find(D > 1/beta);
    D = diag(D(ind) - 1/beta);
    L = U(:,ind) * D * VT(ind,:);
    dL = L - dL;

    %% stopping criterion
    RelChg = norm([dS,dL],'fro') / (1 + nrmLS);
    %fprintf('Iter %d, RelChg %4.2e \n',iter,RelChg);
    if RelChg < tol, break; end

    %% Update Lambda (dual variable)
    Lambda = Lambda - beta * (S + L - M);
end
end

author: andrewssobral

repo: lrslibrary


text:
function Prefs = ReadPreferences()
%ReadPreferences Reads Preferences text file
%
% ////////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED         ///
% ////////////////////////////////////////////
%
% INPUTS:
%   None
%
% OUTPUTS:
%  Prefs - Structure of Preferences
%
%VERSION:
%   1.0
%     - Tim Cox 20110525
%     - initial version.

Prefs = [];

%construct file name
[pathstr, name, ext] = fileparts(which('Taser'));
filename = sprintf('%s/TaserPreferences.txt',pathstr);

%determine if file exists
fid = fopen(filename,'r');

%return if fopen was unsucessful
if (fid < 3)
    return;
end

%parse file

%ImageOverviewFlag
tline = fgetl(fid);
Prefs.ImageOverview = str2double(tline(strfind(tline,'=')+1:end));
%SICD Location
tline = fgetl(fid);
Prefs.SICDLocation = strtrim(tline(strfind(tline,'=')+1:end));
%Max Resolution
tline = fgetl(fid);
Prefs.MaxResolution = str2double(tline(strfind(tline,'=')+1:end));
%Use Overview
tline = fgetl(fid);
Prefs.UseOverview = str2double(tline(strfind(tline,'=')+1:end));
%PHD NumPulses
tline = fgetl(fid);
Prefs.NumPulses = str2double(tline(strfind(tline,'=')+1:end));
%PHD ProcessingPulses
tline = fgetl(fid);
Prefs.ProcessingPulses = str2double(tline(strfind(tline,'=')+1:end));
%PHD dBCheck
tline = fgetl(fid);
Prefs.dBCheck = str2double(tline(strfind(tline,'=')+1:end));
%Complex Default Template Check
tline = fgetl(fid);
Prefs.ComplexDefaultCheck = str2double(tline(strfind(tline,'=')+1:end));
%Complex Template File
tline = fgetl(fid);
Prefs.ComplexTemplate = strtrim(tline(strfind(tline,'=')+1:end));
%PHD Default Template Check
tline = fgetl(fid);
Prefs.PHDDefaultCheck = str2double(tline(strfind(tline,'=')+1:end));
%PHD Template File
tline = fgetl(fid);
Prefs.PHDTemplate = strtrim(tline(strfind(tline,'=')+1:end));

fclose(fid);

% ////////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED         ///
% ////////////////////////////////////////////

author: ngageoint


text:
function test_readcompresseddata

% MEM 4gb
% WALLTIME 00:10:00
% DEPENDENCY inflate_file ft_read_data ft_read_header test_readcompresseddata

% test these data sets
datasets = {
  dccnpath('/home/common/matlab/fieldtrip/data/test/original/meg/ctf151/Subject01.ds.zip')
  dccnpath('/home/common/matlab/fieldtrip/data/test/original/eeg/brainvision/MischaCompressed.zip')
};

for k = 1:numel(datasets)
  cfg = [];
  cfg.continuous = 'yes';
  cfg.dataset = datasets{k};
  data = ft_preprocessing(cfg);
end

% and this mri
mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/test/latest/mri/dicom/dicomzipped.zip'));

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function a = moler1_inverse ( alpha, n )

%*****************************************************************************80
%
%% MOLER1_INVERSE returns the inverse of the MOLER1 matrix.
%
%  Formula:
%
%    if ( I = J )
%      A(I,J) =              min ( N-I, N-J ) * ALPHA^2 + 1
%    else
%      A(I,J) = (-1)^(I+J) * min ( N-I, N-J ) * ALPHA^2 + ALPHA
%
%  Example:
%
%    ALPHA = 2, N = 5
%
%     17 -14  10 -6  2
%    -14  13 -10  6 -2
%     10 -10   9 -6  2
%     -6   6  -6  5 -2
%      2  -2   2 -2  1
%
%  Properties:
%
%    The matrix is symmetric.
%
%    Successive elements of each diagonal decrease by ALPHA**2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the scalar that defines the inverse
%    Moler matrix.
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  v = zeros ( n, 1 );
  v(1) = 1.0;
  v(2) = - alpha;
  for i = 3 : n
    v(i) = - ( alpha - 1.0 ) * v(i-1);
  end

  for i = 1 : n
    for j = 1 : n
      if ( i <= j )
        a(i,j) = v(1+j-i:n+1-i)' * v(1    :n+1-j);
      else
        a(i,j) = v(1    :n+1-i)' * v(1+i-j:n+1-j);
      end
    end
  end

  return
end

repo: jburkardt-m


text:
function [U,S,V]=BL_SVD(A,u0,v0,k)

% Block lanczos method to compute SVD of matrix A with k steps;
% Suppose A is m*n, where m>n; A0=U0*S0*V0', where U0 is m*r, V0 is n*r;
% This is equivalent to compute block lanczos on C=[0,A;A',0];
% The initial lanczos block should set to be uv0=[U0;V0]/sqrt(2);
% where r is the rank of A0 and c is the slack positive interger variable;
%
% A - m * n matrix (required input)
% [u0; v0]/sqrt(2) - Initial subspace (required input)
% k - the number of blocks
%
% U S V - partial singular value decomposition
%
% Reference: Zhouchen Lin and Siming Wei, A Block Lanczos with Warm Start
% Technique for Accelerating Nuclear Norm Minimization Algorithms,
% arxiv: 1012.0365.
%
% Bug report: zclin2000@hotmail.com
%

if nargin < 3
    error('Too few arguments') ;
end

if nargin < 4
    k = 1;
elseif k == -1
    k = 1;
end

[m,n]=size(A);
mark=0;
if m<n % Let m>=n;
    mark=1;
    A=A';
    t=m;
    m=n;
    n=t;
    uv0=[v0; u0]/sqrt(2);
else
    uv0=[u0; v0]/sqrt(2);
end
r=size(uv0,2);
% Initialization;

% tridiagolization;
X(:,:,1)=uv0;
% AX=multiC(A,uv0,m,n);
AX = [A*uv0(m+1:m+n,:);A'*uv0(1:m,:)];
M(:,:,1)=uv0'*AX;
Q=X(:,:,1);
Bdiag=[];
Mdiag=M(:,:,1);
l=1;
stop=0;
tol=1.e-3;

while stop~=1
    if l==1
    R(:,:,l)=AX-uv0*M(:,:,1);
    else
    R(:,:,l)=AX-X(:,:,l)*M(:,:,l)-X(:,:,l-1)*B(:,:,l-1)';
    end
    R_max=max(max(abs(R(:,:,l))));

    if R_max<= tol
        stop=1;
    else
    [X(:,:,l+1),B(:,:,l)]=qr(R(:,:,l),0);
%     AX=multiC(A,X(:,:,l+1),m,n);
%     Y(:,:) = X(:,:,l+1);
    AX = [A*X(m+1:m+n,:,l+1);A'*X(1:m,:,l+1)];
    M(:,:,l+1)=X(:,:,l+1)'*AX;
    Q=[Q,X(:,:,l+1)];

    % Record the tridiagonal elments;
    Bdiag=blkdiag(Bdiag,B(:,:,l));
    Mdiag=blkdiag(Mdiag,M(:,:,l+1));
    l=l+1;
         if l>=k
             stop=1;
         end
    end
end

% Now T=Q'CQ is block tridiagonal;
B0=[zeros(r,l*r);Bdiag,zeros(l*r-r,r)];
T=Mdiag+B0+B0';
T=(T+T')/2;

% Compute EVD of T: T=Z*S*Z';
[Z,S]=eig(T);

% Therefore (QZ)'*C*(QZ)=S
add=0;
QZ=Q*Z(:,l*r:-1:l*r-r+1-add);
U_temp=QZ(1:m,:)*sqrt(2);
V_temp=QZ(m+1:m+n,:)*sqrt(2);
dS=diag(S);
S=diag(dS(l*r:-1:l*r+1-r-add));
if mark==1
    U=V_temp;
    V=U_temp;
else
    U=U_temp;
    V=V_temp;
end


text:
function d = mg_q1cd_supg(xy,ev,expe,eph,epw)
%mg_q1cd_supg  streamline diffusion matrix generator for GMG
%   d = mg_q1cd_supg(xy,ev,expe,eph,epw)
%   input
%           xy         vertex coordinate vector
%           ev         element mapping matrix
%           expe       element peclet numbers
%           eph        flow specific element lengths
%           epw        centroid evaluated wind
%   output
%           d          discrete streamline diffusion operator
%
%   IFISS function: DJS; 26 February 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

%   Analogous to femq1_cd_supg.m
x=xy(:,1); y=xy(:,2);
nvtx=length(x);
nel=length(ev(:,1));
%
% find the elements where streamline diffusion is active
acte = find(isfinite(expe)); nacte=length(acte);
%
% initialise global matrices
      d = sparse(nvtx,nvtx);
%
% set up 2x2 Gauss points
      gpt=1.0e0/sqrt(3.0e0);
      s(1) = -gpt;  t(1) = -gpt;
      s(2) =  gpt;  t(2) = -gpt;
      s(3) =  gpt;  t(3) =  gpt;
      s(4) = -gpt;  t(4) =  gpt;
%
% loop over active and inactive elements
%      for iact = 1:nacte
% ielem = acte(iact);
% inner loop over elements
        for ivtx = 1:4
        xl_v(:,ivtx) = x(ev(:,ivtx));
        yl_v(:,ivtx) = y(ev(:,ivtx));
		end
        de = zeros(nel,4,4);
% loop over 2x2 Gauss points
         for igpt = 1:4
         sigpt=s(igpt);
         tigpt=t(igpt);
%  evaluate derivatives etc
         [jac,invjac,phi,dphidx,dphidy] = deriv(sigpt,tigpt,xl_v,yl_v);
         [flowx,flowy] = gauss_transprt(sigpt,tigpt,xl_v,yl_v);
		 for j = 1:4
               for i = 1:4
    de(:,i,j) = de(:,i,j) + flowx(:).*dphidx(:,i).*flowx(:).*dphidx(:,j).*invjac(:);
    de(:,i,j) = de(:,i,j) + flowy(:).*dphidy(:,i).*flowx(:).*dphidx(:,j).*invjac(:);
    de(:,i,j) = de(:,i,j) + flowx(:).*dphidx(:,i).*flowy(:).*dphidy(:,j).*invjac(:);
    de(:,i,j) = de(:,i,j) + flowy(:).*dphidy(:,i).*flowy(:).*dphidy(:,j).*invjac(:);
               end
	    end
% end of Gauss point loop
         end
%
% scale with the appropriate parameter
      acte = find(isfinite(expe));
      factor = expe(acte); flow_h=eph(acte); flow_l2=epw(acte);
	  lpe =zeros(nel,1); lpe(acte)= factor.*(flow_h./flow_l2);
		 for j = 1:4
               for i = 1:4
               de(:,i,j) = lpe(:) .* de(:,i,j);
		   end
	   end
%
% perform assembly of global matrix  and source vector
      for krow=1:4
	  nrow=ev(:,krow);
          for kcol=1:4
		  ncol=ev(:,kcol);
          d = d + sparse(nrow,ncol,de(:,krow,kcol),nvtx,nvtx);
          end
      end
%
return


text:
function subpak_test07 ( )

%*****************************************************************************80
%
%% TEST07 tests BMI_ENGLISH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  test_num = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST07\n' );
  fprintf ( 1, '  BMI_ENGLISH computes the Body Mass Index\n' );
  fprintf ( 1, '  given body measurements in English Units.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      Weight         Height            BMI\n' );
  fprintf ( 1, '       (LB)     (FT           IN)\n' );
  fprintf ( 1, '\n' );

  seed = 123456789;

  for test = 1 : test_num

    b = 100.0;
    c = 250.0;

    [ w, seed ] = r8_uniform ( b, c, seed );

    b = 4.0;
    c = 6.75;

    [ h, seed ] = r8_uniform ( b, c, seed );

    h_ft = round ( h );
    h_in = round ( 12.0 * ( h - h_ft ) );

    bmi = bmi_english ( w, h_ft, h_in );

    fprintf ( 1, '  %8f  %8f  %8f  %8f\n', w, h_ft, h_in, bmi );

  end

  return
end


text:
function [Vs]=trisurfsmoothHC(F,V,IND_V,alp,bet,tol,n,disp_on)

%Humphreys Classes Smoothening

% alpha [0..1] influence of original/previous points 0-> full previous, 1-> full original
% beta  [0..1] e.g. > 0.5

if isempty(IND_V)
    [~,IND_V]=patchIND(F,V);
end

L=IND_V>0;
Xp=NaN(size(IND_V));  Yp=NaN(size(IND_V));  Zp=NaN(size(IND_V));

p=V; cc=0; i=1;
while cc==0
    q=p;
    Xp(L)=p(IND_V(L),1); Yp(L)=p(IND_V(L),2); Zp(L)=p(IND_V(L),3);
    p=[gnanmean(Xp,2) gnanmean(Yp,2) gnanmean(Zp,2)]; %Mean of neighbourhood, Laplacian operation
    SSQD_new=gnansum((p(:)-q(:)).^2);
    if i>1
        SSQD_ratio=SSQD_new./SSQD_old;
        if disp_on
            disp(['Iteration ',num2str(i),', SSQD ratio: ',num2str(SSQD_ratio)]);
        end
        if abs(1-SSQD_ratio)<=tol
            cc=1; %Convergence tolerance achieved
            if disp_on
                disp('Convergence tolerance on SSQD ratio reached!');
            end
        end
    else
        if disp_on
            disp(['Iteration ',num2str(i),', SSQD initial: ',num2str(SSQD_new)]);
        end
    end

    if i>=n
        cc=1;
        if disp_on
            disp('Maximum number of iteration reached!');
        end
    end
    SSQD_old=SSQD_new;

    b=p-((alp.*V)+((1-alp).*q)); %Difference at centre vertex
    Xp(L)=b(IND_V(L),1); Yp(L)=b(IND_V(L),2); Zp(L)=b(IND_V(L),3);
    c=[gnanmean(Xp,2) gnanmean(Yp,2) gnanmean(Zp,2)]; %Mean of difference at centre vertex of neighbourhood
    p=p-((bet.*b)+(1-bet).*(c));
    i=i+1;
end

Vs=p;

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% test_surface_param.m

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013-2014 University of Oxford
% Version: 0.6.3
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Open surface for 4 cardiac valves
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('data/valve-annula-points-009.mat');

% plot points
hold off
subplot(2, 2, 1)
plot3(x(:, 1), x(:, 2), x(:, 3), '*')
axis equal

% compute simple XY parametrization
uv0 = surface_param(x, 'xy');
d0 = dmatrix(uv0', uv0', 'euclidean');

% plot parametrization
subplot(2, 2, 2)
hold off
plot(uv0(:, 1), uv0(:, 2), '.')
axis equal
title('XY')

% compute PCA parametrization
uv = surface_param(x, 'pca');

% plot parametrization, aligned with the XY case
[~, uv] = procrustes(uv0, uv, 'Scaling', false);
subplot(2, 2, 3)
hold off
plot(uv(:, 1), uv(:, 2), '.')
axis equal
title('PCA')

% compute Isomap parametrization
param.type = 'isomap';
param.neigh = 'epsilon';
param.size = .0020;
param.d = dmatrix(x', x', 'euclidean');
uv = surface_param(x, param);

% plot parametrization
[~, uv] = procrustes(uv0, uv, 'Scaling', false);
subplot(2, 2, 4)
hold off
plot(uv(:, 1), uv(:, 2), '.')
axis equal
title('Isomap')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Open planar triangular mesh to compare isomap and MDS map
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% quadrangular mesh split into asymmetric triangles

% create rectangle with triangular mesh
[tri, x] = surface_tridomain('rect', 'step', 1/10, [0 0], [1 1/5*3]);
x(:, end+1) = 0;

% plot original mesh and parametrization
con = dmatrix_mesh(tri);
hold off
gplot(con, x(:, 1:2))

%% Isomap on a triangular mesh
% note that we need to compute the local distances on the mesh, and then
% apply Dijkstra. If we don't provide the distance matrix, by default
% surface_param() will compute distances between all pairs of vertices
clear param
param.type = 'isomap';
param.d = dmatrix_mesh(tri, x);
param.d = dijkstra(param.d, 1:length(param.d));
uv = surface_param(x, param);

% rigid registration to overlap the solution with the original data
[~, uv] = procrustes(x, uv, 'Scaling',false);

% add parametrization to the plot
hold on
gplot(con, uv, 'r')

%% Open classic MDSmap on the same triangular mesh
clear param
param.type = 'cmdsmap';
param.options.constraint_map = 0.4 * ones(size(x, 1), 1); % trigger warning
param.options.end_points = [1 2 3]; % trigger warning
param.options.nb_iter_max = 4; % trigger warning
param.tri = tri;
uv = surface_param(x, param);

% rigid registration to overlap the solution with the original data
[~, uv] = procrustes(x, uv, 'Scaling',false);

% add parametrization to the plot
hold on
gplot(con, uv, 'g')

%% check that MDSmap with Dijkstra is the same as Isomap
clear param
param.type = 'cmdsmap';
param.options.constraint_map = 0.4 * ones(1, size(x, 1)); % trigger warning
param.options.end_points = [1 2 3]; % trigger warning
param.options.nb_iter_max = 4; % trigger warning
[~, param.d] = dmatrix_mesh(tri, x, 'dijkstra');
uv = surface_param(x, param);

% rigid registration to overlap the solution with the original data
[~, uv] = procrustes(x, uv, 'Scaling',false);

% add parametrization to the plot
% it is expected that this perfectly overlaps the Isomap result computed
% above
hold on
gplot(con, uv, 'k')

%% lmdscale (MDSmap with local neighbourhood)
clear param
param.type = 'lmdscale';
param.tri = tri;
param.options.constraint_map = 0.4 * ones(size(x, 1), 1);
param.options2.MaxIter = 100;
param.options2.MaxInc = .01;
[uv, out] = surface_param(x, param);
out.stopCondition

% rigid registration to overlap the solution with the original data
[~, uv] = procrustes(x, uv, 'Scaling',false);

% add parametrization to the plot
% it is expected that this perfectly overlaps the Isomap result computed
% above
hold on
gplot(con, uv, 'c')

% plot error
figure
plot(out.err.maxinc)
hold on
plot(out.err.medinc)

% plot stress
hold off
plot(out.err.stress1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Closed surface from scattered point set
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% create a closed surface triangular mesh of a left ventricle

% seed the random generator so that we always get the same points
rng(0)

% load low resolution segmentation of LV, with detail removed
scimat = scimat_load('data/008-lv-resampled-0_31.mha');

% keep perimeter only
scimat.data = bwperim(scimat.data);

% coordinates of perimeter points
idx = find(scimat.data);
[r, c, s] = ind2sub(size(scimat.data), idx);
x = scimat_index2world([r c s], scimat);

% subsample boundary points
x = x(randi(size(x, 1), round(size(x, 1)/5), 1), :);

% order a bit the vertices, so that nearby vertices are close
x = sortrows(x);

% plot point set
subplot(2, 2, 1)
hold off
plot3(x(:, 1), x(:, 2), x(:, 3), '.')
axis equal

% mesh the points using an alpha-shape
[~, s] = alphavol(x, scimat.axis(1).spacing * 40);
tri = s.bnd;

% check that we don't have non-manifold vertices
idx = tri_find_nonmanifold_vertex(tri, x, scimat);
if (nnz(idx))
    error('Assertion fail: mesh has non-manifold vertices')
end

% remove points not connected to any triangle
[tri, x] = tri_squeeze(tri, x);

% plot mesh
subplot(2, 2, 2)
hold off
trisurf(tri, x(:, 1), x(:, 2), x(:, 3));
axis equal

% PCA-normalization of mesh
x = pca_normalize(x);

% plot mesh
subplot(2, 2, 3)
hold off
trisurf(tri, x(:, 1), x(:, 2), x(:, 3));
axis equal

% median radius of the points
sphrad = median(sqrt(sum(x.^2, 2)))

%% sphproj: simple projection on a sphere around the centroid
clear param
param.type = 'sphproj';
[latlon, out] = surface_param(x, param);
out.medrad

% rigid registration to overlap the solution with the original data
[xx, yy, zz] = sph2cart(latlon(:, 2), latlon(:, 1), sphrad);
[~, x1] = procrustes(x, [xx, yy, zz], 'Scaling',false);

% plot result
subplot(1, 2, 1)
hold off
trisurf(tri, x(:, 1), x(:, 2), x(:, 3));
axis equal
subplot(1, 2, 2)
hold off
trisurf(tri, x1(:, 1), x1(:, 2), x1(:, 3))
title('sphproj')
axis equal

%% CALD

% % using the SPHARM-MAT toolbox GUI
% faces = tri;
% vertices = x;
% save('/tmp/bar_obj.mat', 'faces', 'vertices')
% SPHARM_MAT
% load('/tmp/bar_CALD_smo.mat')

clear param
param.type = 'cald';
param.tri = tri;
param.options.MeshGridSize = 50;
param.options.MaxSPHARMDegree = 6;
param.options.Tolerance = 2;
param.options.Smoothing = 2;
param.options.Iteration = 100;
param.options.LocalIteration = 10;
latlon = surface_param(x, param);

% check for self-intersections
bad = cgal_check_self_intersect(param.tri, x);
nnz(bad>0)

% rigid registration to overlap the solution with the original data
[xx, yy, zz] = sph2cart(latlon(:, 2), latlon(:, 1), sphrad);
[~, x1] = procrustes(x, [xx, yy, zz], 'Scaling',false);

% plot parametrization
subplot(1, 2, 2)
hold off
trisurf(tri, x1(:, 1), x1(:, 2), x1(:, 3))
title('CALD')
axis equal

%% smdscale

% compute spherical MDSmap parametrization, no need the constrain the
% distance matrix more
clear param
param.type = 'smdscale';
param.tri = tri;
param.dmethod = 'fastmarching';
% param.options.constraint_map = .2 * sphrad * ones(size(x, 1), 1);
param.sphrad = sqrt(sum(cgal_trifacet_area(tri, x))/4/pi);
param.init = 'random';
param.options2.MaxIter = 15;
param.options2.MaxAlpha = .01 / 180 * pi;
[uv, out] = surface_param(x, param);
out.stopCondition
lat = uv(:, 1);
lon = uv(:, 2);

% plot parametrization
subplot(1, 2, 2)
hold off
[xx, yy, zz] = sph2cart(lon, lat, param.sphrad);
[~, x1] = procrustes(x, [xx, yy, zz], 'Scaling', false);
trisurf(param.tri, x1(:, 1), x1(:, 2), x1(:, 3));
title('smdscale')
axis equal

% check for self-intersections
bad = cgal_check_self_intersect(param.tri, x1);
nnz(bad>0)

% plot error
figure
subplot(2, 2, 1)
hold off
plot(out.err.rawstress)
title('Raw stress')
subplot(2, 2, 2)
hold off
plot(out.err.stress1)
title('Stress1')
subplot(2, 1, 2)
hold off
plot(out.err.maxalpha)
hold on
plot(out.err.medalpha)
title('Alpha')

% TODO:
% %% compare spherical Isomap and CALD
%
% % compute great circle distances for CALD
% dsphCALD = zeros(size(out.dsph));
% for I = 1:length(dsphCALD)
%     dsphCALD(:, I) = distance(latCALD(I), lonCALD(I), latCALD, lonCALD, ...
%         [out.sphrad 0], 'radians');
% end
%
% % plot distances to optimize vs spherical distances for CALD
% subplot(2, 1, 1)
% hold off
% plot(d(:), dsphCALD(:), '.')
% hold on
% plot([0 max(d(:))], [0 max(d(:))], 'r', 'LineWidth', 2)
% title('CALD')
% xlabel('Distance on the mesh')
% ylabel('Distance on the spherical parametrization')
% axis equal
%
%
% % plot boxplots with the distance ratios
% N = numel(d);
% subplot(1, 1, 1)
% hold off
% boxplot([d(:)-dsphCALD(:) ; d(:)-out.dsph(:)], ...
%     [ones(N, 1); 2*ones(N, 1)], 'labels', {'CALD', 'Sph Isomap'}, ...
%     'notch', 'on')
% title('Target distance - distance on spherical parametrization')

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function z = v_polynomial_zeros ( n )

%*****************************************************************************80
%
%% V_POLYNOMIAL_ZEROS returns zeroes of the Chebyshev polynomial V(n,x).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the polynomial.
%
%    Output, real Z(N), the zeroes.
%
  for i = 1 : n
    angle = ( 2 * n - 2 * i + 1 ) * pi / ( 2 * n + 1 );
    z(i) = cos ( angle );
  end

  return
end


text:
function varargout = errorMessage(varargin)
% ERRORMESSAGE M-file for errorMessage.fig
%      ERRORMESSAGE, by itself, creates a new ERRORMESSAGE or raises the existing
%      singleton*.
%
%      H = ERRORMESSAGE returns the handle to a new ERRORMESSAGE or the handle to
%      the existing singleton*.
%
%      ERRORMESSAGE('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in ERRORMESSAGE.M with the given input arguments.
%
%      ERRORMESSAGE('Property','Value',...) creates a new ERRORMESSAGE or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before errorMessage_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to errorMessage_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help errorMessage
% Last Modified by GUIDE v2.5 09-Jul-2013 23:40:34
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Doyen Sahoo
% Contributors: Bin LI, Steven C.H. Hoi
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @errorMessage_OpeningFcn, ...
                   'gui_OutputFcn',  @errorMessage_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before errorMessage is made visible.
function errorMessage_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to errorMessage (see VARARGIN)

% Choose default command line output for errorMessage
handles.output = hObject;
    message = varargin{1};
    set(handles.message,'String', message);
% Update handles structure
guidata(hObject, handles);

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

% --- Outputs from this function are returned to the command line.
function varargout = errorMessage_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
close all; clearvars; clc;
N = 256;
Phi = spx.dict.simple.dirac_hadamard_mtx(N);
figure;
imagesc(Phi);
colorbar;
export_fig images/demo_dirac_hadamard_1.png -r120 -nocrop;

mu1 = spx.dict.babel(Phi);
figure;
plot(mu1);
grid on;
export_fig images/demo_dirac_hadamard_babel.png -r120 -nocrop;

author: indigits

repo: sparse-plex


text:
function c = crossProduct3d(a,b)
%CROSSPRODUCT3D Vector cross product faster than inbuilt MATLAB cross.
%
%   C = crossProduct3d(A, B)
%   returns the cross product of the two 3D vectors A and B, that is:
%       C = A x B
%   A and B must be N-by-3 element vectors. If either A or B is a 1-by-3
%   row vector, the result C will have the size of the other input and will
%   be the  concatenation of each row's cross product.
%
%   Example
%     v1 = [2 0 0];
%     v2 = [0 3 0];
%     crossProduct3d(v1, v2)
%     ans =
%         0   0   6
%
%
%   Class support for inputs A,B:
%      float: double, single
%
%   See also DOT.

%   Sven Holcombe

% HISTORY
% 2017-11-24 rename from vectorCross3d to crossProduct3d

% size of inputs
sizeA = size(a);
sizeB = size(b);

% Initialise c to the size of a or b, whichever has more dimensions. If
% they have the same dimensions, initialise to the larger of the two
switch sign(numel(sizeA) - numel(sizeB))
    case 1
        c = zeros(sizeA);
    case -1
        c = zeros(sizeB);
    otherwise
        c = zeros(max(sizeA, sizeB));
end

c(:) = bsxfun(@times, a(:,[2 3 1],:), b(:,[3 1 2],:)) - ...
       bsxfun(@times, b(:,[2 3 1],:), a(:,[3 1 2],:));

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function K = gaussian_correlation ( s, t )

%*****************************************************************************80
%
%% GAUSSIAN_CORRELATION evaluates the Gaussian correlation function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Petter Abrahamsen,
%    A Review of Gaussian Random Fields and Correlation Functions,
%    Norwegian Computing Center, 1997.
%
%  Parameters:
%
%    Input, real S(*), T(*), pairs of argument values.
%
%    Output, real K(*), the correlation function values
%
  K = exp ( - ( s - t ).^2 );

  return
end


text:
function linpack_d_test31 ( )

%*****************************************************************************80
%
%% TEST31 tests DTRSL.
%
%  Discussion:
%
%    DTRSL solves triangular linear systems.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 5;
  lda = n;

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST31\n' );
  fprintf ( 1, '  For a triangular matrix,\n' );
  fprintf ( 1, '  DTRSL solves a linear system.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Lower triangular matrix A.
%
  [ a, seed ] = r8mat_uniform_01 ( n, n, seed );

  for i = 1 : n
    for j = i+1 : n
      a(i,j) = 0.0;
    end
  end

  for i = 1 : n
    x(i,1) = i;
  end

  b(1:n,1) = a(1:n,1:n) * x(1:n,1);

  fprintf ( 1, '\n' );
  fprintf ( 1, '  For a lower triangular matrix A,\n' );
  fprintf ( 1, '  solve A * x = b\n' );

  job = 00;

  [ b, info ] = dtrsl ( a, lda, n, b, job );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The solution (should be 1,2,3,4,5):\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %6d  %14f\n', i, b(i,1) );
  end

  b(1:n,1) = a(1:n,1:n)' * x(1:n,1);

  fprintf ( 1, '\n' );
  fprintf ( 1, '  For a lower triangular matrix A,\n' );
  fprintf ( 1, '  solve A'' * x = b\n' );

  job = 10;

  [ b, info ] = dtrsl ( a, lda, n, b, job );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The solution (should be 1,2,3,4,5):\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %6d  %14f\n', i, b(i,1) );
  end
%
%  Upper triangular matrix A.
%
  [ a, seed ] = r8mat_uniform_01 ( n, n, seed );

  for i = 1 : n
    for j = 1 : i - 1
      a(i,j) = 0.0;
    end
  end

  for i = 1 : n
    x(i,1) = i;
  end

  b(1:n,1) = a(1:n,1:n) * x(1:n,1);

  fprintf ( 1, '\n' );
  fprintf ( 1, '  For an upper triangular matrix A,\n' );
  fprintf ( 1, '  solve A * x = b\n' );

  job = 01;

  [ b, info ] = dtrsl ( a, lda, n, b, job );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The solution (should be 1,2,3,4,5):\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %6d  %14f\n', i, b(i,1) );
  end

  b(1:n,1) = a(1:n,1:n)' * x(1:n,1);

  fprintf ( 1, '\n' );
  fprintf ( 1, '  For an upper triangular matrix A,\n' );
  fprintf ( 1, '  solve A'' * x = b\n' );

  job = 11;

  [ b, info ] = dtrsl ( a, lda, n, b, job );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The solution (should be 1,2,3,4,5):\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %6d  %14f\n', i, b(i,1) );
  end

  return
end


text:
% Example_4_PRVI: Mesh the west North Atlantic Ocean with high resolution
% around Puerto Rico and US Virgin Islands.
% User beware: this is a very intensive mesh generation example, takes
% around 1 hr.
clc; clearvars

addpath('..')
addpath(genpath('../utilities/'))
addpath(genpath('../datasets/'))
addpath(genpath('../m_map/'))

%% The constant parameters for all domains
wl        = 30;         % elements to resolve M2 wavelength
dt        = 0;          % use automatic timestep
grade     = 0.25; 		% mesh grade in decimal percent.
R         = -5;    		% number of elements to resolve feature width.
slp       = 15;         % 2*pi/number of elements to resolve slope
fl        = -50;        % use filter equal to Rossby radius divided by 50

%%  For relatively coarse resolution west North Atlantic Ocean
bbox = [-100 -53   	    % lon_min lon_max
         5    52.5];    % lat_min lat_max
min_el = 1000;  		% minimum resolution in meters.
max_el = 10e3;       	% maximum resolution in meters.

coastline = 'GSHHS_f_L1';
dem       = 'SRTM15+.nc';
gdat{1} = geodata('shp',coastline,...
                  'dem',dem,...
                  'bbox',bbox,...
                  'h0',min_el);

fh{1} = edgefx('geodata',gdat{1},...
               'fs',R,...
               'wl',wl,...
               'slp',slp,'fl',fl,...
               'max_el',max_el,...
               'dt',dt,...
               'g',grade);

%% For High Resolution around Puerto Rico and US Virgin Islands
coastline = ["pr_1s_0m_contour","usvi_0m_contour","sj_0contour_closed"];

dems      = ["pr_1s.nc","usvi_1_mhw_2014.nc", "san_juan_19_prvd02_2015.nc"];

for ii = 1:length(dems)
    % use same parameters as coarse mesh, just change min_el
    if ii == length(dems)
        min_el    = 10;  % minimum resolution in meters.
    else
        min_el    = 30;  % minimum resolution in meters.
    end
    % bbox is taken automatically from the DEM
    gdat{ii+1} = geodata('shp',coastline{ii},...
                         'dem',dems{ii},...
                         'h0',min_el);
    fh{ii+1} = edgefx('geodata',gdat{ii+1},...
                      'fs',R,...
                      'wl',wl,...
                      'slp',slp,'fl',fl,...
                      'max_el',max_el,...
                      'dt',dt,...
                      'g',grade);
end

%% Pass your edgefx class objects along with some meshing options
%% and build the mesh...
% (note that the nested edgefxs will be smoothed together with this call)
mshopts = meshgen('ef',fh,'bou',gdat,'plot_on',1,'itmax',50,...
                  'proj','lam');

% now build the mesh with your options and the edge function.
mshopts = mshopts.build;

% Get out the msh class from meshgen
m = mshopts.grd;

%% Interpolate on the bathy and gradients (automatically loops over all data)
m = interp(m,gdat,'mindepth',1); % interpolate bathy to the mesh with minimum depth of 1 m

%% Make the nodestrings
m = make_bc(m,'auto',gdat{1}); % make the nodestring boundary conditions

%% Plot and save the msh class object/write to fort.14
plot(m,'type','bd'); % plot triangulation with boundary conditions
plot(m,'type','b');  % plot the bathy
plot(m,'type','reso','colormap',[10 0 10e3]) % plot the resolution
% Save as a msh class
save('PRVI_msh.mat','m');
% Write an ADCIRC fort.14 compliant file to disk.
write(m,'PRVI_mesh')

author: CHLNDDEV


text:
function [doseHandle,cMap,window] = matRad_plotDoseSlice3D(axesHandle,ct,doseCube,plane,slice,threshold,alpha,cMap,window)
% matRad function that generates a dose plot of a selected slice in 3D view
%
% call
%   [doseHandle,cMap,window] = matRad_plotDose3DSlice(axesHandle, doseCube,plane,slice,threshold)
%   [doseHandle,cMap,window] = matRad_plotDose3DSlice(axesHandle, doseCube,plane,slice,threshold,alpha)
%   [doseHandle,cMap,window] = matRad_plotDose3DSlice(axesHandle, doseCube,plane,slice,threshold,cMap)
%   [doseHandle,cMap,window] = matRad_plotDose3DSlice(axesHandle, doseCube,plane,slice,threshold,window)
%   [doseHandle,cMap,window] = matRad_plotDose3DSlice(axesHandle, doseCube,plane,slice,threshold,alpha,cMap)
%   [doseHandle,cMap,window] = matRad_plotDose3DSlice(axesHandle, doseCube,plane,slice,threshold,alpha,window)
%   [doseHandle,cMap,window] = matRad_plotDose3DSlice(axesHandle, doseCube,plane,slice,threshold,cMap,window)
%   [doseHandle,cMap,window] = matRad_plotDose3DSlice(axesHandle, doseCube,plane,slice,threshold,alpha,cMap,window)
%
% input
%   axesHandle  handle to axes the slice should be displayed in
%   ct          matRad CT struct which contains resolution
%   doseCube    3D array of the dose to select the slice from
%   plane       plane view (coronal=1,sagittal=2,axial=3)
%   slice       slice in the selected plane of the 3D cube
%   threshold   threshold above which the dose shall be displayed
%               for negative values (i.e. difference maps), also the values
%               smaller than the negative threshold will be displayed
%               if empty, no threshold will be applied
%   alpha       optional argument defining the alpha value, default is 0.6.
%               To use the default when providing a custom culormap, put in
%               an empty array by [].
%   cMap        optional argument defining the colormap, default is jet
%               if you want to use the default map with the window argument
%               you can use an empty array []
%   window      optional argument defining the displayed range. default is
%               [min(doseCube(:)) max(doseCube(:))]
%
% output
%   doseHandle: handle of the plotted dose axes
%   cMap        used colormap (same as input if set)
%   window      used window (same as input if set)
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

matRad_cfg = MatRad_Config.instance();

%Use default colormap?
if nargin < 8 || isempty(cMap)
    cMap = jet(64);
end
if nargin < 7 || isempty(alpha)
    alpha = 0.6;
end
if nargin < 9 || isempty(window)
    window = [min(doseCube(:)) max(doseCube(:))];
end

cMapScale = size(cMap,1) - 1;
maxDose = max(doseCube(:));

%Create the coordinates
coords{1} = ct.resolution.x * (1:ct.cubeDim(2));
coords{2} = ct.resolution.y * (1:ct.cubeDim(1));
coords{3} = ct.resolution.z * (1:ct.cubeDim(3));

if plane == 1  % Coronal plane
    [xMesh,zMesh] = meshgrid(coords{2},coords{3});
    yMesh = slice*ct.resolution.x*ones(size(xMesh));
    %dose_slice = uint8(cMapScale*(squeeze(doseCube(slice,:,:)) - window(1))/(window(2)-window(1)));
    doseSlice = permute(squeeze(doseCube(slice,:,:)),[2 1]);
elseif plane == 2 % sagittal plane
    [yMesh,zMesh] = meshgrid(coords{1},coords{3});
    xMesh = slice*ct.resolution.y*ones(size(yMesh));
    %dose_slice = uint8(cMapScale*(squeeze(doseCube(:,slice,:)) - window(1))/(window(2)-window(1)));
    dose_slice = permute(squeeze(doseCube(:,slice,:)),[2 1]);
elseif plane == 3 % Axial plane
    [xMesh,yMesh] = meshgrid(coords{2},coords{1});
    zMesh = slice*ct.resolution.z*ones(size(xMesh));
    %dose_slice = uint8(cMapScale*(squeeze(doseCube(:,:,slice)) - window(1))/(window(2)-window(1)));
    dose_slice = squeeze(doseCube(:,:,slice));
end

if ~isempty(threshold)
    dose_mask = alpha * (dose_slice < window(2) & dose_slice > window(1) & dose_slice > threshold*maxDose);
else
    dose_mask = alpha * (dose_slice < window(2) & dose_slice > window(1));
end

dose_slice = uint8(cMapScale* (dose_slice - window(1))/(window(2)-window(1)));

%This circumenvents a bug in Octave when the index in the image hase the maximum value of uint8
if matRad_cfg.isOctave
	dose_slice(dose_slice == 255) = 254;
end

dose_rgb = ind2rgb(dose_slice,cMap);
% slice plot with surface(...), colormapping can be done by texture
% mapping, this is why we use surface instead of slice
doseHandle = surface('XData',xMesh, 'YData',yMesh, 'ZData',zMesh,'AlphaData',dose_mask, ...
        'CData',dose_rgb, 'CDataMapping','direct', ...
        'EdgeColor','none', 'FaceColor','texturemap', 'BackFaceLighting','unlit','FaceLighting','flat','FaceAlpha','texture','Parent',axesHandle);

end

repo: matRad


text:
function F = in_fread_edf(sFile, sfid, SamplesBounds, ChannelsRange)
% IN_FREAD_EDF:  Read a block of recordings from a EDF/BDF file
%
% USAGE:  F = in_fread_edf(sFile, sfid, SamplesBounds, ChannelsRange)
%         F = in_fread_edf(sFile, sfid, SamplesBounds)               : Read all channels
%         F = in_fread_edf(sFile, sfid)                              : Read all channels, all the times

% @=============================================================================
% 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, 2012-2023

%% ===== PARSE INPUTS =====
nChannels      = sFile.header.nsignal;
iChanAnnot     = find(strcmpi({sFile.header.signal.label}, 'EDF Annotations'));
iBadChan       = find(sFile.channelflag == -1);
iChanSignal    = setdiff(1:nChannels, iChanAnnot);
iChanWrongRate = find([sFile.header.signal(iChanSignal).sfreq] ~= max([sFile.header.signal(setdiff(iChanSignal,iBadChan)).sfreq]));
iChanSkip      = union(iChanAnnot, iChanWrongRate);
% By default: Read all the channels
if (nargin < 4) || isempty(ChannelsRange)
    ChannelsRange = [1, nChannels];
end
% By default: Read all the file, at the main file frequency
if (nargin < 3) || isempty(SamplesBounds)
    SamplesBounds = round(sFile.prop.times .* sFile.prop.sfreq);
end

% Block of times/channels to extract
nReadChannels = double(ChannelsRange(2) - ChannelsRange(1) + 1);
% Read annotations instead of real data ?
isAnnotOnly = ~isempty(iChanAnnot) && (ChannelsRange(1) == ChannelsRange(2)) && ismember(ChannelsRange(1), iChanAnnot);
isBdfStatus = strcmpi(sFile.format, 'EEG-BDF') && (ChannelsRange(1) == ChannelsRange(2)) && strcmpi(sFile.header.signal(ChannelsRange(1)).label, 'Status');
% Data channels to read
if isAnnotOnly
    % Read only on annotation channel
    iChanF = 1;
else
    % Remove all the annotation channels from the list of channels to read
    iChanF = setdiff(ChannelsRange(1):ChannelsRange(2), iChanSkip) - ChannelsRange(1) + 1;
%     if any(diff(iChanF) ~= 1)
%         error('All the data channels to read from the file must be contiguous (EDF Annotation channels must be at the end of the list).');
%     end
    if any(diff(iChanF) ~= 1)
        iChanLast = find(diff(iChanF) ~= 1, 1);
        iChanF = iChanF(1:iChanLast);
    else
        iChanLast = length(iChanF);
    end
    ChannelsRange = [iChanF(1), iChanF(iChanLast)] + ChannelsRange(1) - 1;
end
% Cannot read channels with different sampling rates at the same time
if (ChannelsRange(1) ~= ChannelsRange(2))
    allFreq = [sFile.header.signal(ChannelsRange(1):ChannelsRange(2)).nsamples];
    if any(allFreq ~= allFreq(1))
        error(['Cannot read channels with mixed sampling rates at the same time.' 10 ...
               'Mark as bad channels with different sampling rates than EEG.' 10 ...
               '(right-click on data file > Good/bad channels > Edit good/bad channels)']);
    end
end
% Number of time samples to read
nTimes = round(sFile.header.reclen * sFile.header.signal(ChannelsRange(1)).sfreq);
iTimes = SamplesBounds(1):SamplesBounds(2);

%% ===== READ ALL NEEDED EPOCHS =====
% Detect which epochs are necessary for the range of data selected
epochRange = floor(SamplesBounds ./ nTimes);
epochsToRead = epochRange(1) : epochRange(2);
% Initialize block of data to read
if isAnnotOnly
    F = zeros(nReadChannels, 2 * length(iTimes));
else
    F = zeros(nReadChannels, length(iTimes));
end
% Marker that we increment when we add data to F
iF = 1;
% Read all the needed epochs
for i = 1:length(epochsToRead)
    % Find the samples to read from this epoch
    BoundsEpoch = nTimes * epochsToRead(i) + [0, nTimes-1];
    BoundsRead  = [max(BoundsEpoch(1), iTimes(1)), ...
                   min(BoundsEpoch(2), iTimes(end))];
    iTimeRead = BoundsRead(1):BoundsRead(2);
    % Convert this samples into indices in this very epoch
    iTimeRead = iTimeRead - nTimes * epochsToRead(i);
    % New indices to read
    if isAnnotOnly
        iNewF = iF:(iF + 2*length(iTimeRead) - 1);
    else
        iNewF = iF:(iF + length(iTimeRead) - 1);
    end
    % Read epoch (full or partial)
    F(iChanF,iNewF) = edf_read_epoch(sFile, sfid, epochsToRead(i), iTimeRead, ChannelsRange, isAnnotOnly, isBdfStatus);
    % Increment marker
    iF = iF + length(iTimeRead);
end

end

%% ===== READ ONE EPOCH =====
function F = edf_read_epoch(sFile, sfid, iEpoch, iTimes, ChannelsRange, isAnnotOnly, isBdfStatus)
    % ===== COMPUTE OFFSETS =====
    nTimes        = sFile.header.reclen * [sFile.header.signal.sfreq];
    nReadTimes    = length(iTimes);
    nReadChannels = double(ChannelsRange(2) - ChannelsRange(1) + 1);
    iChannels     = ChannelsRange(1):ChannelsRange(2);
    % Check that all the channels selected have the same freq rate
    if any(nTimes(iChannels) ~= nTimes(iChannels(1)))
        error('Cannot read at the same signals with different sampling frequency.');
    end
    % Size of one value
    if strcmpi(sFile.format, 'EEG-BDF')
        % BDF: int24 => 3 bytes
        bytesPerVal = 3;
        % Reading status or regular channel
        if isBdfStatus
            dataClass = 'ubit24';
        else
            dataClass = 'bit24';
        end
    else
        % EDF: int16 => 2 bytes
        bytesPerVal = 2;
        dataClass = 'int16';
        isBdfStatus = 0;
    end
    % Offset of the beginning of the recordings in the file
    offsetHeader    = round(sFile.header.hdrlen);
    % Offset of epoch
    offsetEpoch     = round(iEpoch * sum(nTimes) * bytesPerVal);
    % Channel offset
    offsetChannel   = round(sum(nTimes(1:ChannelsRange(1)-1)) * bytesPerVal);
    % Time offset at the beginning and end of each channel block
    offsetTimeStart = round(iTimes(1) * bytesPerVal);
    offsetTimeEnd   = round((nTimes(ChannelsRange(1)) - iTimes(end) - 1) * bytesPerVal);
    % ALL THE "ROUND" CALLS WERE ADDED AFTER DISCOVERING THAT THERE WERE SOMETIMES ROUNDING ERRORS IN THE MULTIPLICATIONS

    % Where to start reading in the file ?
    % => After the header, the number of skipped epochs, channels and time samples
    offsetStart = offsetHeader + offsetEpoch + offsetChannel + offsetTimeStart;
    % Number of time samples to skip after each channel
    offsetSkip = offsetTimeStart + offsetTimeEnd;

    % ===== READ DATA BLOCK =====
    % Position file at the beginning of the trial
    fseek(sfid, offsetStart, 'bof');
    % Read annotation data (char)
    if isAnnotOnly
        dataClass = 'uint8';
        nReadTimes = bytesPerVal * nReadTimes;  % 1 byte instead of 2
    end
    % Read trial data
    % => WARNING: CALL TO FREAD WITH SKIP=0 DOES NOT WORK PROPERLY
    if (offsetSkip == 0)
        F = fread(sfid, [nReadTimes, nReadChannels], dataClass)';
    elseif (bytesPerVal == 2)
        precision = sprintf('%d*%s', nReadTimes, dataClass);
        F = fread(sfid, [nReadTimes, nReadChannels], precision, offsetSkip)';
    % => WARNING: READING USING ubit24 SOMETIMES DOESNT WORK => DOING IT MANUALLY
    elseif (bytesPerVal == 3)
        % Reading each bit independently
        precision = sprintf('%d*%s', 3*nReadTimes, 'uint8');
        F = fread(sfid, [3*nReadTimes, nReadChannels], precision, offsetSkip)';
        % Grouping the 3 bits together
        F = F(:,1:3:end) + F(:,2:3:end)*256 + F(:,3:3:end)*256*256;
        % 2-Complement (negative value indicated by most significant bit)
        if strcmpi(dataClass, 'bit24')
            iNeg = (F >= 256*256*128);
            F(iNeg) = F(iNeg) - 256*256*256;
        end
    end
    % Check that data block was fully read
    if (numel(F) < nReadTimes * nReadChannels)
        % Number of full time samples that were read
        nTimeTrunc = max(0, floor(numel(F) / nReadChannels) - 1);
        % Error message
        disp(sprintf('EDF> ERROR: File is truncated (%d time samples were read instead of %d). Padding with zeros...', nTimeTrunc, nReadTimes));
        % Pad data with zeros
        Ftmp = zeros(nReadTimes, nReadChannels);
        F = F';
        Ftmp(1:numel(F)) = F(:);
        F = Ftmp';
    end

    % Processing for BDF status file
    if isBdfStatus
        % Mask to keep only the first 16 bits (Triggers bits)
        % Bit 16    : High when new Epoch is started
        % Bit 17-19 : Speed bits 0 1 2
        % Bit 20 	: High when CMS is within range
        % Bit 21 	: Speed bit 3
        % Bit 22 	: High when battery is low
        % Bit 23    : High if ActiveTwo MK2
        F = bitand(F, hex2dec('00FFFF'));
    % Processing for real data
    elseif ~isAnnotOnly
        % Convert to double
        F = double(F);
        % Apply gains
        F = bst_bsxfun(@rdivide, F, [sFile.header.signal(iChannels).gain]');
% IN THEORY: THIS OFFSET SECTION IS USEFUL, BUT IN PRACTICE IT LOOKS LIKE THE VALUES IN ALL THE FILES ARE CENTERED ON ZERO
%         % Add offset
%         if isfield(sFile.header.signal, 'offset') && ~isempty(sFile.header.signal(1).offset)
%             % ...
%         end
    end
end

author: brainstorm-tools


text:
function noise = heavisideNoiseExpandParam(noise, params)

% HEAVISIDENOISEEXPANDPARAM Expand heaviside noise structure from param vector.

% IVM

noise.eta = 0.5*sigmoidBound(params(1));
noise.bias = params(2:end);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [b_value_ma, a_value_ma, b_error_ma] = ma(mag,cumnum)

    % compute the slope, intercept, and their errors using
    % "Major Axis" regression.
    % the inputs are two column vectors.
    % e.g. [b,a,e] = ma(X,Y) were X and Y are column vectors of the same length.
    % outputs are the b and a values and the standard error.

    report_this_filefun();

    X = [mag cumnum];
    m = mean(X);
    ssx = (X(:,1)-m(1))'*(X(:,1)-m(1));
    ssy = (X(:,2)-m(2))'*(X(:,2)-m(2));
    ssxy = (X(:,1)-m(1))'*(X(:,2)-m(2));

    % MAJOR AXIS REGRESSION
    [V,D] = eig(cov(X));
    % the slope of the line in the direction of max. variance
    b_value_ma = V(2,2)/V(1,2);
    a_value_ma = m(2) - (b_value_ma * m(1));

    % the standard error in the slopes
    r = ssxy/(sqrt(ssx.*ssy));
    b_error_ma = abs(b_value_ma)*sqrt((1-r^2)/length(X));
end

author: CelsoReyes


text:
function node_h = h_coef ( node_num, node_xy )

%*****************************************************************************80
%
%% H_COEF evaluates the coefficient H(X,Y) of DEL U in the Poisson equation.
%
%  Discussion:
%
%    The equation is
%
%      - Del H(X,Y) DEL U + K(X,Y) * U = RHS(X,Y)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM),
%    the coordinates of the points.
%
%    Output, real NODE_H(NODE_NUM,1),
%    the value of the coefficient of DEL U .
%
  node_h(1:node_num,1) = 1.0;

  return
end


text:
function [models, info] = slproglearn(source, getter, learnfunctor, varargin)
%SLPROGLEARN Performs Progressive Learning from sample source
%
% $ Syntax $
%   - models = slproglearn(source, getter, learnfunctor, ...)
%   - [models, info] = slproglearn(source, getter, learnfunctor, ...)
%
% $ Arguments $
%   - source:       the source from which the data are fetched
%   - getter:       the functor to fetch data from source, in the form:
%                   data= f(source, ...)
%   - learnfunctor: the functor to construct/update models from data,
%                   it is given in the following form:
%                   models = f(models, data, ...)
%                   If the learning is recorded, it is like:
%                   [models, record] = f(models, data, ...)
%                   On initial construction, the input models is [].
%   - models:       The constructed models
%
% $ Description $
%   - models = slproglearn(source, getter, learnfunctor, ...)
%     construct models based on source, from which the samples are
%     fetched. It is assumed that the source can continuously offer
%     infinite number of samples. The construction is controlled by
%     the following properties.
%     \*
%     \t   Table. Properties of Progressive Learning
%     \h       name        &            description             \\
%          'isrecorded'    &  Whether each iteration is recorded
%                             (default = false)                 \\
%          'gtuner'        &  The functor to tune the getter, in the form:
%                             getter = f(getter, models, ...)
%                             (default = {})                      \\
%          'cmpfunctor'    &  The functor to compare two set of models
%                             and judge whether discrepancies meet the
%                             criteria of convergence.
%                             It is given in the following form:
%                             isconverged = f(models1, models2, ...)
%                             This properties must be specified.     \\
%          'iter'          &  The iteration control properties for
%                             sliterproc, default = {}.
%          'verbose'       &  whether to show progress information
%                             (default = true)                       \\
%          'initmodels'    &  The initial models, (default = [])
%     \*
%
%   - [models, info] = slproglearn(source, getter, learnfunctor, ...) also
%     returns the information of iterative process.
%
% $ History $
%   - Created by Dahua Lin, on Aug 31, 2006
%

%% Parse and verify input arguments

if nargin < 3
    raise_lackinput('slproglearn', 3);
end

opts.isrecorded = false;
opts.gtuner = [];
opts.cmpfunctor = [];
opts.iter = {};
opts.verbose = true;
opts.initmodels = [];
opts = slparseprops(opts, varargin{:});

if isempty(opts.cmpfunctor)
    error('sltoolbox:invalidarg', ...
        'You should specify a models comparison functor');
end

%% Main Skeleton

slsharedisp_attach('slproglearn', 'show', opts.verbose);

slsharedisp(opts, 'Progressive Learning from source');

objects = {source, opts.initmodels, getter};
iterfunctor = {@proglearn_iter, learnfunctor, opts};
cmpfunctor = {@proglearn_cmp, opts};
if nargout < 2
    objects = ...
        sliterproc(objects, iterfunctor, cmpfunctor, opts.isrecorded, opts.iter{:});
else
    [objects, info] = ...
        sliterproc(objects, iterfunctor, cmpfunctor, opts.isrecorded, opts.iter{:});
end

models = objects{2};

slsharedisp_detach();

%% Core Iteration function

function varargout = proglearn_iter(objects, learnfunctor, opts)
% objects = {source, models, getter}

% take input
source = objects{1};
models = objects{2};
getter = objects{3};

% fetch data
data = slevalfunctor(getter, source);

% learn models
if ~opts.isrecorded
    models = slevalfunctor(learnfunctor, models, data);
else
    [models, rec] = slevalfunctor(learnfunctor, models, data);
end

% tune getter
if ~isempty(opts.gtuner)
    getter = slevalfunctor(opts.gtuner, getter, models);
end

% make output
objects = {source, models, getter};
if ~opts.isrecorded
    varargout = {objects};
else
    varargout = {objects, rec};
end

function isconverged = proglearn_cmp(objects_prev, objects, opts)

models_prev = objects_prev{2};
models = objects{2};

isconverged = slevalfunctor(opts.cmpfunctor, models_prev, models);

author: lmthang


text:
function f=involute(f,dim);
%INVOLUTE  Involution
%   Usage: finv=involute(f);
%          finv=involute(f,dim);
%
%   `involute(f)` will return the involution of *f*.
%
%   `involute(f,dim)` will return the involution of *f* along dimension *dim*.
%   This can for instance be used to calculate the 2D involution::
%
%     f=involute(f,1);
%     f=involute(f,2);
%
%   The involution *finv* of *f* is given by::
%
%     finv(l+1)=conj(f(mod(-l,L)+1));
%
%   for $l=0,\ldots,L-1$.
%
%   The relation between conjugation, Fourier transformation and involution
%   is expressed by::
%
%     conj(dft(f)) == dft(involute(f))
%
%   for all signals *f*. The inverse discrete Fourier transform can be
%   expressed by::
%
%     idft(f) == conj(involute(dft(f)));
%
%   See also:  dft, pconv

% Assert correct input.

%   AUTHOR : Peter L. Søndergaard
%   TESTING: TEST_INVOLUTE
%   REFERENCE: OK

complainif_argnonotinrange(nargin,1,2,mfilename);

if nargin==1
  dim=[];
end;

L=[];
[f,L,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,L,dim,'INVOLUTE');

% This is where the calculation is performed.
% The reshape(...,size(f) ensures that f will keep its
% original shape if it is multidimensional.
f=reshape(conj([f(1,:); ...
	  flipud(f(2:L,:))]),size(f));

f=assert_sigreshape_post(f,dim,permutedsize,order);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function bigclu(var1)
    %bigclu.m                     A.Allmann
    %function to select only clusters with spezial number values
    %
    % Last modification 8/95
   global mess te1  sys
    global tmp1 tmp2 cluslength
    global freq_field1 freq_field2 close_button go_button
    global backbgevent original equi bgevent backequi clu newclcat backcat
    global equi_button bg_button
    global plot1_h plot2_h file1 clust

    if var1==1
        figure_w_normalized_uicontrolunits(mess)
        clf
        set(gca,'visible','off');
        set(gcf,'visible','off');
        set(gcf,'Name','Number Selection');
        cltemp=cluslength(equi(:,10));
        tmp2=min(cltemp);
        freq_field1= uicontrol('Style','edit',...
            'Position',[.70 .60 .17 .10],...
            'Units','normalized','String',num2str(tmp2),...
            'Callback','tmp2=str2double(get(freq_field1,''String'')); set(freq_field1,''String'',num2str(tmp2));');

        tmp1=max(cltemp);
        freq_field2=uicontrol('Style','edit',...
            'Position',[.70 .40 .17 .10],...
            'Units','normalized','String',num2str(tmp1),...
            'Callback','tmp1=str2double(get(freq_field2,''String'')); set(freq_field2,''String'',num2str(tmp1));');

        close_button=uicontrol('Style','Pushbutton',...
            'Position', [.60 .05 .15 .15 ],...
            'Units','normalized','Callback','welcome;done','String','Cancel');
        go_button=uicontrol('Style','Pushbutton',...
            'Position',[.25 .05 .15 .15 ],...
            'Units','normalized',...
            'Callback','welcome;done;bigclu(3);',...
            'String','Go');

        txt1 = text(...
            'Color',[0 0 0 ],...
            'EraseMode','normal',...
            'Position',[0. 0.65 0 ],...
            'Rotation',0 ,...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold' ,...
            'String','Minimum Events in Cluster:');

        txt2 = text(...
            'Color',[0 0 0 ],...
            'EraseMode','normal',...
            'Position',[0. 0.40 0 ],...
            'Rotation',0 ,...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold' ,...
            'String','Maximum Events in Cluster:');

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

    elseif var1==3
        figure_w_normalized_uicontrolunits(clu);
        if isempty(newclcat)  &&  isempty(backcat)   %no selection before
            backequi=equi;
            backbgevent=bgevent;
        end
        tt= find(cluslength>=tmp2 & cluslength<=tmp1);
        for j= tt
            tt1=find(equi(:,10)==j);
            if isempty(tt1)
                tt1=0;
            end
            tt2(j)=tt1;
        end
        tmp=find(tt2);
        equi=backequi(tmp,:);
        bgevent=backbgevent(tmp,:);

        set(equi_button,'value',1)
        st1=get(equi_button,'Callback');
        eval(st1);
        pause(2);
        tmpcat=clust(:,tmp);
        newclcat=original(tmpcat(find(clust(:,tmp))),:);
        plot1_h=[];plot2_h=[];
        cluoverl(7);

        strib=[' Polygon of  ' file1];
        hold on
        title(strib,'FontWeight','bold',...
            'FontSize',ZmapGlobal.Data.fontsz.l,'Color','r')

    end


text:
%--- help for generic/log_prior_density ---
%
%  Computes the probability density function of the prior corresponding to
%  the parameter values
%
%  Note:
%     In effort to make RISE modular, this function is available so that one
%     can use a different sampler if needed, but most likely, one should
%     just use available stock samplers.
%
%  ::
%
%     [lnprior, retcode] = log_prior_density(model)
%     [lnprior, retcode] = log_prior_density(model, param)
%     [lnprior, retcode] = log_prior_density(model, param,filtration)
%
%  Args:
%     model (dsge | rise object): model object
%
%     param (column vector): parameter values
%
%     filtration (empty|struct): results from model filtration that can
%     potentially be used for endogenizing priors.
%
%
%  Returns:
%     : [lnprior, retcode]
%
%     - **lnprior** (double): log of prior density function
%     - **retcode**: return code
%
%  See also:
%     - :func:`log_posterior_kernel <dsge.log_posterior_kernel>`
%
%

author: jmaih


text:
function func_YUV_split
% Objective: extract YUV components from the CIF 4:2:0 video file, that is,
% split CIF file into three files, i.e., .Y, .U, .V.
%
% Jing Tian
% Contact: scuteejtian@hotmail.com
% This program is written in 2005 during my postgraduate studying in
% NTU, Singapore.

in_file_name = 'test.cif'
nFrame = 10;

[fid1 message]= fopen(in_file_name,'rb');

if length(strfind(in_file_name, '.qcif')) == 0
    nRow = 288;
    nColumn = 352;
    in_file_name2 = strrep(in_file_name, '.cif', '.Y');
    [fid2 message]= fopen(in_file_name2,'w');
    in_file_name3 = strrep(in_file_name, '.cif', '.U');
    [fid3 message]= fopen(in_file_name3,'w');
    in_file_name4 = strrep(in_file_name, '.cif', '.V');
    [fid4 message]= fopen(in_file_name4,'w');
else
    nRow = 288 / 2;
    nColumn = 352 / 2;
    in_file_name2 = strrep(in_file_name, '.qcif', '.Y');
    [fid2 message]= fopen(in_file_name2,'w');
    in_file_name3 = strrep(in_file_name, '.qcif', '.U');
    [fid3 message]= fopen(in_file_name3,'w');
    in_file_name4 = strrep(in_file_name, '.qcif', '.V');
    [fid4 message]= fopen(in_file_name4,'w');
end

for i = 1: nFrame
    %reading Y component
	img_y = fread(fid1, nRow * nColumn, 'uchar');
    img_y = reshape(img_y, nColumn, nRow);
    img_y = img_y';
    %reading U component
    img_u = fread(fid1, nRow * nColumn / 4, 'uchar');
    img_u = reshape(img_u, nColumn/2, nRow/2);
    img_u = img_u';
    %reading V component
    img_v = fread(fid1, nRow * nColumn / 4, 'uchar');
    img_v = reshape(img_v, nColumn/2, nRow/2);
    img_v = img_v';
    %writing file
	count = fwrite(fid2, img_y', 'uchar');
	count = fwrite(fid3, img_u', 'uchar');
	count = fwrite(fid4, img_v', 'uchar');
end

fclose(fid1);
fclose(fid2);
fclose(fid3);
fclose(fid4);

disp('ok');

author: Sable

repo: mcbench-benchmarks


text:
function tf=isNonFailure(obj)
% return whether this test was a non-failure
%
% tf=isSuccess(obj)
%
% Input:
%   obj                     MoxUnitPassedTestOutcome object
%
% Output:
%   tf                      set to true, because a skipped test is not a
%                           failure
%
% See also: @MoxUnitTestOutcome

    tf=true;

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%Test weka's Naive Bayes implementation on the iris data.

load fisheriris;    %built in to matlab

%Shuffle the data
rand('twister',0);
perm = randperm(150);
meas = meas(perm,:);
species = species(perm,:);

featureNames = {'sepallength','sepalwidth','petallength','petalwidth','class'};

%Prepare test and training sets.
data = [num2cell(meas),species];
train = data(1:120  ,:);
test  = data(121:end,:);

classindex = 5;

%Convert to weka format
train = matlab2weka('iris-train',featureNames,train,classindex);
test =  matlab2weka('iris-test',featureNames,test);

%Train the classifier
nb = trainWekaClassifier(train,'bayes.NaiveBayes');

%Test the classifier
predicted = wekaClassify(test,nb);

%The actual class labels (i.e. indices thereof)
actual = test.attributeToDoubleArray(classindex-1); %java indexes from 0

errorRate = sum(actual ~= predicted)/30


text:
classdef MOEADSTM < ALGORITHM
% <multi/many> <real/integer>
% MOEA/D with stable matching

%------------------------------- Reference --------------------------------
% K. Li, Q. Zhang, S. Kwong, M. Li, and R. Wang, Stable matching-based
% selection in evolutionary multiobjective optimization, IEEE Transactions
% on Evolutionary Computation, 2014, 18(6): 909-923.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate the weight vectors
            [W,Problem.N] = UniformPoint(Problem.N,Problem.M);
            % Size of neighborhood
            T  = ceil(Problem.N/10);

            %% Detect the neighbours of each solution
            B = pdist2(W,W);
            [~,B] = sort(B,2);
            B = B(:,1:T);

            %% Generate random population
            Population = Problem.Initialization();
            z          = min(Population.objs,[],1);
            % Utility for each subproblem
            Pi = ones(Problem.N,1);
            % Old Tchebycheff function value of each solution on its subproblem
            oldObj = max(abs((Population.objs-repmat(z,Problem.N,1))./W),[],2);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                for subgeneration = 1 : 5
                    % Choose I
                    Bounday = find(sum(W<1e-3,2)==Problem.M-1)';
                    I = [Bounday,TournamentSelection(10,floor(Problem.N/5)-length(Bounday),-Pi)];
                    % Generate an offspring for each solution in I
                    P = zeros(length(I),3);
                    for i = 1 : length(I)
                        % Choose the parents
                        if rand < 0.9
                            P(i,:) = B(I(i),randperm(size(B,2),3));
                        else
                            P(i,:) = randperm(Problem.N,3);
                        end
                    end
                    Offspring = OperatorDE(Problem,Population(P(:,1)),Population(P(:,2)),Population(P(:,3)));
                    z         = min([z;Offspring.objs],[],1);

                    % STM selection
                    Population = STM([Population,Offspring],W,z,max(Population.objs,[],1));
                end
                if ~mod(ceil(Problem.FE/Problem.N),10)
                    % Update Pi for each solution
                    newObj    = max(abs((Population.objs-repmat(z,Problem.N,1))./W),[],2);
                    DELTA     = oldObj - newObj;
                    Temp      = DELTA < 0.001;
                    Pi(~Temp) = 1;
                    Pi(Temp)  = (0.95+0.05*DELTA(Temp)/0.001).*Pi(Temp);
                    oldObj    = newObj;
                end
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%{
nohup matlab -nodesktop -nodisplay < main002_instSeg_v1_ftAbsEucMM_epoch83.m 1>main002_instSeg_v1_ftAbsEucMM_epoch83.log &
%}
clear
% close all
clc;

addpath(genpath('../libs'))
path_to_matconvnet = '/home/skong2/scratch/matconvnet-1.0-beta23_modifiedDagnn';

run(fullfile(path_to_matconvnet, 'matlab', 'vl_setupnn'));
addpath(genpath(fullfile('dependencies', 'matconvnet','examples')));
%% read matconvnet model
load('imdb_toydata_v3_from_mnist.mat');
imdb.path = './toydata_v3';
imdb.path_to_dataset = './toydata_v3';
% set GPU
gpuId = 4;
gpuDevice(gpuId);
flagSaveFig = true; % {true false} whether to store the result

% saveFolder = 'semantic_main003_seg_v4_largeLR_l2norm_toyDigitV3';
% modelName = 'softmax_net-epoch-63.mat';

saveFolder = 'main007_instSeg_v1_absEucMM';
modelName = 'softmax_net-epoch-83.mat';

netbasemodel = load( fullfile('./exp', saveFolder, modelName) );
netbasemodel = netbasemodel.net;

netbasemodel.layers(136).block = rmfield(netbasemodel.layers(136).block, 'ignoreAverage');
netbasemodel.layers(135).block = rmfield(netbasemodel.layers(135).block, 'ignoreAverage');

netbasemodel = dagnn.DagNN.loadobj(netbasemodel);
%% 1st mean-shift grouping loop
keepLayerName = sprintf('obj_instSeg_reg');
netbasemodel.layers(netbasemodel.getLayerIndex(keepLayerName)).block.lastLayerName = 'res7_cosSim';
keepLayerName = sprintf('obj_instSeg_MM');
netbasemodel.layers(netbasemodel.getLayerIndex(keepLayerName)).block.lastLayerName = 'res7_cosSim';

weight_for_losses = {'obj_instSeg_reg', 1, 'obj_instSeg_MM', 1};
sName_l2norm = 'res7_l2norm';
gt_name =  sprintf('gt_ins');
GaussianBandwidth = 0.1;
randSampleRatio = 0.2;
for loopIdx = 1:5
    [netbasemodel, sName, sName_l2norm] = addOneLoop_forMeanShiftGrouping(...
        netbasemodel, sName_l2norm, loopIdx,...
        GaussianBandwidth, randSampleRatio);

    % add regression loss
    obj_name = sprintf('loop%d_instSeg_reg', loopIdx);
    netbasemodel.addLayer(obj_name, ...
        InstanceSegRegLoss_randSample('loss', 'cosinesimilarityabsregloss', 'lastLayerName', sName), ... softmaxlog logistic
        {sName, gt_name}, obj_name);
    weight_for_losses{end+1} = obj_name;
    weight_for_losses{end+1} = 1;

    % add max-margin loss
    obj_name = sprintf('loop%d_instSeg_MM', loopIdx);
    netbasemodel.addLayer(obj_name, ...
        InstanceSegMMLoss_randSample('loss', 'cosinesimilaritymmloss', 'marginAlpha_', 0.07, 'adaptiveMM', false, 'lastLayerName', sName), ...
        {sName, gt_name}, obj_name)
    weight_for_losses{end+1} = obj_name;
    weight_for_losses{end+1} = 1;
end
%% show learning rates for all layers
for ii = 1:numel(netbasemodel.layers)
    curLayerName = netbasemodel.layers(ii).name;
    if strfind(curLayerName, 'bn')
        fprintf('%03d, %s\n', ii, curLayerName);
%         net.layers(ii).block.bnorm_moment_type_trn = 'batch';
%         net.layers(ii).block.bnorm_moment_type_tst = 'global';
        netbasemodel.params(netbasemodel.layers(ii).paramIndexes(3)).learningRate = 0.1;
    end
end
netbasemodel.params(netbasemodel.getParamIndex('res6_conv_f')).learningRate = 1;
netbasemodel.params(netbasemodel.getParamIndex('res6_conv_b')).learningRate = 1;
%{
for i = 125:numel(netbasemodel.params)
    if ~isempty(strfind(netbasemodel.params(i).name, '_f')) && isempty(strfind(netbasemodel.params(i).name, '_bn'))
        tmp = netbasemodel.params(i).value;
        tmp = single(randn(size(tmp), 'single')*sqrt(2/size(tmp,4)));
        netbasemodel.params(i).value = tmp;
        fprintf('%s \n', netbasemodel.params(i).name);
    elseif ~isempty(strfind(netbasemodel.params(i).name, '_b')) && isempty(strfind(netbasemodel.params(i).name, '_bn'))
        tmp = netbasemodel.params(i).value;
        netbasemodel.params(i).value = zeros(1, size(tmp,2), 'single');
        fprintf('%s \n', netbasemodel.params(i).name);
    elseif ~isempty(strfind(netbasemodel.params(i).name, '_bn')) && ~isempty(strfind(netbasemodel.params(i).name, '_w'))
        tmp = netbasemodel.params(i).value;
        tmp = ones(size(tmp), 'single');
        netbasemodel.params(i).value = tmp;
        fprintf('%s \n', netbasemodel.params(i).name);
    elseif ~isempty(strfind(netbasemodel.params(i).name, '_bn')) && ~isempty(strfind(netbasemodel.params(i).name, '_b'))
        tmp = netbasemodel.params(i).value;
        tmp = zeros(size(tmp), 'single');
        netbasemodel.params(i).value = tmp;
        fprintf('%s \n', netbasemodel.params(i).name);
    elseif ~isempty(strfind(netbasemodel.params(i).name, '_bn')) && ~isempty(strfind(netbasemodel.params(i).name, '_m'))
        tmp = netbasemodel.params(i).value;
        tmp = zeros(size(tmp), 'single');
        netbasemodel.params(i).value = tmp;
        fprintf('%s \n', netbasemodel.params(i).name);
    else
        fprintf('%s ERROR!!!\n', netbasemodel.params(i).name);
    end
end
%}
for i = 1:numel(netbasemodel.params)
    fprintf('%d\t%25s, \t%.2f',i, netbasemodel.params(i).name, netbasemodel.params(i).learningRate);
    fprintf('\tsize: %dx%dx%dx%d\n', size(netbasemodel.params(i).value,1), size(netbasemodel.params(i).value,2), size(netbasemodel.params(i).value,3), size(netbasemodel.params(i).value,4));
end
%% configure training environment
batchSize = 1;
totalEpoch = 100;
learningRate = 1:totalEpoch;
learningRate = (1.0e-4) * (1-learningRate/totalEpoch).^0.9;

weightDecay=0.0005; % weightDecay: usually use the default value

opts.batchSize = batchSize;
opts.learningRate = learningRate;
opts.weightDecay = weightDecay;
opts.momentum = 0.9 ;

opts.expDir = fullfile('./exp', 'main002_instSeg_v1_ftAbsEucMM_epoch83');
if ~isdir(opts.expDir)
    mkdir(opts.expDir);
end

opts.withSemanticSeg = false ;
opts.withInstanceSeg = true ;
opts.withWeights = false ;

opts.numSubBatches = 1 ;
opts.continue = true ;
opts.gpus = gpuId ;
%gpuDevice(opts.train.gpus); % don't want clear the memory
opts.prefetch = false ;
opts.sync = false ; % for speed
opts.cudnn = true ; % for speed
opts.numEpochs = numel(opts.learningRate) ;
opts.learningRate = learningRate;

for i = 1:2
    curSetName = imdb.sets.name{i};
    idxList = find(imdb.set==i);
%     if i == 1
%         idxList = idxList(1:100);
%     end
    curList = imdb.imgList(idxList);
    opts.(curSetName) = curList;
end

opts.checkpointFn = [];
mopts.classifyType = 'softmax';

rng(777);
bopts = netbasemodel.meta.normalization;
bopts.numThreads = 12;
bopts.imdb = imdb;
%% train
fn = getBatchWrapper4toyDigitV2(bopts);

opts.backPropDepth = inf; % could limit the backprop
prefixStr = [mopts.classifyType, '_'];
opts.backPropAboveLayerName = 'conv1_conv';

trainfn = @cnnTrain;
[netbasemodel, info] = trainfn(netbasemodel, prefixStr, imdb, fn, 'derOutputs', ...
    weight_for_losses, ...
    opts);

%% leaving blank
%{
%}

author: aimerykong


text:
function y= DAB(T,P,MA,MB,sigmA,sigmB,epsA,epsB)
%DAB Calculates gas-phase diffusivity using Wilke-Lee equation, p. 17 Text.
%   DAB(T,P,MA,MB,sigmA,sigmB,epsA,epsB)
%   T = absolute temperature in K
%   P = pressure in bar
%   MA, MB = molecular weights
%   sigmA, sigmB, epsA, epsB = Lennard-Jones parameters
%   in Angstroms and Kelvin, respectively
%   DAB = diffusivity in square cm/sec.
sigmAB=(sigmA+sigmB)/2;
epsAB=sqrt(epsA*epsB);
x=T/epsAB;
a=1.06036;b=0.15610;c=0.19300;d=0.47635;
e1=1.03587;f=1.52996;g=1.76474;h=3.89411;
MAB=2*(1/MA+1/MB)^-1;
omega=a/x^b+c/exp(d*x)+e1/exp(f*x)+g/exp(h*x);
y=0.001*(3.03-0.98/sqrt(MAB))*T^1.5/(P*sigmAB^2*omega*sqrt(MAB));


text:
function chx = bbp_pi ( d )

%*****************************************************************************80
%
%% BBP_PI implements the Bailey-Borwein-Plouffe algorithm.
%
%  Discussion:
%
%    The BBP algorithm can be used to compute a a few hex digits of pi starting
%    at any position.
%
%    In particular, bbp_pi ( d ) is a char string of hex digits d+1 through
%    d+13 in the hexadecimal expansion of pi.
%
%    This function does not require extended precision arithmetic or
%    symbolic computation.
%
%    The results are usually accurate to about 11 or 12 of the 13 digits.
%
%    This program is derived from a C program by David H. Bailey dated
%    2006-09-08, http://www.experimentalmath.info/bbp-codes/piqpr8.c
%
%    For many other references: Google "BBP Pi".
%
%  Licensing:
%
%    Copyright (c) 2011, The MathWorks, Inc.
%    All rights reserved.
%
%    Redistribution and use in source and binary forms, with or without
%    modification, are permitted provided that the following conditions are
%    met:
%
%        * Redistributions of source code must retain the above copyright
%          notice, this list of conditions and the following disclaimer.
%        * Redistributions in binary form must reproduce the above copyright
%          notice, this list of conditions and the following disclaimer in
%          the documentation and/or other materials provided with the distribution
%        * Neither the name of the The MathWorks, Inc. nor the names
%          of its contributors may be used to endorse or promote products derived
%          from this software without specific prior written permission.
%
%    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%    AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%    IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%    ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%    LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%    POSSIBILITY OF SUCH DAMAGE.
%
%  Modified:
%
%    27 February 2012
%
%  Author:
%
%    Cleve Moler
%
%  Reference:
%
%    Cleve Moler,
%    Cleve's Corner, "Computing Pi",
%    http://www.mathworks.com/company/newsletters/news_notes/2011/
%
%  Parameters:
%
%    Input, integer D, the number of decimal digits desired.
%
%    Output, symbolic P, the value of pi to D digits.
%
   s1 = series ( 1, d );
   s2 = series ( 4, d );
   s3 = series ( 5, d );
   s4 = series ( 6, d );
   P = 4. * s1 - 2. * s2 - s3 - s4;
   P = P - floor ( P ) + 1.;
   chx = hexchar ( P, 13 );

  return
end
function s = series ( m, d )

%*****************************************************************************80
%
%% SERIES
%
% s = series(m,d) = sum_k 16^(d-k)/(8*k+m)
% using the modular powering technique.
%
  s = 0;
  k = 0;
  t = Inf;

  while k < 13 || t > eps
    ak = 8 * k + m;
    if k < d
      t = powmod(16, d - k, ak) / ak;
     else
       t = 16^(d - k) / ak;
     end
     s = mod(s + t, 1);
     k = k + 1;
   end

  return
end
function r = powmod ( b, p, a )

%*****************************************************************************80
%
%% POWMOD computes mod(b^p,a) without computing b^p.
%
  persistent twop

  if isempty ( twop )
    twop = 2.^(0:25)';
  end

  if a == 1
    r = 0;
    return
  end

  n = find ( p <= twop, 1, 'first' );

  pt = twop ( n );
  r = 1;
  for j = 1 : n
    if p >= pt
      r = mod ( b * r, a );
      p = p - pt;
    end
    pt = 0.5 * pt;
    if pt >= 1
      r = mod ( r * r, a );
    end
  end

  return
end
function chx = hexchar ( s, n )

%*****************************************************************************80
%
%% HEXCHAR returns hex digits.
%
% chx(s,n) = string of the first n hex digits of s.
%
  hx = '0123456789ABCDEF';

  s = abs ( s );

  for j = 1 : n
    s = 16. * mod ( s, 1 );
    chx(j) = hx ( floor ( s ) + 1 );
  end

  return
end

author: johannesgerer


text:
function c = crossProduct3d(a,b)
%CROSSPRODUCT3D Vector cross product faster than inbuilt MATLAB cross.
%
%   C = crossProduct3d(A, B)
%   returns the cross product of the two 3D vectors A and B, that is:
%       C = A x B
%   A and B must be N-by-3 element vectors. If either A or B is a 1-by-3
%   row vector, the result C will have the size of the other input and will
%   be the  concatenation of each row's cross product.
%
%   Example
%     v1 = [2 0 0];
%     v2 = [0 3 0];
%     crossProduct3d(v1, v2)
%     ans =
%         0   0   6
%
%
%   Class support for inputs A,B:
%      float: double, single
%
%   See also DOT.

% ------
% Author: Sven Holcombe
% e-mail: N/A

% HISTORY
% 2017-11-24 rename from vectorCross3d to crossProduct3d

% size of inputs
sizeA = size(a);
sizeB = size(b);

% Initialise c to the size of a or b, whichever has more dimensions. If
% they have the same dimensions, initialise to the larger of the two
switch sign(numel(sizeA) - numel(sizeB))
    case 1
        c = zeros(sizeA);
    case -1
        c = zeros(sizeB);
    otherwise
        c = zeros(max(sizeA, sizeB));
end

c(:) = bsxfun(@times, a(:,[2 3 1],:), b(:,[3 1 2],:)) - ...
       bsxfun(@times, b(:,[2 3 1],:), a(:,[3 1 2],:));

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function obj = updateImplicitFunctionSurface(obj, surfaceIndex)

%-AXIS INDEX-%
axIndex = obj.getAxisIndex(obj.State.Plot(surfaceIndex).AssociatedAxis);

%-CHECK FOR MULTIPLE AXES-%
[xsource, ysource] = findSourceAxis(obj,axIndex);

%-SURFACE DATA STRUCTURE- %
image_data = get(obj.State.Plot(surfaceIndex).Handle);
figure_data = get(obj.State.Figure.Handle);

%-AXIS DATA-%
eval(['xaxis = obj.layout.xaxis' num2str(xsource) ';']);
eval(['yaxis = obj.layout.yaxis' num2str(ysource) ';']);

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

%-surface xaxis-%
obj.data{surfaceIndex}.xaxis = ['x' num2str(xsource)];

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

%-surface yaxis-%
obj.data{surfaceIndex}.yaxis = ['y' num2str(ysource)];

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

%-surface type-%
obj.data{surfaceIndex}.type = 'surface';

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

%-getting x,y,z surface data-%

strf = func2str(image_data.Function);
ind1 = strfind(strf, '('); ind1 = ind1(1)+1;
ind2 = strfind(strf, ')'); ind2 = ind2(1)-1;
vars = split(strf(ind1:ind2), ',');

strf = [strf(ind2+2:end) '==0'];
strf = replace(strf, vars{1}, 'Xx');
strf = replace(strf, vars{2}, 'Yy');
strf = replace(strf, vars{3}, 'Zz');

syms Xx Yy Zz;
f = eval(strf);
s = solve(f, Zz);

x = image_data.XRange;
y = image_data.YRange;
z = image_data.ZRange;
N = 400;

[Xx,Yy] = meshgrid(linspace(x(1),x(2),N), linspace(y(1),y(2),N));
X = []; Y = []; Z = [];

for n = 1:length(s)
    X = [X; Xx];
    Y = [Y; Yy];
    Z = [Z; eval(s(n))];
end

clear Xx Yy Zz;
Z(Z < z(1)) = nan; Z(Z > z(2)) = nan;
X(Z < z(1)) = nan; X(Z > z(2)) = nan;
Y(Z < z(1)) = nan; Y(Z > z(2)) = nan;

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

%-surface x-%
obj.data{surfaceIndex}.x = X;

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

%-surface y-%
obj.data{surfaceIndex}.y = Y;

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

%-surface z-%
obj.data{surfaceIndex}.z = Z;

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

%- setting grid mesh by default -%
% x-direction
mden = image_data.MeshDensity;
xsize = (x(2) - x(1)) / mden;
obj.data{surfaceIndex}.contours.x.start = x(1);
obj.data{surfaceIndex}.contours.x.end = x(2);
obj.data{surfaceIndex}.contours.x.size = xsize;
obj.data{surfaceIndex}.contours.x.show = true;
obj.data{surfaceIndex}.contours.x.color = 'black';
% y-direction
ysize = (y(2) - y(1)) / mden;
obj.data{surfaceIndex}.contours.y.start = y(1);
obj.data{surfaceIndex}.contours.y.end = y(2);
obj.data{surfaceIndex}.contours.y.size = ysize;
obj.data{surfaceIndex}.contours.y.show = true;
obj.data{surfaceIndex}.contours.y.color = 'black';
% z-direction
zsize = (z(2) - z(1)) / mden;
obj.data{surfaceIndex}.contours.z.start = z(1);
obj.data{surfaceIndex}.contours.z.end = z(2);
obj.data{surfaceIndex}.contours.z.size = zsize;
obj.data{surfaceIndex}.contours.z.show = true;
obj.data{surfaceIndex}.contours.z.color = 'black';

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

%-image colorscale-%

cmap = figure_data.Colormap;
len = length(cmap)-1;

for c = 1: length(cmap)
    col = 255 * cmap(c, :);
    obj.data{surfaceIndex}.colorscale{c} = { (c-1)/len , ['rgb(' num2str(col(1)) ',' num2str(col(2)) ',' num2str(col(3)) ')'  ]  };
end

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

%-image surfacescale-%
obj.data{surfaceIndex}.surfacecolor = Z;

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

%-surface name-%
obj.data{surfaceIndex}.name = image_data.DisplayName;

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

%-surface showscale-%
obj.data{surfaceIndex}.showscale = false;

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

%-surface visible-%
obj.data{surfaceIndex}.visible = strcmp(image_data.Visible,'on');

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

leg = get(image_data.Annotation);
legInfo = get(leg.LegendInformation);

switch legInfo.IconDisplayStyle
    case 'on'
        showleg = true;
    case 'off'
        showleg = false;
end

obj.data{surfaceIndex}.showlegend = showleg;

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

end

author: plotly


text:
function varargout = pulldown(varargin)

%PULLDOWN Standalone pulldown menus with title.
%   [A, B, C, ...] = PULLDOWN(TITLE, X, Y, Z, ...) will invoke simple
%   pulldown menus with title TITLE, and lists of choices in cell arrays,
%   X, Y, Z, ..., and provides the choice index for list X in A, the choice
%   index for list Y in B, etc.
%
%   CHOICES = PULLDOWN(TITLE, X, Y, Z, ...) does the same thing, except
%   CHOICES is a numerical array containing all the individual choice
%   indices from X, Y, Z, etc.
%
%   [...] = PULLDOWN(..., V) does the same thing, but will initialize the
%   pulldown menus to the values in vector V.  V may have fewer or more
%   choices than the number of pulldown menus.
%
%   Examples
%   --------
%       [letter, number] = pulldown('Test', {'A', 'B', 'C'}, num2cell(1:5))
%
%       choices = pulldown('Mission', {'117', '116', '115', '121', '114'})
%
%       % initialize to choices 'B' and 4
%       choices = pulldown('Test', {'A', 'B', 'C'}, num2cell(1:5), [2 4])
%
%   See also UICONTROL, UIWAIT, GET, SET, VARARGIN.
%
%   Version 1.2 Kevin Crosby

% DATE      VER  NAME          DESCRIPTION
% 01-23-07  1.0  K. Crosby     First Release.  Promoted from sub-function.
% 07-11-07  1.1  K. Crosby     Added initial settings for pulldown.
% 07-23-07  1.2  K. Crosby     Changed to modal.

% Contact: Kevin.L.Crosby@gmail.com

m = zeros(1, nargin-1); % maximums
n = 0; % number of pulldowns
while n < nargin-1 && iscell(varargin{n+2})
  arg = n + 2;
  n = n + 1;
  m(n) = length(varargin{arg});
  if ~iscellstr(varargin{arg})
    varargin{arg} = cellstr(num2str(cell2mat(varargin{arg})'))';
  end % if ~iscellstr(varargin{arg})
end % while n
m(n+1:end) = [];

u = ones(1, n); % defaults
k = nargin - n - 1; % number of default overrides
switch k
  case 0
    % do nothing
  case 1
    l = length(varargin{end});
    u(1:min(n, l)) = varargin{end}(1:min(n, l));
  otherwise
    u(1:min(n, k)) = [varargin{n+2:end}];
end % switch k

% ensure within bounds
u = round(u); % ensure integers
u(u < 1) = 1; % force minimum to be 1
u(u > m) = m(u > m); % force maximum to be m

h = 1 / (n+1);
x = .15;      % figure width
y = .02*(n+1); % figure height

popup = zeros(1, n);

fig = figure('MenuBar', 'none', 'Resize', 'off', 'CloseRequestFcn', '', ...
  'Name', varargin{1}, 'NumberTitle', 'off', 'ToolBar', 'none', ...
  'WindowStyle', 'modal', ...
  'Units', 'normalized', 'Position', [.5-x/2 .5-y/2 x y]);
for p = 1:n
  popup(p) = uicontrol('Parent', fig, 'Style', 'popupmenu', ...
    'Units', 'normalized', 'Position', [0 (n+1-p)*h 1 h], ...
    'String', varargin{p+1}, 'Value', u(p), 'Callback', {@callbackfcn, p});
end
uicontrol('Parent', fig, 'Style', 'pushbutton', ...
  'Units', 'normalized', 'Position', [0 0 1 h], ...
  'String', 'OK', 'Callback', 'closereq');
set(0, 'UserData', u); % default choices
uiwait(fig);
choices = get(0, 'UserData');

% the following was added while promoting to a function
switch nargout
  case {0, 1}
    varargout = num2cell(choices, 2);
  otherwise
    varargout = num2cell(choices);
end % switch nargout

% call back function
function callbackfcn(src, eventdata, p) % src, eventdata, varargin ...

if isempty(eventdata) || ~isstruct(eventdata)
  eventdata.Key = '';
end % isempty(eventdata) || ~isstruct(eventdata)

u = get(0, 'UserData');
u(p) = get(src, 'Value');
set(0, 'UserData', u);