text:
function mesh = mfmAdjustSpacing(mesh, unfoldMesh, params, unfolded2D, insideNodes);
%
%  mesh = mfmAdjustSpacing(mesh, unfoldMesh, params, unfolded2D,
%  insideNodes);
%
% Author: Winawer
% Purpose:
%   Adjust the spacing of the points in a flattened mesh so they don't
%   bunch up too much.  The method is to find a cartesian grid within the
%   data, find a Delaunay triangulation of this grid, and then use a series
%   of affine transformations to transform each of the triangles to an
%   equal area representation with the proper grid topology.  This is
%   explained in more detail in flatAdjustSpacing
%
%   Sub-routine derived from Alex's unfoldMeshFromGUI code.
%
% See Also:  unfoldMeshFromGUI

% Get all the variables we may need
meshFileName     = params.meshFileName;
grayFileName     = params.grayFileName;
flatFileName     = params.flatFileName;
startCoords      = params.startCoords;
scaleFactor      = params.scaleFactor;
perimDist        = params.perimDist;
statusHandle     = params.statusHandle;
busyHandle       = params.busyHandle;
spacingMethod    = params.spacingMethod;
adjustSpacing    = params.adjustSpacing;
gridSpacing      = params.gridSpacing;
showFigures      = params.showFigures;
saveExtra        = params.saveExtra;
truePerimDist    = params.truePerimDist;
hemi             = params.hemi;
nperims          = params.NPERIMS;
saveIntermeidate = params.SAVE_INTERMEDIATE;
numberOfSteps    = params.NUMBEROFSTEPS;

if (adjustSpacing)
    str = sprintf('Spacing points (method = %s)',spacingMethod);
    statusStringAdd(statusHandle,str);

    unfoldMesh.maxFractionDist = 0.8;
    unfoldMesh.gridSpacing = gridSpacing;
    unfoldMesh.locs2d = unfolded2D(:,1:2);
    unfoldMesh.startCoords = startCoords;
    unfoldMesh.scaleFactor = scaleFactor;

    % The user can select Cartesian or Polar spacing methods.  Only
    % Cartesian is working now, though, I think -- BW
    [newLocs2d,goodIdx] = flatAdjustSpacing(unfoldMesh,spacingMethod);
end

if (showFigures)
    statusStringAdd(statusHandle,['Displaying unfolded mesh']);
    h = unfoldMeshFigure(unfoldMesh);

	[p f] = fileparts(flatFileName);
    ensureDirExists(p);
	savePath = fullfile(p, [f,'_Unfolded Mesh.png']);
	saveas(gcf, savePath);
	fprintf('[%s]: Saved %s.\n', mfilename, savePath);
end

% Finally - the mapping of grey to mesh points takes place using the entire mesh.
% Therefore, we need to generate X for the mesh as well as the unfold mesh

% insideNodes is an array of indices into the original (non-bounded) mesh.
% Each entry in insideNodes relates a point in the unfoldMesh to a point in
% the original mesh

% Recover the perimeter and internal points
mesh.X=zeros(length(mesh.uniqueVertices),2);

% In order to deal with a cropped mesh (as generated by flatAdjustSpacing)
% we need to compute the ... Alex?
if (adjustSpacing)

    mesh.X(insideNodes(goodIdx),:)=newLocs2d;
    hasCoords=[insideNodes(goodIdx)];

else

    unfoldToOrigPerimeter=insideNodes(unfoldMesh.orderedUniquePerimeterPoints);
    unfoldToOrigInside=insideNodes(unfoldMesh.internalNodes);

    mesh.X(unfoldToOrigPerimeter,:)=unfoldMesh.X_zero;
    mesh.X(unfoldToOrigInside,:)=unfoldMesh.X;
    hasCoords=[unfoldToOrigPerimeter(:);unfoldToOrigInside(:)];

end

coords=mesh.X(hasCoords,:);
dists=mesh.dist(hasCoords);

% use griddata to image the distance map
warning off MATLAB:griddata:DuplicateDataPoints;
mesh.distMap=makeMeshImage(coords,dists,128);
warning on MATLAB:griddata:DuplicateDataPoints;

ZI=mesh.distMap;

if (showFigures),
	h = unfoldDistFigure(mesh);

	[p f ext] = fileparts(flatFileName);
	savePath = fullfile(p, [f ' Unfold Distance.png']);
	saveas(h, savePath);
	fprintf('[%s]: Saved %s.\n', mfilename, savePath);
end

% Record which nodes in the big mesh are in the unfold
mesh.insideNodes=insideNodes;

if (saveIntermeidate)
    statusStringAdd(statusHandle,'Saving intermediate data.');
    save ('meshOutTemp.mat','mesh');
end

return

%----------------------------------
% ---- SUBROUTINES
%----------------------------------
%----------------------------------
function h = unfoldMeshFigure(unfoldMesh);
h = figure;
hold off;
gplot(unfoldMesh.N,unfoldMesh.X);
title ('Unfolded mesh'); axis equal; axis off; zoom on
return;

%----------------------------------
function h = unfoldDistFigure(mesh)
h = figure;
imagesc(mesh.distMap);
axis image; colormap hot; title('Manifold distance map'); colorbar;

% Added the abs to cope with the fact that we sometimes get -ve numbers here...
% imagesc(abs(log(areaErrorMap)));
% axis image;
% colorbar;
% title('Log 3D/2D area');
% colormap hot;

return;

author: vistalab

repo: vistasoft


text:
function [Y,W,SetupStruc] = Process_ICA_Sawada(s,Transfer,SetupStruc)
K = SetupStruc.ICA_Sawada.K;
hop = SetupStruc.ICA_Sawada.hop;
win = hanning(K,'periodic');
win = win/sqrt(sum(win(1:hop:K).^2));
SetupStruc.ICA_Sawada.win = win;  % Preserve 'win' in 'SetupStruc'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = size(s,2);
for i = 1:N
    X(:,:,i) = fft(enframe(s(:,i),win,hop)');
end
frame_N = size(X,2);
K_m = K/2+1;
Num = size(Transfer,3);
Y = zeros((frame_N-1)*hop+K,Num);
Y_f = zeros(size(X,1),size(X,2),Num);
Y_P = zeros(frame_N,Num,K_m);
%%%%%%%%%%%%%%%%%%%%%%%%%% Obtain processing matrix 'W'
theta = 10^-4;
W = zeros(Num,N,K_m);
A = zeros(1001,K/2)-1; %%%% Show the decrease of non-linear correlation, ICA max iterations 1000
for i = 2:K_m
    X_f = permute(X(i,:,:),[3 2 1]);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% PCA and ICA processing
    [E,D] = PCA(X_f,1,Num);
    V = sqrt(D)\E';
    X_f = V*X_f;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Steer = permute(Transfer(i,:,:),[2 3 1]);
%     Ori = V*Steer;
%     if rcond(Ori)<theta
%         Ori = Ori+eye(Num)*min(diag(Ori));
%     end
%     [Y_,W_ICA,A] = FDICA(X_f,inv(Ori),A,i);  %%% 'A', 'i' record the decrease for observation
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [Y_,W_ICA,A] = FDICA(X_f,eye(Num),A,i);  %%% 'A', 'i' record the decrease for observation
    W(:,:,i) = W_ICA*V;
    Y_P(:,:,i) = Y_.';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Process the ambiguity of permutation and amplitude
P = Permu_Sawada(W,Y_P,SetupStruc,'all');  %%%% Options: 'DOA','cor', 'all'
for i = 2:K_m
    W(:,:,i) = P(:,:,i)*W(:,:,i);
    Y_ = permute(Y_P(:,:,i),[2 1 3]);
    Y_ = P(:,:,i)*Y_;
    Y_f(i,:,:) = Y_.';
    if(i~=K_m)
        Y_f(K+2-i,:,:) = Y_';
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Recover signals
if(K/hop==2)
    win = ones(K,1);
end
for i = 1:Num
    Y(:,i) = overlapadd(real(ifft(Y_f(:,:,i)))',win,hop);
end
return;

repo: Voice-Separation-and-Enhancement


text:
function [ obblist, labellist ] = visualizeReloffset( leafreps, labeldefs, dim )
%VISUALIZERELPOS
% this function only considers the artificial dataset: exact one bed, no
% more than two objects in each category
classnum = 7;
if(dim==2)
    % turn the relpos into 3D
    nleafreps = zeros(classnum*(dim+1)*2+classnum,size(leafreps,2));
    for i = 1:classnum
        nleafreps((i-1)*(dim+1)*2+1,:) = leafreps((i-1)*dim*2+1,:);
        nleafreps((i-1)*(dim+1)*2+3,:) = leafreps((i-1)*dim*2+2,:);
        nleafreps((i-1)*(dim+1)*2+4,:) = leafreps((i-1)*dim*2+3,:);
        nleafreps((i-1)*(dim+1)*2+6,:) = leafreps((i-1)*dim*2+4,:);
    end
    nleafreps(end-classnum+1:end,:) = leafreps(end-classnum+1:end,:);
    leafreps = nleafreps;
end

offsetnum = 6;
objnum = size(leafreps,2);
relvecs = leafreps(1:classnum*offsetnum,:);
labelvecs = leafreps(classnum*offsetnum+1:end,:);

% recon the size
objsizelist = zeros(offsetnum/2,objnum);
for i = 1:objnum
    rv = relvecs(:,i);

    sizelist = [];
    for j = 1:classnum
        rvj = rv((j-1)*offsetnum+1:j*offsetnum);
        if(~all(rvj==0))
            sizelist = [sizelist,rvj(1:offsetnum/2)-rvj(offsetnum/2+1:end)];
        end
    end
    si = mean(sizelist');

    objsizelist(:,i)=si(:);
end

% assume there's only one bed and located in [0,0]
bed_label = labeldefs(:,1);
dist = zeros(objnum,1);
for i = 1:objnum
    label = labelvecs(:,i);
    dis = norm(label-bed_label,2);
    dist(i) = dis;
end
[~,I] = min(dist);
bed_index = I;

obblist = zeros(12,objnum);
info = relvecs(1:offsetnum,:);% we only use the relpos of each object relative to the bed, based on the assumption
% actually, the redundancy will cause mess.
bedminp = [0;0;0]-objsizelist(:,bed_index);
anchor = [bedminp/2;bedminp/2];
for i = 1:objnum
    bbox = info(:,i)+anchor;
    maxp = bbox(1:offsetnum/2);
    minp = bbox(offsetnum/2+1:end);

%     maxp = [maxp(1);0.1;maxp(2)];
%     minp = [minp(1);0;minp(2)];
    cent = (maxp+minp)/2;
    orient = [1;0;0;0;1;0];
    s = maxp-minp;
    obb = [cent;orient;s];
    obblist(:,i) = obb;
end

% recon the labellist
tlabellist = {'bed','stand','lamp','rug','ottoman','person','floor'};
labellist = {};
for i = 1:objnum
    label = labelvecs(:,i);
    dist = zeros(classnum,1);
    for j = 1:classnum
        clabel = labeldefs(:,j);
        dis = norm(label-clabel,2);
        dist(j) = dis;
    end
    [~,I] = min(dist);
    labellist{length(labellist)+1} = tlabellist{I};
end

end

repo: GRAINS


text:
%kkgentb 'kgentb: Creates a toolbox'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kgentb.pane file
%
% Parameters:
% String: tb 'Toolbox Name', default: ' ': 'toolbox to create'
% String: path 'Toolbox Path', default: ' ': 'toolbox path'
% String: title 'Toolbox Title ', default: ' ': 'toolbox title'
% String: info 'Toolbox Info ', default: ' ': 'toolbox description file'
% String: author 'Toolbox Author', default: ' ': 'toolbox author'
% String: email 'Author Email ', default: ' ': 'toolbox author email'
% Toggle: add 'Add Toolbox Instead of Generate', default: 0: 'add toolbox instead of generate?'
%
% Example: kkgentb( {'tb',' ';'path',' ';'title',' ';'info',' ';'author',' ';'email',' ';'add',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kgentb - Creates a toolbox from scratch.
%
%  DESCRIPTION
% This program takes the specifications provided, and creates a new toolbox
% object.  If a toolbox by the name given (or located at the path given)
% already exists an error message will be given and the program will exit.
%
% If the [-add] flag is specified, the appropriate toolbox entry will be
% added to the toolbox table, without generating the toolbox object (in this
% case, it is understood that the toolbox in question already exists).
% This allows users to add newly-acquired toolboxes to their list of available
% toolboxes automatically.  Note that the entry will not be added if the
% specified path is not the top level path to a real (existing) toolbox.
%
%
%
%  EXAMPLES
% None
%
%  "SEE ALSO"
% None
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkgentb(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kkgentb(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'tb', ' ';'path', ' ';'title', ' ';'info', ' ';'author', ' ';'email', ' ';'add', 0};
maxval={0,0,0,0,0,0,0};
minval={0,0,0,0,0,0,0};
istoggle=[0,0,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'String','String','String','String','String','String','Toggle'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'kgentb"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [header,tracks] = trk_read(filePath)
%TRK_READ - Load TrackVis .trk files
%TrackVis displays and saves .trk files in LPS orientation. After import, this
%function attempts to reorient the fibers to match the orientation of the
%original volume data.
%
% Syntax: [header,tracks] = trk_read(filePath)
%
% Inputs:
%    filePath - Full path to .trk file [char]
%
% Outputs:
%    header - Header information from .trk file [struc]
%    tracks - Track data structure array [1 x nTracks]
%      nPoints - # of points in each streamline
%      matrix  - XYZ coordinates (in mm) and associated scalars [nPoints x 3+nScalars]
%      props   - Properties of the whole tract (ex: length)
%
% Example:
%    exDir           = '/path/to/along-tract-stats/example';
%    subDir          = fullfile(exDir, 'subject1');
%    trkPath         = fullfile(subDir, 'CST_L.trk');
%    [header tracks] = trk_read(trkPath);
%
% Other m-files required: none
% Subfunctions: get_header
% MAT-files required: none
%
% See also: http://www.trackvis.org/docs/?subsect=fileformat
%           http://github.com/johncolby/along-tract-stats/wiki/orientation

% Author: John Colby (johncolby@ucla.edu)
% UCLA Developmental Cognitive Neuroimaging Group (Sowell Lab)
% Mar 2010

% Parse in header
fid    = fopen(filePath, 'r');
header = get_header(fid);

% Check for byte order
if header.hdr_size~=1000
    fclose(fid);
    fid    = fopen(filePath, 'r', 'b'); % Big endian for old PPCs
    header = get_header(fid);
end

if header.hdr_size~=1000, error('Header length is wrong'), end

% Check orientation
[tmp ix] = max(abs(header.image_orientation_patient(1:3)));
[tmp iy] = max(abs(header.image_orientation_patient(4:6)));
iz = 1:3;
iz([ix iy]) = [];

% Fix volume dimensions to match the reported orientation.
header.dim        = header.dim([ix iy iz]);
header.voxel_size = header.voxel_size([ix iy iz]);

% Parse in body
if header.n_count > 0
	max_n_trks = header.n_count;
else
	% Unknown number of tracks; we'll just have to read until we run out.
	max_n_trks = inf;
end

% It's impossible to preallocate the "tracks" variable because we don't
% know the number of points on each curve ahead of time; we find out by
% reading the file.  The line below suppresses preallocation warnings.
%#ok<*AGROW>

iTrk = 1;
while iTrk <= max_n_trks
	pts = fread(fid, 1, 'int');
	if feof(fid)
		break;
	end
    tracks(iTrk).nPoints = pts;
    tracks(iTrk).matrix  = fread(fid, [3+header.n_scalars, tracks(iTrk).nPoints], '*float')';
    if header.n_properties
        tracks(iTrk).props = fread(fid, header.n_properties, '*float');
    end

    % Modify orientation of tracks (always LPS) to match orientation of volume
    coords = tracks(iTrk).matrix(:,1:3);
    coords = coords(:,[ix iy iz]);
    if header.image_orientation_patient(ix) < 0
        coords(:,ix) = header.dim(ix)*header.voxel_size(ix) - coords(:,ix);
    end
    if header.image_orientation_patient(3+iy) < 0
        coords(:,iy) = header.dim(iy)*header.voxel_size(iy) - coords(:,iy);
    end
    tracks(iTrk).matrix(:,1:3) = coords;
	iTrk = iTrk + 1;
end

if header.n_count == 0
	header.n_count = length(tracks);
end

fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function header = get_header(fid)

header.id_string                 = fread(fid, 6, '*char')';
header.dim                       = fread(fid, 3, 'short')';
header.voxel_size                = fread(fid, 3, 'float')';
header.origin                    = fread(fid, 3, 'float')';
header.n_scalars                 = fread(fid, 1, 'short')';
header.scalar_name               = fread(fid, [20,10], '*char')';
header.n_properties              = fread(fid, 1, 'short')';
header.property_name             = fread(fid, [20,10], '*char')';
header.vox_to_ras                = fread(fid, [4,4], 'float')';
header.reserved                  = fread(fid, 444, '*char');
header.voxel_order               = fread(fid, 4, '*char')';
header.pad2                      = fread(fid, 4, '*char')';
header.image_orientation_patient = fread(fid, 6, 'float')';
header.pad1                      = fread(fid, 2, '*char')';
header.invert_x                  = fread(fid, 1, 'uchar');
header.invert_y                  = fread(fid, 1, 'uchar');
header.invert_z                  = fread(fid, 1, 'uchar');
header.swap_xy                   = fread(fid, 1, 'uchar');
header.swap_yz                   = fread(fid, 1, 'uchar');
header.swap_zx                   = fread(fid, 1, 'uchar');
header.n_count                   = fread(fid, 1, 'int')';
header.version                   = fread(fid, 1, 'int')';
header.hdr_size                  = fread(fid, 1, 'int')';

author: yetianmed

repo: subcortex

save_path: github-repos/MATLAB/yetianmed-subcortex


text:
function adj = tet_mesh_order10_adj_set ( node_num, tetra_num, tetra_node, ...
  adj_num, adj_row )

%*****************************************************************************80
%
%% TET_MESH_ORDER10_ADJ_SET sets the nodal adjacency matrix.
%
%  Discussion:
%
%    A compressed format is used for the nodal adjacency matrix.
%
%    It is assumed that we know ADJ_NUM, the number of adjacency entries
%    and the ADJ_ROW array, which keeps track of the list of slots
%    in ADJ where we can store adjacency information for each row.
%
%    We essentially repeat the work of TET_MESH_ORDER4_ADJ_COUNT, but
%    now we have a place to store the adjacency information.
%
%    A copy of the ADJ_ROW array is useful, as we can use it to keep track
%    of the next available entry in ADJ for adjacencies associated with
%    a given row.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer TETRA_NUM, the number of tetrahedrons.
%
%    Input, integer TETRA_NODE(10,TETRA_NUM), the nodes that make up the
%    tetrahedrons.
%
%    Input, integer ADJ_NUM, the total number of adjacency relationships,
%
%    Input, integer ADJ_ROW(NODE_NUM+1), the ADJ pointer array.
%
%    Output, integer ADJ(ADJ_NUM), the adjacency information.
%

%
%  Each order 10 tetrahedron defines 45 adjacency pairs.
%
  k = 0;
  for i = 1 : 9
    for j = i + 1 : 10
      pair(1,k*tet_num+1:(k+1)*tet_num) = tet_node(i,1:tet_num);
      pair(2,k*tet_num+1:(k+1)*tet_num) = tet_node(j,1:tet_num);
      k = k + 1;
    end
  end
%
%  Force the nodes of each pair to be listed in ascending order.
%
  pair_num = 45 * tetra_num;
  pair = i4col_sort2_a ( 2, pair_num, pair );
%
%  Rearrange the columns in ascending order.
%
  pair = i4col_sort_a ( 2, pair_num, pair );
%
%  Mark all entries of ADJ so we will know later if we missed one.
%
  adj(1:adj_num) = -1;
%
%  Copy the ADJ_ROW array and use it to keep track of the next
%  free entry for each row.
%
  adj_row_copy(1:node_num) = adj_row(1:node_num);
%
%  Now set up the ADJ_ROW counts.
%
  for k = 1 : pair_num

    if ( 1 < k )
      if ( pair(1,k-1) == pair(1,k) && pair(2,k-1) == pair(2,k) )
        continue
      end
    end

    i = pair(1,k);
    j = pair(2,k);

    adj(adj_row_copy(i)) = j;
    adj_row_copy(i) = adj_row_copy(i) + 1;
    adj(adj_row_copy(j)) = i;
    adj_row_copy(j) = adj_row_copy(j) + 1;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [hAx]=axesArea(varargin)
%function [hAx]=axesArea(hAx,varargin)
%
% Set the margine of the axis by specifying a vector of distance to the figure edges.
%
% Input:
%  [varargin=hAx] = axis handle
%    [varargin=p] = position spec: 1, 2 or 4 element vector that specify the distance
%                   from the edges of the figure by a percentage number between (0-49).
%                   If 1 element it is used for all margins.
%                   If 2 elements, p=[x-margins, y-margins].
%                   If 4 elements, p=[left, lower, right, upper] margins.
% Output:
%  [hAx] = axis handle of the axes.
%
% See also: axis, axes, ishandle, set
%
% By: Michael Wu  --  michael.wu@lithium.com (Mar 2009)
%
%====================

% Check if 1st input is axis handle
%--------------------
if ishghandle(varargin{1},'axes')
	hAx=varargin{1};
	p=varargin{2};
else
	hAx=gca;
	p=varargin{1};
end

% Process input arguments
%--------------------
p(p<0)=0;
p(p>49)=49;
p=p/100;

% Compute position property to be set
%--------------------
switch length(p)
	case 1
		xmin=p;
		ymin=xmin;
		xlen=1-2*p;
		ylen=xlen;
	case 2
		xmin=p(1);
		ymin=p(2);
		xlen=1-2*p(1);
		ylen=1-2*p(2);
	case 4
		xmin=p(1);
		ymin=p(2);
		xlen=1-p(1)-p(3);
		ylen=1-p(2)-p(4);
	otherwise
		% Default Matlab position setting
		%--------------------
		xmin=0.13;
		ymin=0.11;
		xlen=0.775;
		ylen=0.815;
end

% Set new position property of the axes
%--------------------
set(hAx,'position',[xmin ymin xlen ylen]);

author: Sable

repo: mcbench-benchmarks


text:
function [area] = triarea(pp,tt)
%TRIAREA calc. triangle areas for a 2-simplex triangulation
%embedded in the two-dimensional plane.
%   [AREA] = TRIAREA(VERT,TRIA) returns the signed triangle
%   areas, where AREA is a T-by-1 vector, VERT is a V-by-2
%   array of XY coordinates, and TRIA is a T-by-3 array of
%   vertex indexing, where each row defines a triangle, such
%   that VERT(TRIA(II,1),:), VERT(TRIA(II,2),:) and VERT(
%   TRIA(II,3),:) are the coordinates of the II-TH triangle.
%
%   See also TRISCR2, TRIANG2, TRIBAL2

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

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

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

    nnod = size(pp,1) ;

%---------------------------------------------- basic checks
    if (min(min(tt(:,1:3))) < +1 || ...
            max(max(tt(:,1:3))) > nnod )
        error('triarea:invalidInputs', ...
            'Invalid TRIA input array.') ;
    end

%--------------------------------------- compute signed area
    ev12 = pp(tt(:,2),:)-pp(tt(:,1),:) ;
    ev13 = pp(tt(:,3),:)-pp(tt(:,1),:) ;

    switch (size(pp,2))
        case +2

        area = ev12(:,1).*ev13(:,2) ...
             - ev12(:,2).*ev13(:,1) ;
        area = 0.5 * area;

        case +3

        avec = cross(ev12,ev13);
        area = sqrt(sum(avec.^2,2)) ;
        area = 0.5 * area;

        otherwise
        error('Unsupported dimension.') ;
    end

end

author: dengwirda


text:
function plot_accuracy_keyframe

color = {'r', 'y', 'g', 'b', 'm'};
leng = {'PoseCNN', 'PoseCNN+ICP', 'PoseCNN+Multiview', 'PoseCNN+ICP+Multiview', ...
    '3D Coordinate Regression'};
aps = zeros(5, 1);
lengs = cell(5, 1);
close all;

% load results
object = load('results_keyframe.mat');
distances_sys = object.distances_sys;
distances_non = object.distances_non;
rotations = object.errors_rotation;
translations = object.errors_translation;
cls_ids = object.results_cls_id;

index_plot = [4, 2, 5, 3, 1];

% read class names
fid = fopen('classes.txt', 'r');
C = textscan(fid, '%s');
classes = C{1};
classes{end+1} = 'All 21 objects';
fclose(fid);

hf = figure('units','normalized','outerposition',[0 0 1 1]);
font_size = 24;
max_distance = 0.1;

% for each class
for k = 1:numel(classes)
    index = find(cls_ids == k);
    if isempty(index)
        index = 1:size(distances_sys,1);
    end

    % distance symmetry
    subplot(2, 2, 1);
    for i = index_plot
        D = distances_sys(index, i);
        D(D > max_distance) = inf;
        d = sort(D);
        n = numel(d);
        accuracy = cumsum(ones(1, n)) / n;
        plot(d, accuracy, color{i}, 'LineWidth', 4);
        aps(i) = VOCap(d, accuracy);
        lengs{i} = sprintf('%s (%.2f)', leng{i}, aps(i) * 100);
        hold on;
    end
    hold off;
    %h = legend('network', 'refine tranlation only', 'icp', 'stereo translation only', 'stereo full', '3d coordinate');
    %set(h, 'FontSize', 16);
    h = legend(lengs(index_plot), 'Location', 'southeast');
    set(h, 'FontSize', font_size);
    h = xlabel('Average distance threshold in meter (symmetry)');
    set(h, 'FontSize', font_size);
    h = ylabel('accuracy');
    set(h, 'FontSize', font_size);
    h = title(classes{k}, 'Interpreter', 'none');
    set(h, 'FontSize', font_size);
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', font_size)

    % distance non-symmetry
    subplot(2, 2, 2);
    for i = index_plot
        D = distances_non(index, i);
        D(D > max_distance) = inf;
        d = sort(D);
        n = numel(d);
        accuracy = cumsum(ones(1, n)) / n;
        plot(d, accuracy, color{i}, 'LineWidth', 4);
        aps(i) = VOCap(d, accuracy);
        lengs{i} = sprintf('%s (%.2f)', leng{i}, aps(i) * 100);
        hold on;
    end
    hold off;
    %h = legend('network', 'refine tranlation only', 'icp', 'stereo translation only', 'stereo full', '3d coordinate');
    %set(h, 'FontSize', 16);
    h = legend(lengs(index_plot), 'Location', 'southeast');
    set(h, 'FontSize', font_size);
    h = xlabel('Average distance threshold in meter (non-symmetry)');
    set(h, 'FontSize', font_size);
    h = ylabel('accuracy');
    set(h, 'FontSize', font_size);
    h = title(classes{k}, 'Interpreter', 'none');
    set(h, 'FontSize', font_size);
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', font_size)

    % rotation
    subplot(2, 2, 3);
    for i = index_plot
        D = rotations(index, i);
        d = sort(D);
        n = numel(d);
        accuracy = cumsum(ones(1, n)) / n;
        plot(d, accuracy, color{i}, 'LineWidth', 4);
        hold on;
    end
    hold off;
    %h = legend('network', 'refine tranlation only', 'icp', 'stereo translation only', 'stereo full', '3d coordinate');
    %set(h, 'FontSize', 16);
    h = legend(leng(index_plot), 'Location', 'southeast');
    set(h, 'FontSize', font_size);
    h = xlabel('Rotation angle threshold');
    set(h, 'FontSize', font_size);
    h = ylabel('accuracy');
    set(h, 'FontSize', font_size);
    h = title(classes{k}, 'Interpreter', 'none');
    set(h, 'FontSize', font_size);
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', font_size)

    % translation
    subplot(2, 2, 4);
    for i = index_plot
        D = translations(index, i);
        D(D > max_distance) = inf;
        d = sort(D);
        n = numel(d);
        accuracy = cumsum(ones(1, n)) / n;
        plot(d, accuracy, color{i}, 'LineWidth', 4);
        hold on;
    end
    hold off;
    h = legend(leng(index_plot), 'Location', 'southeast');
    set(h, 'FontSize', font_size);
    h = xlabel('Translation threshold in meter');
    set(h, 'FontSize', font_size);
    h = ylabel('accuracy');
    set(h, 'FontSize', font_size);
    h = title(classes{k}, 'Interpreter', 'none');
    set(h, 'FontSize', font_size);
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', font_size)

    filename = sprintf('plots/%s.png', classes{k});
    hgexport(hf, filename, hgexport('factorystyle'), 'Format', 'png');
end

function ap = VOCap(rec, prec)

index = isfinite(rec);
rec = rec(index);
prec = prec(index)';

mrec=[0 ; rec ; 0.1];
mpre=[0 ; prec ; prec(end)];
for i = 2:numel(mpre)
    mpre(i) = max(mpre(i), mpre(i-1));
end
i = find(mrec(2:end) ~= mrec(1:end-1)) + 1;
ap = sum((mrec(i) - mrec(i-1)) .* mpre(i)) * 10;

author: yuxng


text:
function val = get(f, propName)
%GET   GET method for CHEBFUN3T class.
%   P = GET(F, PROP) returns the property P specified in the string PROP
%   from the CHEBFUN F. Valid entries for the string PROP are:
%    'COEFFS'
%    'VSCALE'
%    'DOMAIN'

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

% Loop through an array of CHEBFUN3 objects.
if ( numel(f) > 1 )
    val = cell(numel(f));
    for k = 1:numel(f)
        val{k} = get(f(k), propName);
    end
    return
end

% Get the properties.
switch ( propName )
    case 'coeffs'
        val = f.coeffs;
    case 'vscale'
        val = f.vscale;
    case 'domain'
        val = f.domain;
    otherwise
        error('CHEBFUN:CHEBFUN3T:get:propName', ...
            [propName,' is not a valid CHEBFUN3T property.'])
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% EEGTHRESH -  reject trials with out-of-bounds channel values within a
%                specified epoch time range.
%
% Usage:
%   >> [Iin, Iout, newsignal, elec] = eegthresh( signal, frames, ...
%                elecs, negthresh, posthresh, timerange, starttime,endtime);
%
% Required inputs:
%   signal     - 2-D data matrix [channels, frames*sweeps],
%                or 3-D data matrix [channels, frames, sweeps]
%   frames     - number of points per epoch
%   elecs      - [int vector] electrode indices to reject on
%   negthresh  - minimum rejection threshold(s) in uV. This can be an array
%                of values for individual electrodes. If fewer values than
%                electrodes, the last value is used for the remaining
%                electrodes.
%   posthresh  - maximum rejection threshold(s) in uV (same syntax as for
%                negthresh)
%   timerange  - [mintime maxtime] time range limits of the signal
%   starttime  - Starting limit (in seconds or Hz) of range to perform
%                rejection (same syntax  as negthresh)
%   endtime    - Ending limit (in seconds or Hz) of range to perform
%                rejection (same syntax  as negthresh)
%
% Outputs:
%   Iin        - Indexes of epochs accepted
%   Iout       - Indexes of epochs rejected
%   newsignal  - input data after epoch rejection
%   elec       - electrode that triggered the rejections (array of 0s
%                and 1s with the same number of columns as Iout
%                and number of rows = number of electrodes).
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: POP_EEGTHRESH, EEGLAB

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [Iin, Iout, newsignal, elec] = eegthresh( signal, pnts, electrodes, negthresh, posthresh, timerange, starttime, endtime)

if nargin < 7
	help eegthresh;
	return;
end

if starttime < timerange(1)
	disp('eegthresh: starting point out of range, adjusted');
    starttime = timerange(1);
end
if endtime > timerange(2)
	disp('eegthresh: ending point out of range, adjusted');
    endtime = timerange(2);
end
if isempty(electrodes)
    electrodes = 1:size(signal,1);
end

% complete thresholds values if necessary
%----------------------------------------
if size(posthresh,2) < size(electrodes,2)
	posthresh = [ posthresh posthresh(end)*ones(1,size(electrodes,2)-size(posthresh,2))];
end
if size(negthresh,2) < size(electrodes,2)
	negthresh = [ negthresh negthresh(end)*ones(1,size(electrodes,2)-size(negthresh,2))];
end

% complete timeselect values if necessary
%----------------------------------------
if size(starttime,2) < size(electrodes,2)
	starttime = [ starttime starttime(end)*ones(1,size(electrodes,2)-size(starttime,2))];
end
if size(endtime,2) < size(electrodes,2)
	endtime = [ endtime endtime(end)*ones(1,size(electrodes,2)-size(endtime,2))];
end

% find the maximum for each trial
%--------------------------------
sweeps = size(signal(1,:),2)/pnts;
signal = reshape(signal(electrodes,:), size(electrodes(:),1), pnts, sweeps);

% reject the selected trials
%---------------------------
elec = zeros(size(electrodes(:),1), sweeps);
allelec = zeros(1, sweeps);
for indexe = 1:size(electrodes(:),1)
	% transform the time range
	% ------------------------
	framelowlimit  = max(1,floor((starttime(indexe)-timerange(1))/(timerange(2)-timerange(1))*(pnts-1))+1);
	framehighlimit = floor((endtime(indexe)  -timerange(1))/(timerange(2)-timerange(1))*(pnts-1))+1;

	% remove outliers
	% ---------------
	sigtmp = squeeze(signal(indexe,framelowlimit:framehighlimit,:));
    if size(signal,3) == 1, sigtmp = sigtmp'; end
	sigmax = max(sigtmp, [], 1);
	sigmin = min(sigtmp, [], 1);
	elec(indexe,:) = ( sigmin < negthresh(indexe) ) | ( sigmax > posthresh(indexe) );
	allelec = allelec | elec(indexe,:);
end
Iout = find( allelec == 1 );
Iin  = find( allelec == 0 );
elec = elec(:,Iout);

% reject the selected trials
%---------------------------
newsignal = reshape(signal, size(signal,1), pnts, sweeps);
if ~isempty(Iin)
	newsignal = newsignal(:,:,Iin);
	if ndims(signal) == 2
		newsignal = newsignal(:,:);
	end
else
	newsignal = [];
end
return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [update] = polynomial_expansion_registration(...
    moving, fixed, multiModal, transformationModel, varargin)
% POLYNOMIAL_EXPANSION_REGISTRATION Estimates a displacement field using polynomial expansion
%
% INPUT ARGUMENTS
% moving                    - Moving image
% fixed                     - Fixed image
% multiModal                - Set whether to perform multi-modal or
%                             uni-modal image registration
%                             false/true
% transformationModel       - Transformation model for estimating the
%                             displacement field
%                             'translation', 'affine', 'non-rigid'
%
% OPTIONAL INPUT ARGUMENTS
% 'signalModel'             - Local signal model to use when computing the
%                             polynomial expansion transformation
%                             'linear' (deafult), 'quadratic'
%
% 'numberOfChannels'        - Number of channels to use in when computing
%                             the entropy (based on channel coding). This
%                             is only relevant if multiModal is set to
%                             true.
%                             Default value is 8
%
% OUTPUT ARGUMENTS
% update
%   displacementUpdate      - Estimated update field
%   certaintyUpdate         - Certainty related to the estimated update field
%   transformationMatrix    - Estimate transformation matrix (only if

% Copyright (c) 2012 Daniel Forsberg
% danne.forsberg@outlook.com
%
% 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/>.

%% Set up default parameters
signalModel = 'linear';
numberOfChannels = 8;

% Overwrites default parameter
for k=1:2:length(varargin)
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

%%
dims = ndims(moving);

if multiModal && strcmp(signalModel,'quadratic')
	error(['Multi-modal image registration using polynomial expansion',...
		' is only available for the linear signal model'])
end

switch dims
    case 2
		[update] = ...
			polynomial_expansion_registration2d(...
			moving, fixed,...
			'signalModel',signalModel,...
			'transformationModel',transformationModel,...
			'multiModal',multiModal,...
			'numberOfChannels',numberOfChannels);
    case 3
		[update] = ...
			polynomial_expansion_registration3d(...
			moving, fixed,...
			'signalModel',signalModel,...
			'transformationModel',transformationModel,...
			'multiModal',multiModal,...
			'numberOfChannels',numberOfChannels);
end

author: fordanic

repo: image-registration


text:
%kkcube 'Output = Input**3'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kcube.pane file
%
% Parameters:
% InputFile: i1 'Input ', required: 'Input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kkcube(i1, {'i1','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kcube - Output = Input**3
%
%  DESCRIPTION
% Output = Input*Input*Input
%
% The "Cube" operator returns the cube (value raised to the power of 3)
% of each data
% point in the input data object, \fBInput".
%
% Executing "Cube" runs the program \fIkarith2\fP
% with the -pow flag and -real 3.0 -imag 0.0.
%
%
%  "Data Type " 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "Map Data" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_1input
%
%  "Validity Mask" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "Explicit Location and Time Data" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "Failure Modes"
% .cI $DATAMANIP/repos/shared/man/sections/fail_1input
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::karith2
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkcube(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkcube(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'o', '__output'};
maxval={0,0};
minval={0,0};
istoggle=[0,0];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'karith2"  -pow'],Inputs,narglist);


text:
function tnorm = tensor_array_norm( X )
% X is a cell array of tensors

tnorm = 0;
for i = 1:length(X)
    tnorm = tnorm + norm(X{i})^2;
end
tnorm = sqrt(tnorm);

end

author: andrewssobral

repo: lrslibrary


text:
% BWBLKSLV Solves block sparse upper-triangular system.
%    y = bwblkslv(L,b) yields the same result as
%              y(L.perm,:) = L.L'\b
%    However, BWBLKSLV is faster than the built-in operator "\",
%    because it uses dense linear algebra and loop-unrolling on
%    supernodes.
%
%    Typical use, with X sparse m x m positive definite and b is m x n:
%            L = sparchol(symbchol(X),X);
%            y = bwblkslv(L,fwblkslv(L,b));
%    Then y solves X*y=b.
%
% See also symbchol, fwblkslv, mldivide, mrdivide

function y = bwblkslv(L,b) %#ok

 %
 %   This file is part of CholTool 1.00
 %   Copyright (C) 1998 Jos F. Sturm
 %   CRL, McMaster University, Canada.
 %   Supported by the Netherlands Organization for Scientific Research (NWO).
 %
 %   This program is free software; you can redistribute it and/or modify
 %   it under the terms of the GNU General Public License as published by
 %   the Free Software Foundation; either version 2 of the License, or
 %   (at your option) any later version.
 %
 %   This program is distributed in the hope that it will be useful,
 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 %   GNU General Public License for more details.
 %
 %   You should have received a copy of the GNU General Public License
 %   along with this program; if not, write to the Free Software
 %   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 %

error('At OS prompt, type "make" to create cholTool mex-files.')

repo: REST


text:
function [f,s] = is_planar(V,epsilon)
  % IS_PLANAR Return whether the 3D mesh is planar
  %
  % [f,s] = is_planar(V)
  % [f,s] = is_planar(V,epsilon)
  %
  % Inputs:
  %   V  #V x 3 matrix of vertex coordinates
  %   epsilon  optional parameter used as threshold for smallest eigen value of
  %     pca, {1e-10}
  % Outputs:
  %   f  flag whether mesh is planar
  %   s  smallest principle component variance
  %
  %

  if ~exist('epsilon','var')
    epsilon = 1e-10;
  end

  if size(V,2) == 2
    f = true;
    s = 0;
    return;
  end

  % This seems to have changed during matlab. For large V this is
  % infeasible.
  [c,l] = pcacov(V);

  s = min(l);

  f = s < epsilon;

end

author: alecjacobson

repo: gptoolbox


text:
% ASL  Read an AMPL .NL file and convert to Matlab Values
%
%   NOTE - This API has substantially changed from OPTI v1.79
%
% THIS IS A LOW LEVEL FUNCTION - USE amplRead() INSTEAD!
%
% asl uses the Netlib AMPL Solver Library to open and parse the file. The
% MEX file is a little different in that it stays open after a call to
% asl(), in order to evaluate nonlinear functions later. You must use the
% close command in order to close the mex file.
%
%   [varargout] = asl(command,varargin)
%
%   Input arguments:
%       command - a string containing the command to execute (see below)
%       varargin - command specific arguments
%
%   Return arguments:
%       varagout - command specific return variables
%
%   Command: 'open' [Open an AMPL file and read basic information + check if LP/QP/QCQP]
%       [aslprob,sizes] = asl('open',path,isNLP)
%
%       path - full path to the file on your PC
%       isNLP - 1 to skip identification of LP/QP/QCQP (assume (MI)NLP)
%
%       aslprob structure fields:
%           H  - Quadratic Objetive H matrix (sparse, if QP)
%           f  - Linear Objective f vector (if LP/QP)
%           lb - Decision variable lower bounds
%           ub - Decision variable upper bounds
%           A  - Linear Constraints A matrix (sparse, if LP/QP)
%           cl - Constraints lower bounds
%           cu - Constraints upper bounds
%           Q  - Cell array of quadratic constraint Q matrices (sparse, if QCQP)
%           l  - Matrix of quadratic constraint linear vectors (if QCQP)
%           qcind - Vector of indices of Quadratic Constraints (to remove from A and cl,cu, if QCQP)
%           x0 - initial decision variable guess
%           v0 - initial dual variable guess
%           sense - 1: minimization, -1: maximization
%           conlin - vector of constraint linearity (0 constant, 1 linear, 2 quadratic, 3 general nonlinear)
%
%       sizes - vector of problem sizes (see amplRead.m)
%
%   Command: 'isopen' [Check if AMPL ASL interface has been opened]
%       asl('isopen')
%
%   Command: 'close' [Close an opened AMPL file and destory MEX file memory]
%       asl('close')
%
%   Command: 'fun' [Evaluate objective function at x]
%       f = asl('fun',x)
%
%   Command: 'grad' [Evaluate objective gradient at x]
%       g = asl('grad',x)
%
%   Command: 'con' [Evaluate constraint function at x]
%       c = asl('con',x)
%
%   Command: 'jac' [Evaluate constraint Jacobian at x]
%       j = asl('jac',x)
%       js = asl('jac',x,1) %return sparse Jacobian
%
%   Command: 'jacstr' [Evaluate constraint Jacobian Sparsity Structure]
%       s = asl('jacstr')
%
%   Command: 'hess' [Evaluate Hessian of the Lagrangian (L = sigma*Hobj + sum(v*Hcon))]
%       h = asl('hess',x,sigma,lambda)
%       hs = asl('hess',x,sigma,lambda,1) %return sparse Hessian
%
%   Command: 'hessstr' [Evaluate Hessian Sparsity Structure]
%       s = asl('hessstr')
%
%   Copyright (C) 2011-2013 Jonathan Currie (IPL)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function planC = nrrd2cerr(filename,scanName,planC,save_flag)
% nii2cerr.m
% http://teem.sourceforge.net/nrrd/format.html#space
%
% Import .nrrd data to CERR
% Usage: planC = nrrd2cerr(filename,scanName,planC,save_flag);
%
% AI 4/10/19 for nii
% EL 06/17/21 adapt for nrrd

if ~exist('planC','var') || isempty(planC)
    planC = [];
end

if ~exist('save_flag','var') || isempty(save_flag)
    save_flag = 0;
end

% read nrrd file
% [vol3M,infoS] = nifti_read_volume(filename);
% nii = load_nii(filename);
headerInfo = nhdr_nrrd_read(filename,1);

% vol3M = nii.img;
vol3M = headerInfo.data;
infoS.Offset = [headerInfo.spaceorigin(1) headerInfo.spaceorigin(2) headerInfo.spaceorigin(3)];
infoS.PixelDimensions = [abs(headerInfo.spacedirections_matrix(1,1)) abs(headerInfo.spacedirections_matrix(2,2)) abs(headerInfo.spacedirections_matrix(3,3))];
infoS.Dimensions = headerInfo.sizes;

scanOffset = 0;
volMin = min(vol3M(:));
if volMin<0
    scanOffset = -volMin;
end

if headerInfo.spacedirections_matrix(1,1) > 0
    vol3M = flip(vol3M,1);
end

if headerInfo.spacedirections_matrix(2,2) > 0
    vol3M = flip(vol3M,2);
end

% infoS.Offset = [0 0 0];
% infoS.PixelDimensions = infoS.pixdim;
% infoS.Dimensions = infoS.dimension(2:4);

planC = mha2cerr(infoS,vol3M,scanOffset,scanName,planC,save_flag);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% Description:
%
%     Remove one term at a time from a MANCOVAN model until all p-values are
%     below the specified threshold.
%
% Syntax:
%
%     [ T, p, stats ] = mStepwise(Y, groups, [], threshold)
%     [ T, p, stats ] = mStepwise(Y, [], covariates, threshold)
%     [ T, p, stats ] = mStepwise(Y, groups, covariates, threshold)
%     [ T, p, stats ] = mStepwise(Y, groups, covariates, threshold, options)
%
%     [ T, p, stats ] = mStepwise(Y, X, terms, threshold, options)
%
% Inputs:
%
%     Y          - [ N x M ] (double) - multivariate response
%     groups     - [ N x G ] (int)    - qualitative variables
%     covariates - [ N x C ] (double) - quantitative variables
%     threshold  - [ 1 x 1 ] (double) - p-value at which to stop dropping terms
%     options    - [ 1 x P ] (cell)   - see Options
%
%     X     - [ N x T ] (double) - design matrix
%     terms - [ 1 x T ] (cell)   - model terms
%
% Outputs:
%
%     T - [ (G + C) x 1 ] (double)
%     p - [ (G + C) x 1 ] (double)
%
%     stats.U - [ N x N ] (double) - U resulting from the SVD (option 'SVD').
%     stats.S - [ N x N ] (double) - S resulting from the SVD (option 'SVD').
%     stats.V - [ M x N ] (double) - V resulting from the SVD (option 'SVD').
%
%     stats.BIC - [ 1 x P ] (double) - values of the Bayesian information
%         criterion (BIC) associated with the SVD (option 'SVD').
%
%     stats.PVE - [ N x 1 ] (double) - the percentage of variance in Y explained
%         by each consecutive column of U (option 'SVD').
%
%     stats.Terms - [ 1 x B ] (cell) - full model terms numbering groups from
%         one through size(groups, 2) and covariates from size(groups, 2) + 1
%         through size(groups, 2) + size(covariates, 2) + 1, and interactions
%         according to combinations of these numbers.
%
%     stats.X - [ N x B ] (double) - full model design matrix, including columns
%         for the grand mean, groups, covariates, and interactions.
%
%     stats.Y - [ N x M ] (double) - the multivariate response used for the
%         computation.
%
%     stats.B - [ B x P ] (double) - regression coefficients associated with the
%         full model after stepwise removal of insignificant terms.
%
%     stats.SSE - [ P x P ] (int) - full model sum of squared errors associated
%         with each column of Y.
%
%     stats.DFE - [ 1 x M ] (int) - degrees of freedom used in the computation
%         of stats.MSE.
%
%     stats.MSE - [ P x P ] (int) - full model mean squared errors associated
%         with each column of Y.
%
% Details:
%
% Options:
%
%     'group-group'         - include group-group interactions
%     'covariate-covariate' - include covariate-covariate interactions
%     'group-covariate'     - include group-covariate interactions
%     'over-determined'     - use over-determined coding for the design matrix
%     'sigma-restricted'    - use sigma-restricted coding for the design matrix
%     'SVD'                 - reduce the dimensionality of Y using an SVD
%     'verbose'             - display extra information to the command window
%
% Examples:
%
%     The following example uses a simple additive model with covariates, but no
%     interactions and avoids using the Statistics Toolbox:
%
%         n          = 100;
%         groups     = round(2 * rand(n, 2) + 0.5);
%         covariates = 10 * randn(n, 2);
%         Y          = groups + covariates + randn(n, 2);
%
%         [ T, p, stats ] = mStepwise(Y, groups, covariates, 0.10, ...
%             { 'group-group' 'covariate-covariate' 'group-covariate' 'verbose' });
%
%     For other examples, refer to mancovan.m.  The setup and syntax is exactly
%     the same in every case except for the addition of a threshold.
%
% Notes:
%
% Author(s):
%
%     William Gruner (williamgruner@gmail.com)
%
% References:
%
%     Refer to the references listed in mancovan.m.
%
% Acknowledgements:
%
%     Many thanks to Dr. Erik Erhardt and Dr. Elena Allen of the Mind Research
%     Network (www.mrn.org) for their continued collaboration.
%
% Version:
%
%     $Author: williamgruner $
%     $Date: 2010-04-15 07:46:14 -0600 (Thu, 15 Apr 2010) $
%     $Revision: 496 $

function [ T, p, stats ] = mStepwise(Y, groups, covariates, threshold, options)

    if nargin == 0
        T = BIT(); return
    end

    if ~exist('options', 'var')
        options = cell(0);
    end

    if iscell(covariates)
        X     = groups;
        terms = covariates;
    else
        [ X, terms ] = mX(groups, covariates, options);
    end

    if ~isempty(strmatch('SVD', options, 'exact'))

        [ BIC, U, S, V ] = mSVD(Y, options);

        stats.U   = U;
        stats.S   = S;
        stats.V   = V;
        stats.BIC = BIC;
        stats.PVE = cumsum(diag(S) ./ sum(diag(S)));

        b = find(diff(BIC) > 0);
        b = b(1);

        U = U(  :, 1:b);
        S = S(1:b, 1:b);
        V = V(  :, 1:b);

    else

        U = Y;

    end

    if ~isempty(strmatch('verbose', options, 'exact'))
        fprintf('\n')
    end

    while length(terms) > 1

        if ~isempty(strmatch('over-determined', options, 'exact'))
            M = X * pinv(X' * X) * X';
        else
            M = X * inv(X' * X) * X';
        end

        B = mUnique(terms(2:end));
        T = repmat(NaN, length(B), 1);
        p = repmat(NaN, length(B), 1);

        for i = 1 : length(B)

            I0 = mTerms(B{i}, terms);
            X0 = X(:, I0);

            if ~isempty(strmatch('over-determined', options, 'exact'))
                M0 = X0 * pinv(X0' * X0) * X0';
            else
                M0 = X0 * inv(X0' * X0) * X0';
            end

            if ~isempty(strmatch('verbose', options, 'exact'))
                mDispModels([ {0} B ], mUnique(terms(I0)))
            end

            [ T(i), p(i) ] = mLHT(U, X, X0, M, M0, options);

            if ~isempty(strmatch('verbose', options, 'exact'))
                fprintf('\n')
            end

        end

        if max(p) < threshold
            break
        end

        [ q, I ] = sort(p, 'descend');

        for i = 1 : length(I)
            if mIsInteractionTerm(B{I(i)})
                break
            elseif mIsMainTerm(B{I(i)})
                if ~ismember(B{I(i)}, cat(2, B{mFindInteractionTerms(B)}))
                    break
                end
            end
        end

        if q(i) < threshold

            break

        else

            if ~isempty(strmatch('verbose', options, 'exact'))

                if mIsInteractionTerm(B{I(i)})
                    fprintf('Removing %d%d ...\n', B{I(i)}(1), B{I(i)}(2))
                elseif mIsMainTerm(B{I(i)})
                    fprintf('Removing %d ...\n', B{I(i)})
                end

                fprintf('\n')

            end

            X    (:, mFindTerms(B{I(i)}, terms)) = [];
            terms(:, mFindTerms(B{I(i)}, terms)) = [];

        end

    end

    if ~isempty(strmatch('over-determined', options, 'exact'))
        stats.B = pinv(X' * X) * X' * U;
    else
        stats.B = inv(X' * X) * X' * U;
    end

    stats.Terms     = terms;
    stats.X         = X;
    stats.Y         = Y;
    stats.SSE       = U' * (eye(size(M)) - M) * U;
    stats.DFE       = size(X, 1) - rank(X);
    stats.MSE       = stats.SSE / stats.DFE;
    stats.Residuals = U - stats.X * stats.B;

function b = BIT()

    b = true;

    % Compare the results to those obtained from the previous version.

    s = load('mStepwise-BIT-1.mat');

    [ T, p, stats ] = mStepwise(s.Y, s.groups, s.covariates, 0.10, ...
        { 'sigma-restricted' 'group-group' 'covariate-covariate' 'group-covariate' 'verbose' });

    e = s.T - T;

    if any(abs(e(:)) > sqrt(eps))
        b = false;
    end

    e = s.p - p;

    if any(abs(e(:)) > sqrt(eps))
        b = false;
    end

    % Compare the results to those obtained from the previous version.

    s = load('mStepwise-BIT-2.mat');

    [ T, p, stats ] = mStepwise(s.Y, s.groups, s.covariates, 0.01, ...
        { 'sigma-restricted' 'group-group' 'covariate-covariate' 'group-covariate' 'verbose' });

    e = s.T - T;

    if any(abs(e(:)) > sqrt(eps))
        b = false;
    end

    e = s.p - p;

    if any(abs(e(:)) > sqrt(eps))
        b = false;
    end


text:
function length = circle01_length ( )

%*****************************************************************************80
%
%% CIRCLE01_LENGTH: length of the circumference of the unit circle in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real LENGTH, the length.
%
  r = 1.0;
  length = 2.0 * pi * r;

  return
end


text:
function W=computeW(geo,angles,gpuids)

    geoaux=geo;
    geoaux.sVoxel([1 2])=geo.sDetector(1)*1.1; % a bit bigger, to avoid numerical division by zero (small number)
    geoaux.sVoxel(3)=max(geo.sDetector(2),geo.sVoxel(3)); % make sure lines are not cropped. One is for when image is bigger than detector and viceversa
    geoaux.nVoxel=[2,2,2]'; % accurate enough?
    geoaux.dVoxel=geoaux.sVoxel./geoaux.nVoxel;
    W=Ax(ones(geoaux.nVoxel','single'),geoaux,angles,'Siddon','gpuids',gpuids);
    W(W<min(geo.dVoxel)/4)=Inf;
    W=1./W;
    W(W>0.1)=0.1;

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
hydro = struct();

hydro = readNEMOH(hydro,'../Coer_Comp/');
% hydro = readWAMIT(hydro,'../../WAMIT/Coer_Comp/coer_comp.out',[]);
% hydro = combineBEM(hydro); % Compare to WAMIT
hydro = radiationIRF(hydro,10,[],[],[],[]);
hydro = radiationIRFSS(hydro,[],[]);
hydro = excitationIRF(hydro,10,[],[],[],[]);
writeBEMIOH5(hydro)
plotBEMIO(hydro)

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function [ t, rank ] = subset_lex_successor ( n, t, rank )

%*****************************************************************************80
%
%% SUBSET_LEX_SUCCESSOR computes the subset lexicographic successor.
%
%  Discussion:
%
%    In the original code, there is a last element with no successor.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the number of elements in the master set.
%    N must be positive.
%
%    Input/output, integer T(N), describes a subset.  T(I) is 0 if
%    the I-th element of the master set is not in the subset, and is
%    1 if the I-th element is part of the subset.
%    On input, T describes a subset.
%    On output, T describes the next subset in the ordering.
%    If the input T was the last in the ordering, then the output T
%    will be the first.
%
%    Input/output, integer RANK, the rank.
%    If RANK = -1 on input, then the routine understands that this is
%    the first call, and that the user wishes the routine to supply
%    the first element in the ordering, which has RANK = 0.
%    In general, the input value of RANK is increased by 1 for output,
%    unless the very last element of the ordering was input, in which
%    case the output value of RANK is 0.
%

%
%  Return the first element.
%
  if ( rank == -1 )
    t(1:n) = 0;
    rank = 0;
    return
  end
%
%  Check.
%
  subset_check ( n, t );

  for i = n : -1 : 1

    if ( t(i) == 0 )
      t(i) = 1;
      rank = rank + 1;
      return
    else
      t(i) = 0;
    end

  end

  rank = 0;

  return
end


text:
function y = reshape1d(x)
% y = reshape1d(x)
% Reshapes matrix x into 1-D column vector y.
%
% $Id: reshape1d.m,v 1.1 2004/03/12 07:31:58 sayres Exp $

y = reshape(x, [prod(size(x)) 1]);

return;

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef ContrastiveLoss < dagnn.GenericLoss
  properties
    margin = 1;
  end

  methods
    function outputs = forward(obj, inputs, params)
      switch numel(inputs)
        case 3
          outputs{1} = vl_nncontrloss(inputs{:}, 'margin', obj.margin);
        case 4
          outputs{1} = vl_nncontrloss(inputs{1:3}, 'margin', inputs{4});
        otherwise
          error('Invalid number of inputs.');
      end
      if true
        isPos = inputs{3} == 1;
        % posDist = mean(norm(reshape(inputs{1}(:, :, :, isPos) - inputs{2}(:, :, :, isPos), [], sum(isPos))));
        % negDist = mean(norm(reshape(inputs{1}(:, :, :, ~isPos) - inputs{2}(:, :, :, ~isPos), [], sum(~isPos))));
        posDist = sqrt(sum((reshape(inputs{1}(:, :, :, isPos) - inputs{2}(:, :, :, isPos), [], sum(isPos))).^2,1));
        negDist = sqrt(sum((reshape(inputs{1}(:, :, :, ~isPos) - inputs{2}(:, :, :, ~isPos), [], sum(~isPos))).^2,1));
        eer = vl_eer([-1*ones(1,numel(posDist)) 1*ones(1,numel(negDist))],[posDist negDist]);
        fprintf('PD: %.2f ND: %.2f EER: %.3f ', mean(posDist), mean(negDist) ,eer);
      end

      logoutputs = outputs;
      logoutputs{1} = numel(inputs{3})*eer;
      obj.account(inputs, logoutputs);
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      switch numel(inputs)
        case 3
          [dzdx1, dzdx2] = vl_nncontrloss(inputs{:}, derOutputs{1}, 'margin', obj.margin);
        case 4
          [dzdx1, dzdx2] = vl_nncontrloss(inputs{1:3}, derOutputs{1}, 'margin', inputs{4});
        otherwise
          error('Invalid number of inputs.');
      end
      derInputs = {dzdx1, dzdx2, []};
      derParams = {} ;
    end

    function obj = ContrastiveLoss(varargin)
      obj.load(varargin{:}) ;
    end
  end
end

author: a-nagrani

repo: VGGVox

save_path: github-repos/MATLAB/a-nagrani-VGGVox


text:
function y = norm_pdf(x,mu,sigma)
%NORM_PDF Normal probability density function (pdf).
%
%   Y = NORMPDF(X,MU,SIGMA) Returns the normal pdf with
%   mean, MU, and standard deviation, SIGMA, at the values in X.
%
%   The size of X is the common size of the input arguments. A scalar input
%   functions as a constant matrix of the same size as the other inputs.
%
%   Default values for MU and SIGMA are 0 and 1 respectively.
%
% Copyright (c) 1998-2004 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin < 3,
  sigma = 1;
end

if nargin < 2;
  mu = 0;
end

if nargin < 1,
  error('Requires at least one input argument.');
end

y = exp(-0.5 * ((x-mu)./sigma).^2 -log(sigma) -log(2*pi)/2);

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


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

%*****************************************************************************80
%
%% JED_TO_YMDF_HINDU_SOLAR converts a JED to a Hindu solar YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Reingold, Nachum Dershowitz, Stewart Clamen,
%    Calendrical Calculations, II: Three Historical Calendars,
%    Software - Practice and Experience,
%    Volume 23, Number 4, pages 383-404, April 1993.
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, integer Y, M, D, real F, the YMDF date.
%

%
%  Date = JF.
%
  jf = jed;
%
%  Date = JED_EPOCH + JF
%
  jed_epoch = epoch_to_jed_hindu_solar ( );
  jf = jf - jed_epoch;
%
%  Date = JED_EPOCH + Y years + JF
%
  y = floor ( jf / year_length_hindu_solar ( ) );
  jf = jf - floor ( y * year_length_hindu_solar ( ) );
%
%  Date = JED_EPOCH + Y years + ( M - 1 ) months + JF
%
  m = 1 + floor ( jf / month_length_hindu_solar ( ) );
  jf = jf - floor ( ( m - 1 ) * month_length_hindu_solar ( ) );
%
%  Date = JED_EPOCH + Y years + ( M - 1 ) months + ( D - 1 ) days + f
%
  d = floor ( jf ) + 1;
  f = jf - ( d - 1 );

  return
end


text:
function [gKern1, gKern2, gVarmeans, gVarcovars, gInd] = rbfard2biasVardistPsi2Gradient(rbfardKern, biasKern, vardist, Z, covGrad, learnInducing)

% RBFARD2BIASVARDISTPSI2GRADIENT description.

% VARGPLVM

if nargin < 6
    learnInducing = 1;
end

% variational means
N = size(vardist.means,1);
%  inducing variables
[M Q] = size(Z);

[Psi2, Pnobias, Psi1] = rbfard2biasVardistPsi2Compute(rbfardKern, biasKern, vardist, Z);

% inverse variances
A = rbfardKern.inputScales;

% gradient wrt variance of the rbfard2 kernel
gKernvar = sum(sum(Psi2.*covGrad))/rbfardKern.variance;
% gradient for the bias parameter
gKern2 = sum(sum(Pnobias.*covGrad));
Bnm = biasKern.variance*ones(size(Psi1));
BPsi1Covg = Psi1.*(Bnm*covGrad);

%-- New: preallocation
    gVarmeans = zeros(N,Q);
    gVarcovars = zeros(N,Q);
    gInd = zeros(M,Q);
%---

% compute the gradient wrt lengthscales, variational means and variational variances
for q=1:vardist.latentDimension
%
    S_q = vardist.covars(:,q);
    Mu_q = vardist.means(:,q);
    Z_q = Z(:,q)';

    % B3_q term (without A(q); see report)
    B_q = (repmat(Mu_q,[1 M]) - repmat(Z_q,[N 1]))./repmat(A(q)*S_q + 1, [1 M]);

    % derivatives wrt variational means and inducing inputs
    tmp = (B_q.*BPsi1Covg);

    % variational means: you sum out the columns (see report)
    gVarmeans(:,q) = -A(q)*sum(tmp,2);

    % inducing inputs: you sum out the rows
    if learnInducing
        gInd(:,q) = A(q)*sum(tmp,1)';
    end

    %
    B_q = (B_q.*(repmat(Mu_q,[1 M]) - repmat(Z_q,[N 1])));

    % B1_q term (see report)
    B1_q = (repmat(S_q, [1 M]) + B_q)./repmat((A(q)*S_q + 1), [1 M]);

    % gradients wrt kernel hyperparameters (lengthscales)
    gKernlengcs(q) = -0.5*sum(sum(B1_q.*BPsi1Covg));

    % gradient wrt variational covars (diagonal covariance matrices)
    gVarcovars(:,q) = sum((BPsi1Covg./repmat((A(q)*S_q + 1), [1 M])).*(A(q)*B_q - 1),2);

    %
end
%

gKern1 = [gKernvar 2*gKernlengcs];

% gVarmeans is N x Q matrix (N:number of data, Q:latent dimension)
% this will unfold this matrix column-wise
gVarmeans = 2*gVarmeans(:)';

% gVarcovars is N x Q matrix (N:number of data, Q:latent dimension)
% this will unfold this matrix column-wise
gVarcovars = repmat(A,[N 1]).*gVarcovars;
gVarcovars = gVarcovars(:)';

% gInd is M x Q matrix (M:number of inducing variables, Q:latent dimension)
% this will unfold this matrix column-wise
gInd = 2*gInd(:)';

author: SheffieldML

repo: vargplvm


text:
function plota2south
% set the default plot direction of the x-axis

setMTEXpref('aAxisDirection','south');
setMTEXpref('bAxisDirection','');
setMTEXpref('aStarAxisDirection','');

author: mtex-toolbox

repo: mtex

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


text:
function test_int_test ( )

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

  test_int_test01 ( )
  test_int_test02 ( )
  test_int_test03 ( )
  test_int_test04 ( )
  test_int_test05 ( )
  test_int_test06 ( )
  test_int_test07 ( )
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_INT_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( )

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pot = cg_mom_to_can(pot)
% CG_MOM_TO_CAN Convert a CG potential from moment to canonical form, if necessary.
% pot = cg_mom_to_can(pot)

if pot.subtype ~= 'c'
  for i=1:pot.dsize
    pot.can{i} = mpot_to_cpot(pot.mom{i});
  end
  pot.subtype = 'c';
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%% Lane Marking Identification
% One component of future luxury automobile safety
% systems is warning drivers that they are drifting
% between lanes.  To do this, it is first necessary
% to identify where the lanes are.  In this example,
% we will look at identifying the lanes from a image
% of a road.
%
% Copyright 2007-2013 MathWorks, Inc.
%

%% Setup
clear

%% Initialize Objects
% Video Reader

hvfr = vision.VideoFileReader('viplanedeparture.avi', ...
    'VideoOutputDataType', 'uint8');

%% Configure Video Players
hVidSource = vision.VideoPlayer('Name', 'Source'      );
hVidLanes  = vision.VideoPlayer('Name', 'Verification');

%% Replay Loop
NumLoops = 1;
for loopIdx = 1:NumLoops

    % Loop through video frames
    while ~isDone(hvfr)
        % Read Frame and display
        frameRGB = step(hvfr);
        step(hVidSource, frameRGB);

        % Call algorithm
        laneVertices = lanemarking_Algorithm(frameRGB);

        % Display lanes
        displayLanes(frameRGB, laneVertices);

    % End While
    end

    reset(hvfr);
%% end loop
end

%% Release resources
release(hvfr   );

author: Sable

repo: mcbench-benchmarks


text:
function gpcf = gpcf_cat(varargin)
%GPCF_CAT  Create a categorical covariance function
%
%  Description
%    GPCF = GPCF_CAT('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a categorical covariance function structure in
%    which the named parameters have the specified values. Any
%    unspecified parameters are set to default values.
%
%    GPCF = GPCF_CAT(GPCF,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a covariance function structure with the named
%    parameters altered with the specified values.
%
%    Categorical covariance function returns correlation 1 if input
%    values X_i and X_j are equal and 0 otherwise.
%
%    Parameters for categorical covariance function
%      selectedVariables - a vector defining which inputs are used
%
%  See also
%    GP_SET, GPCF_*, PRIOR_*, MEAN_*
%
% Copyright (c) 2007-2010 Jarno Vanhatalo
% Copyright (c) 2008-2010 Jaakko Riihimäki
% Copyright (c) 2010 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'GPCF_CAT';
  ip.addOptional('gpcf', [], @isstruct);
  ip.addParamValue('selectedVariables',[], ...
                   @(x) isempty(x) || (isvector(x) && all(x>0)));
  ip.parse(varargin{:});
  gpcf=ip.Results.gpcf;

  if isempty(gpcf)
    % Initialize a covariance function structure
    init=true;
    gpcf.type = 'gpcf_cat';
  else
    % Modify a covariance function structure
    if ~isfield(gpcf,'type') && ~isequal(gpcf.type,'gpcf_cat')
      error('First argument does not seem to be a valid covariance function structure')
    end
    init=false;
  end

  if ~ismember('selectedVariables',ip.UsingDefaults)
    if ~isempty(ip.Results.selectedVariables)
      gpcf.selectedVariables = ip.Results.selectedVariables;
    elseif isfield(gpcf,'selectedVariables')
      gpcf=rmfield(gpcf,'selectedVariables');
    end
  end
  if init
    gpcf.fh.pak = @gpcf_cat_pak;
    gpcf.fh.unpak = @gpcf_cat_unpak;
    gpcf.fh.lp = @gpcf_cat_lp;
    gpcf.fh.lpg = @gpcf_cat_lpg;
    gpcf.fh.cfg = @gpcf_cat_cfg;
    gpcf.fh.ginput = @gpcf_cat_ginput;
    gpcf.fh.cov = @gpcf_cat_cov;
    gpcf.fh.trcov  = @gpcf_cat_trcov;
    gpcf.fh.trvar  = @gpcf_cat_trvar;
    gpcf.fh.recappend = @gpcf_cat_recappend;
  end

end

function [w,s,h] = gpcf_cat_pak(gpcf, w)
%GPCF_CAT_PAK  Combine GP covariance function parameters into
%              one vector.
%
%  Description
%    W = GPCF_CAT_PAK(GPCF) takes a covariance function
%    structure GPCF and combines the covariance function
%    parameters and their hyperparameters into a single row
%    vector W. This is a mandatory subfunction used for
%    example in energy and gradient computations.
%
%       w = []
%
%  See also
%    GPCF_CAT_UNPAK

  w = []; s = {}; h=[];
end

function [gpcf, w] = gpcf_cat_unpak(gpcf, w)
%GPCF_CAT_UNPAK  Sets the covariance function parameters into
%                the structure
%
%  Description
%    [GPCF, W] = GPCF_CAT_UNPAK(GPCF, W) takes a covariance
%    function structure GPCF and a parameter vector W, and
%    returns a covariance function structure identical to the
%    input, except that the covariance parameters have been set
%    to the values in W. Deletes the values set to GPCF from W
%    and returns the modified W. This is a mandatory subfunction
%    used for example in energy and gradient computations.
%
%    Assignment is inverse of
%       w = []
%
%  See also
%   GPCF_CAT_PAK

end

function lp = gpcf_cat_lp(gpcf)
%GPCF_CAT_LP  Evaluate the energy of prior of covariance function parameters
%
%  Description
%    LP = GPCF_CAT_LP(GPCF) takes a covariance function
%    structure GPCF and returns log(p(th)), where th collects the
%    parameters.  This is a mandatory subfunction used for example
%    in energy computations.
%
%  See also
%    GPCF_CAT_PAK, GPCF_CAT_UNPAK, GPCF_CAT_LPG, GP_E

  lp = 0;

end

function lpg = gpcf_cat_lpg(gpcf)
%GPCF_CAT_LPG  Evaluate gradient of the log prior with respect
%               to the parameters.
%
%  Description
%    LPG = GPCF_CAT_LPG(GPCF) takes a covariance function
%    structure GPCF and returns LPG = d log (p(th))/dth, where th
%    is the vector of parameters.  This is a mandatory subfunction
%    used for example in gradient computations.
%
%  See also
%    GPCF_CAT_PAK, GPCF_CAT_UNPAK, GPCF_CAT_LP, GP_G

  lpg = [];

end

function DKff = gpcf_cat_cfg(gpcf, x, x2, mask, i1)
%GPCF_CAT_CFG  Evaluate gradient of covariance function
%              with respect to the parameters.
%
%  Description
%    DKff = GPCF_CAT_CFG(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to th (cell array with matrix elements). This is a
%    mandatory subfunction used in gradient computations.
%
%    DKff = GPCF_CAT_CFG(GPCF, X, X2) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X2) with
%    respect to th (cell array with matrix elements). This subfunction
%    is needed when using sparse approximations (e.g. FIC).
%
%    DKff = GPCF_CAT_CFG(GPCF, X, [], MASK) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the diagonal of gradients of covariance matrix
%    Kff = k(X,X2) with respect to th (cell array with matrix
%    elements). This subfunction is needed when using sparse
%    approximations (e.g. FIC).
%
%    DKff = GPCF_CAT_CFG(GPCF, X, X2, MASK, i) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2), or k(X,X) if X2 is empty, with respect to ith
%    hyperparameter. This subfunction is needed when using
%    memory save option in gp_set.
%
%  See also
%    GPCF_CAT_PAK, GPCF_CAT_UNPAK, GPCF_CAT_LP, GP_G

  if nargin==5
    % Use memory save option
    if i1==0
      % Return number of hyperparameters
      DKff=0;
      return
    end
  end
  DKff = {};

end

function [DKff, lpg]  = gpcf_cat_ginput(gpcf, x, x2, i1)
%GPCF_CAT_GINPUT  Evaluate gradient of covariance function with
%                 respect to x.
%
%  Description
%    DKff = GPCF_CAT_GINPUT(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to X (cell array with matrix elements). This subfunction
%    is needed when computing gradients with respect to inducing
%    inputs in sparse approximations.
%
%    DKff = GPCF_CAT_GINPUT(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors
%    and returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to X (cell array with matrix elements).
%    This subfunction is needed when computing gradients with
%    respect to inducing inputs in sparse approximations.
%
%    DKff = GPCF_CAT_GINPUT(GPCF, X, X2, i) takes a covariance
%    function structure GPCF, a matrix X of input vectors
%    and returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2), or k(X,X) if X2 is empty, with respect to ith
%    covariate in X. This subfunction is needed when using
%    memory save option in gp_set.
%
%  See also
%   GPCF_CAT_PAK, GPCF_CAT_UNPAK, GPCF_CAT_LP, GP_G

  [n, m] =size(x);

  if nargin==4
    % Use memory save option
    if i1==0
      % Return number of covariates
      if isfield(gpcf,'selectedVariables')
        DKff=length(gpcf.selectedVariables);
      else
        DKff=m;
      end
      return
    end
  end

  if nargin == 2 || isempty(x2)
    ii1 = 0;
    for j = 1:n
      for i=1:m
        ii1 = ii1 + 1;
        DKff{ii1} = zeros(n);
        lpg(ii1) = 0;
      end
    end
  elseif nargin == 3 || nargin == 4
    ii1 = 0;
    for j = 1:n
      for i=1:m
        ii1 = ii1 + 1;
        DKff{ii1} = zeros(n, size(x2,1));
        lpg(ii1) = 0;
      end
    end
  end
end

function C = gpcf_cat_cov(gpcf, x1, x2, varargin)
%GP_CAT_COV  Evaluate covariance matrix between two input vectors
%
%  Description
%    C = GP_CAT_COV(GP, TX, X) takes in covariance function of a
%    Gaussian process GP and two matrixes TX and X that contain
%    input vectors to GP. Returns covariance matrix C. Every
%    element ij of C contains covariance between inputs i in TX
%    and j in X. This is a mandatory subfunction used for example
%    in prediction and energy computations.
%
%  See also
%    GPCF_CAT_TRCOV, GPCF_CAT_TRVAR, GP_COV, GP_TRCOV

  if isempty(x2)
    x2=x1;
  end
  [n1,m1]=size(x1);
  [n2,m2]=size(x2);

  if m1~=m2
    error('the number of columns of X1 and X2 has to be same')
  end

  C=repmat(true,n1,n2);
  if isfield(gpcf, 'selectedVariables')
    for j = 1:length(gpcf.selectedVariables)
      jj=gpcf.selectedVariables(j);
      C = C & bsxfun(@eq,x1(:,jj),x2(:,jj)');
    end
  else
    for j = 1:m1
      C = C & bsxfun(@eq,x1(:,j),x2(:,j)');
    end
  end
  C=double(C);
  C = sparse(C);
end

function C = gpcf_cat_trcov(gpcf, x)
%GP_CAT_TRCOV  Evaluate training covariance matrix of inputs
%
%  Description
%    C = GP_CAT_TRCOV(GP, TX) takes in covariance function of a
%    Gaussian process GP and matrix TX that contains training
%    input vectors. Returns covariance matrix C. Every element ij
%    of C contains covariance between inputs i and j in TX. This
%    is a mandatory subfunction used for example in prediction and
%    energy computations.
%
%  See also
%    GPCF_CAT_COV, GPCF_CAT_TRVAR, GP_COV, GP_TRCOV

  [n,m]=size(x);

  C=repmat(true,n,n);
  if isfield(gpcf, 'selectedVariables')
    for j = 1:length(gpcf.selectedVariables)
      jj=gpcf.selectedVariables(j);
      C = C & bsxfun(@eq,x(:,jj),x(:,jj)');
    end
  else
    for j = 1:m
      C = C & bsxfun(@eq,x(:,j),x(:,j)');
    end
  end
  C=double(C);
  C = sparse(C);
end

function C = gpcf_cat_trvar(gpcf, x)
%GP_CAT_TRVAR  Evaluate training variance vector
%
%  Description
%    C = GP_CAT_TRVAR(GPCF, TX) takes in covariance function of a
%    Gaussian process GPCF and matrix TX that contains training
%    inputs. Returns variance vector C. Every element i of C
%    contains variance of input i in TX. This is a mandatory
%    subfunction used for example in prediction and energy
%    computations.
%
%  See also
%    GPCF_CAT_COV, GP_COV, GP_TRCOV

  [n,m]=size(x);
  C=ones(n,1);

end

function reccf = gpcf_cat_recappend(reccf, ri, gpcf)
%RECAPPEND  Record append
%
%  Description
%    RECCF = GPCF_CAT_RECAPPEND(RECCF, RI, GPCF) takes a
%    covariance function record structure RECCF, record index RI
%    and covariance function structure GPCF with the current MCMC
%    samples of the parameters. Returns RECCF which contains
%    all the old samples and the current samples from GPCF.
%    This subfunction is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC and GP_MC -> RECAPPEND

  if nargin == 2
    % Initialize the record
    reccf.type = 'gpcf_cat';

    % Set the function handles
    reccf.fh.pak = @gpcf_cat_pak;
    reccf.fh.unpak = @gpcf_cat_unpak;
    reccf.fh.lp = @gpcf_cat_lp;
    reccf.fh.lpg = @gpcf_cat_lpg;
    reccf.fh.cfg = @gpcf_cat_cfg;
    reccf.fh.cov = @gpcf_cat_cov;
    reccf.fh.trcov  = @gpcf_cat_trcov;
    reccf.fh.trvar  = @gpcf_cat_trvar;
    reccf.fh.recappend = @gpcf_cat_recappend;
  else
    % Append to the record
    if isfield(gpcf, 'selectedVariables')
      reccf.selectedVariables = gpcf.selectedVariables;
    end
  end
end


text:
%SpatialMomentum Spatial momentum class
%
% Concrete subclass of SpatialF6 and represents the
% translational and rotational momentum of a rigid-body moving in 3D space.
%
%          SpatialVec6 (abstract handle class)
%            |
%            +--- SpatialM6 (abstract)
%            |     |
%            |     +---SpatialVelocity
%            |     +---SpatialAcceleration
%            |
%            +---SpatialF6 (abstract)
%                 |
%                 +---SpatialForce
%                 +---SpatialMomentum
%
% Methods::
%  SpatialMomentum    ^constructor invoked by subclasses
%  new                construct new concrete class of same type
%  double             ^convert to a 6xN double
%  char               ^convert to string
%  cross              ^^cross product
%  display            ^display in human readable form
%
% Operators::
%  +          ^add spatial vectors of the same type
%  -          ^subtract spatial vectors of the same type
%  -          ^unary minus of spatial vectors
%
% Notes:
%  - ^ is inherited from SpatialVec6.
%  - ^^ is inherited from SpatialM6.
%
% References::
%
%  - Robot Dynamics Algorithms, R. Featherstone, volume 22,
%    Springer International Series in Engineering and Computer Science,
%    Springer, 1987.
%  - A beginner's guide to 6-d vectors (part 1), R. Featherstone,
%    IEEE Robotics Automation Magazine, 17(3):83-94, Sep. 2010.
%
% See also SpatialVec6, SpatialF6, SpatialForce.

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

classdef SpatialMomentum < SpatialF6
        methods
        function n = new(a, val)
            %SpatialMomentum.new Construct a new object of the same type
            %
            % A2 = A.new(X) creates a new object of the same type as A, with the value
            % X (6x1).
            %
            % Notes::
            %  - Serves as a dynamic constructor.
            %  - This method is polymorphic across all SpatialVec6 derived classes, and
            %    allows easy creation of a new object of the same class as an existing
            %    one without needing to explicitly determine its type.
            n = SpatialMomentum(val);
        end
    end
end

author: petercorke

repo: spatialmath-matlab


text:
%% triangleRayIntersection
% Below is a demonstration of the features of the |triangleRayIntersection| function

%%
clear; close all; clc;

%% Syntax
% |[V_intersect,L_intersect,T] = triangleRayIntersection (V_ori,R,V,F,optStruct);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*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:
% Monopole Point Source In A Homogeneous Propagation Medium Example
%
% This example provides a simple demonstration of using k-Wave for the
% simulation and detection of a time varying pressure source within a
% two-dimensional homogeneous propagation medium.  It builds on the
% Homogeneous Propagation Medium and Recording The Particle Velocity
% examples.
%
% author: Bradley Treeby
% date: 2nd December 2009
% last update: 2nd September 2012
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

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

clear all;

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
Nx = 128;           % number of grid points in the x (row) direction
Ny = 128;           % number of grid points in the y (column) direction
dx = 50e-3/Nx;    	% grid point spacing in the x direction [m]
dy = dx;            % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500;  % [m/s]
medium.alpha_coeff = 0.75;  % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;

% create the time array
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);

% define a single source point
source.p_mask = zeros(Nx, Ny);
source.p_mask(end - Nx/4, Ny/2) = 1;

% define a time varying sinusoidal source
source_freq = 0.25e6;   % [Hz]
source_mag = 2;         % [Pa]
source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array);

% filter the source to remove high frequencies not supported by the grid
source.p = filterTimeSeries(kgrid, medium, source.p);

% define a single sensor point
sensor.mask = zeros(Nx, Ny);
sensor.mask(Nx/4, Ny/2) = 1;

% define the acoustic parameters to record
sensor.record = {'p', 'p_final'};

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);

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

% plot the final wave-field
figure;
imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data.p_final + source.p_mask + sensor.mask, [-1 1]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;

% plot the simulated sensor data
figure;
[t_sc, scale, prefix] = scaleSI(max(kgrid.t_array(:)));

subplot(2, 1, 1), plot(kgrid.t_array*scale, source.p, 'k-');
xlabel(['Time [' prefix 's]']);
ylabel('Signal Amplitude');
axis tight;
title('Input Pressure Signal');

subplot(2, 1, 2), plot(kgrid.t_array*scale, sensor_data.p, 'r-');
xlabel(['Time [' prefix 's]']);
ylabel('Signal Amplitude');
axis tight;
title('Sensor Pressure Signal');

repo: Aero-matlab


text:
function chrpak_test ( )

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

  chrpak_test001 ( );
  chrpak_test005 ( );
  chrpak_test006 ( );

  chrpak_test011 ( );
  chrpak_test016 ( );

  chrpak_test021 ( );
  chrpak_test022 ( );
  chrpak_test026 ( );
  chrpak_test029 ( );

  chrpak_test046 ( );

  chrpak_test051 ( );
  chrpak_test054 ( );
  chrpak_test056 ( );
  chrpak_test057 ( );

  chrpak_test065 ( );
  chrpak_test066 ( );
  chrpak_test067 ( );

  chrpak_test083 ( );

  chrpak_test090 ( );
  chrpak_test091 ( );

  chrpak_test1015 ( );
  chrpak_test102 ( );
  chrpak_test105 ( );

  chrpak_test115 ( );
  chrpak_test1155 ( );

  chrpak_test1225 ( );
  chrpak_test1227 ( );
  chrpak_test1255 ( );
  chrpak_test129 ( );

  chrpak_test137 ( );
  chrpak_test138 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHRPAK_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
%%
% Get input file without braces.
% Outputs grid info:
params = read_input_file('input2d.IB2D');
Nx = params{find(strcmp({params{:,1}},'Nx')),2}
Ny = params{find(strcmp({params{:,1}},'Ny')),2}
Lx = params{find(strcmp({params{:,1}},'Lx')),2}
Ly = params{find(strcmp({params{:,1}},'Ly')),2}
supp = params{find(strcmp({params{:,1}},'supp')),2}
d_springs = params{find(strcmp({params{:,1}},'damped_springs')),2}
str_name = params{find(strcmp({params{:,1}},'string_name')),2}

params

%%
% Get input file with braces.
% Outputs grid info
params = read_input_file('input2d.IB2D');
grid_Info = params{find(strcmp({params{:,1}},'GridParameters')),2}

grid_Info

author: nickabattista


text:
function f = ecc2flat(e)
%ECC2FLAT   Convert the eccentricity of an ellipsoid to its flattening
%
%  F = ECC2FLAT(E) returns the flattening given the eccentricity.
%
%   See also FLAT2ECC.

  e2 = e.^2;
  f = e2 ./ (1 + sqrt(1 - e2));
end