text:
%--------------------------------------------------------------------------------------------------------
% The system is created based on the principles described in the following paper
% Jimmy SJ. Ren and Li Xu, "On Vectorization of Deep Convolutional Neural Networks for Vision Tasks",
% The 29th AAAI Conference on Artificial Intelligence (AAAI-15). Austin, Texas, USA, January 25-30, 2015
% email: jimmy.sj.ren@gmail.com
%--------------------------------------------------------------------------------------------------------
addpath applications/image_denoise/
addpath applications/image_denoise/utility/
addpath utils/
addpath cuda/
addpath mem/
addpath layers/
addpath layers_adapters/
addpath optimization/
addpath pipeline/
addpath data/

clearvars -global config;
clearvars -global mem;
clear;
global config mem;
denoise_configure();
init(0);

load('data/denoise/val/val_1');
perm = randperm(size(test_samples, 4));
test_samples = test_samples(:,:,:,perm);
test_labels = test_labels(:,:,:,perm);
test_samples = config.NEW_MEM(test_samples(:,:,:,1:1000));
test_labels = config.NEW_MEM(test_labels(:,:,:,1:1000));

count = 0;
cost_avg = 0;
epoc = 0;
points_seen = 0;
display_points = 1000;
save_points = 10000;
max_grad = 1;
fprintf('%s\n', datestr(now, 'dd-mm-yyyy HH:MM:SS FFF'));
for pass = 1:10
    for p = 1:11
        load(strcat('data/denoise/train/patches_', num2str(p), '.mat'));
        perm = randperm(10000);
        samples = samples(:,:,:,perm);
        labels = labels(:,:,:,perm);
        train_imgs = config.NEW_MEM(samples);
        train_labels = config.NEW_MEM(labels);

        for i = 1:size(train_labels, 4) / config.batch_size
            points_seen = points_seen + config.batch_size;
            in = train_imgs(:,:,:,(i-1)*config.batch_size+1:i*config.batch_size);
            out = train_labels(:,:,:,(i-1)*config.batch_size+1:i*config.batch_size);
            out = out((size(in, 1) - config.output_size(1)) / 2 + 1:(size(in, 1) - config.output_size(1)) / 2 + config.output_size(1), ...
                      (size(in, 2) - config.output_size(2)) / 2 + 1:(size(in, 2) - config.output_size(2)) / 2 + config.output_size(2), :, :);

            % operate the training pipeline
            op_train_pipe(in, out);
            % update the weights
            config.UPDATE_WEIGHTS();

            if(cost_avg == 0)
                cost_avg = config.cost;
            else
                cost_avg = (cost_avg + config.cost) / 2;
            end

            % display point
            if(mod(points_seen, display_points) == 0)
                count = count + 1;
                fprintf('%d ', count);
            end
            % save point
            if(mod(points_seen, save_points) == 0)
                fprintf('\n%s', datestr(now, 'dd-mm-yyyy HH:MM:SS FFF'));
                epoc = epoc + 1;
                test_cost = 0;
                for t = 1:size(test_samples, 4) / config.batch_size
                    t_label = test_labels(:,:,:,(t-1)*config.batch_size+1:t*config.batch_size);
                    t_label = config.NEW_MEM(t_label((size(in, 1) - config.output_size(1)) / 2 + 1:(size(in, 1) - config.output_size(1)) / 2 + config.output_size(1), ...
                                            (size(in, 2) - config.output_size(2)) / 2 + 1:(size(in, 2) - config.output_size(2)) / 2 + config.output_size(2), :));

                    op_test_pipe(test_samples(:,:,:,(t-1)*config.batch_size+1:t*config.batch_size), t_label);
                    test_out = gather(mem.output);
                    test_cost = test_cost + config.cost;
                end
                test_cost = test_cost / size(test_samples, 4);
                fprintf('\nepoc %d, training avg cost: %f, test avg cost: %f\n', epoc, cost_avg, test_cost);

                save_weights(strcat('applications/image_denoise/results/epoc', num2str(epoc), '.mat'));
                cost_avg = 0;
            end
        end
    end
end

author: jimmy-ren


text:
function mergedData = mergeTextDataAndData(textdata, data, headings)
% Merges textdata and data imported from .xls file assuming that the first
% row of textdata is column headings
%
% USAGE:
%
%    mergedData = mergeTextDataAndData(textdata, data, headings)
%
% INPUTS:
%    textdata:      cell array from .xls import
%    data:          matrix with numeric data from .xls import
%
% OPTIONAL INPUT:
%    headings:      {(1), 0}, zero if no column headings
%
% OUTPUT:
%    mergedData:    merged cell array with all data from .xls import
%
% .. Author: - Ronan Fleming 29/10/2008

if ~exist('headings','var'); %check for headings
    headings=1;
    fprintf('%s\n','Assuming there were headings in the original xls file');
end
if headings==1
    start=2;
else
    start=1;
end

%checking for exceptions in otherwise data columns
fprintf('%s\n','Checking for exceptions in otherwise numeric columns, e.g. NaN');
fprintf('%s\n','Checking for exceptions in otherwise numeric columns, e.g. -');
fprintf('%s\n','Checking for exceptions in otherwise numeric columns, e.g. Not calculated');
fprintf('%s\n','Checking for exceptions in otherwise numeric columns, e.g. #N/A');
fprintf('%s\n','Checking for empty cells in otherwise numeric columns, e.g. []');
[ylt,xlt]=size(textdata);
for y=1:ylt
    for x=1:xlt
        if strcmp(textdata{y,x},'NaN')
            textdata{y,x}='';
        end
    end
end
for y=1:ylt
    for x=1:xlt
        if strcmp(textdata{y,x},'-')
            textdata{y,x}='';
        end
    end
end
for y=1:ylt
    for x=1:xlt
        if strcmp(textdata{y,x},'Not calculated')
            textdata{y,x}='';
        end
    end
end
for y=1:ylt
    for x=1:xlt
        if strcmp(textdata{y,x},'#N/A')
            textdata{y,x}='';
        end
    end
end
%replace empty cells with blank string
for y=1:ylt
    for x=1:xlt
        if isempty(textdata{y,x})
            textdata{y,x}='';
        end
    end
end
%preallocate
mergedData=textdata;
dataCol=1;
beginData=0;
for x=1:xlt
    %check for blank column
    if min(strcmp(textdata(start:ylt,x),''))==1
        beginData=1;
        fprintf('%s\n',['Merging numerical data into column ' int2str(x)]);
        for y=start:ylt
            if (y-start+1)>size(data,1)
                %if the last rows are NaN then this is needed as data will
                %be too short
                mergedData{y,x}=NaN;
            else
                mergedData{y,x}=data(y-start+1,dataCol);
            end
        end
        dataCol=dataCol+1;
        if dataCol>size(data,2)
            break
        end
    else
        %data can have empty columns corresponding to textdata columns
        if beginData==1
            dataCol=dataCol+1;
        end
    end
end
%data column might be the last column of original xls file
if dataCol<size(data,2)+1
    p=1;
    while dataCol<size(data,2)+1
        for y=start:ylt
            mergedData{y,x+p}=data(y-start+1,dataCol);
        end
        p=p+1;
        dataCol=dataCol+1;
    end
end

% for x=1:137 gMW{x,1}=iCore.genes{x}; for y=1:5300 if
% strncmp(iCore.genes{x},mergedData{y,1},length(iCore.genes{x})) gMW{x,2}=d(y); end; end; end;

author: opencobra

repo: cobratoolbox


text:
function T = tprod1(S, U, n)
%TPROD1 Tensor Product of a tensor and a matrix
%	T = TPROD1(S, U, n)
%
%	S  - tensor (multidimensional array)
%	U  - matrix compatible with the nth size of S ie. size(S,n)==size(U,1)
%	n  - execute tprod in this dimension
%
%	T  - result of the product
%
%	eg. tprod1(ones(2,3,4), [1 0 0; 0 1 0], 2)
%
%	See also TPROD.

% TODO: n > dimensions of S -> ndim_expand

siz = size(S);
siz(n) = size(U,1);
H = ndim_unfold(S, n);
T = ndim_fold(U*H, n, siz);

author: Sable

repo: mcbench-benchmarks


text:
% getFolder
%
% Similar to uigetdir except that it will remember the location of the
% last opened folder.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use plotAverageWaveforms
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Kian Torab
%   ktorab@blackrockmicro.com
%   Blackrock Microsystems
%   Salt Lake City, UT
%
%   Version 1.0.0.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dataFolder = getFolder

settingFileFullPath = getSettingFileFullPath('getFolder');

%% Opens the getFolder.ini file to see what the last accessed folder was
if exist(settingFileFullPath, 'file') == 2
    settingsFID = fopen(settingFileFullPath, 'r');
    defaultOpenLocation = fscanf(settingsFID, '%200c');
    fclose(settingsFID);
else
    defaultOpenLocation = [];
end

%% Gets a folder by opening the last accessed folder as the last location
dataFolder = uigetdir(defaultOpenLocation);

%% Saves the last opened folder in the getFolder.ini file for later use
if ischar(dataFolder)
    settingsFID = fopen(settingFileFullPath, 'w');
    fprintf(settingsFID, '%s', dataFolder);
    fclose(settingsFID);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function im2ucms_all(database,gt_set)

if nargin==0
   database = 'pascal2012';
end
if nargin<2
   gt_set = 'train2012';
end

res_id = 'sf_mUCM';

% Multiple scales to work on, finest first
scales = [2, 1, 0.5];

% Result dir
res_dir = fullfile(root_dir,'datasets',database,res_id);

% Check that results dir exist
if ~exist(fullfile(res_dir,'multi'),'dir')
    mkdir(fullfile(res_dir,'multi'))
end
for ii=1:length(scales)
    if ~exist(fullfile(res_dir,['scale_',sprintf('%1.2f',scales(ii))]),'dir')
        mkdir(fullfile(res_dir,['scale_',sprintf('%1.2f',scales(ii))]))
    end
end

% Which images to process
im_ids = database_ids(database,gt_set);

% Load pre-trained model for contour detector
model = loadvar(fullfile(root_dir, 'datasets', 'models', 'sf_modelFinal.mat'),'model');

% Sweep all images in parallel
num_workers = 32;
if verLessThan('matlab','8.3') %Pre-2014A code
    if matlabpool('size') == 0
        matlabpool(num_workers)
    end
    parfor ii=1:length(im_ids)
        im = get_image(database,im_ids{ii});

        % Check if it is not computed already
        if ~exist(fullfile(res_dir,'multi',[im_ids{ii} '.mat']),'file')

            % Call the actual code
            [ucm2,ucms] = img2ucms(im, model, scales);

            % Store ucms at each scale separately
            parsave(fullfile(res_dir,'multi',[im_ids{ii} '.mat']),ucm2)
            for jj=1:length(scales)
                parsave(fullfile(res_dir,['scale_',sprintf('%1.2f',scales(jj))],[im_ids{ii} '.mat']),ucms(:,:,jj))
            end
        end
    end
    matlabpool close
else                    %2015A or later
    if isempty(gcp('nocreate'))
        parpool('local', num_workers)
    end
    parfor ii=1:length(im_ids)
        im = get_image(database,im_ids{ii});

        % Check if it is not computed already
        if ~exist(fullfile(res_dir,'multi',[im_ids{ii} '.mat']),'file')

            % Call the actual code
            [ucm2,ucms] = img2ucms(im, model, scales);

            % Store ucms at each scale separately
            parsave(fullfile(res_dir,'multi',[im_ids{ii} '.mat']),ucm2)
            for jj=1:length(scales)
                parsave(fullfile(res_dir,['scale_',sprintf('%1.2f',scales(jj))],[im_ids{ii} '.mat']),ucms(:,:,jj))
            end
        end
    end
    delete(gcp('nocreate'))
end
end

function parsave(res_file, ucm2) %#ok<INUSD>
    save(res_file,'ucm2');
end

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function [c,fc] = FInvSABR4_2(x,x_sabr,y_sabr)
% This function computes the inverse of a SABR cumulative distribution
% Since it is based on FSABR2 it is very fast and works for vectors!
f2 = @(t) (x - FSABR2_2(t,x_sabr,y_sabr));

% Use simple bisection
a = 0.00001 * ones(1,length(x)); % left starting points
b = 1*ones(1,length(x));        % right starting points

eps = 1e-3;                     % accuracy level
k = 0;                          % init iteration counter

c=(b+a)/2;
fc = f2(c);
%Ind = -eps <= fc & fc <= eps;
while(k <= 200)
    k = k+1;                            % increase iteration counter
    if( (-eps <= min(fc)) & (max(fc) <= eps))
        break;                          % end if all values have been found
    else
        b(fc<0) = 0.5*(b(fc<0)+a(fc<0));               % update right boundary
        a(fc>=0) = 0.5*(b(fc>=0)+a(fc>=0));            % update left boundary
    end
    c = 0.5*(b+a);
    fc = f2(c);                         % evaluate function at new c
end

end


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function theStruct = MU_parseXML(dnode)
%DoParseXML(XML_filename)
%convert XML file to a MATLAB structure.

try
   tree = xmlread(dnode);
catch
   error('Failed to read XML file %s.',dnode);
end

% Recurse over child nodes. This could run into problems
% with very deeply nested trees.

theStruct = parseChildNodes(tree);
theStruct.current=1;

% ----- Subfunction PARSECHILDNODES -----
function children = parseChildNodes(theNode)
% Recurse over node children.
children = [];
if theNode.hasChildNodes
   childNodes = theNode.getChildNodes;
   numChildNodes = childNodes.getLength;
   allocCell = cell(1, numChildNodes);
   children = struct('Name', allocCell,...
                     'Attributes', allocCell,...
                     'Data', allocCell,...
                     'Children', allocCell,...
                     'current', allocCell);
    i=0;
    for count = 1:numChildNodes
        theChild = childNodes.item(count-1);
        c = makeStructFromNode(theChild);
        if isempty(c)
            i=i+1;
        else
            children(count-i) = c;
        end
    end
    if i>0
        children(numChildNodes-i+1:numChildNodes)=[];
    end
end

% ----- Subfunction MAKESTRUCTFROMNODE -----
function nodeStruct = makeStructFromNode(theNode)
% Create structure of node info.

PC=parseChildNodes(theNode);
nodeStruct = struct('Name', char(theNode.getNodeName),...
                    'Attributes', parseAttributes(theNode),...
                    'Data', '',...
                    'Children', PC,...
                    'current',0);

if any(strcmp(methods(theNode), 'getData'))
   nodeStruct.Data = char(theNode.getData);
else
   nodeStruct.Data = '';
end
if strcmp(nodeStruct.Name,'#text')
   nodeStruct=[];
end

% ----- Subfunction PARSEATTRIBUTES -----
function attributes = parseAttributes(theNode)
% Create attributes structure.
attributes = [];
if theNode.hasAttributes
   theAttributes = theNode.getAttributes;
   numAttributes = theAttributes.getLength;
   allocCell = cell(1, numAttributes);
   attributes = struct('Name', allocCell,...
                       'Value', allocCell);

   for count = 1:numAttributes
      attrib = theAttributes.item(count-1);
      attributes(count).Name = char(attrib.getName);
      attributes(count).Value = char(attrib.getValue);
   end
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function MRSIStruct = setSize(MRSIStruct, newSize)
    MRSIStruct.sz = newSize;
end

author: CIC-methods

repo: FID-A

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


text:
function [node,elem,bdFlag,HB,tree] = bisect(node,elem,markedElem,bdFlag)
%% BISECT bisect a 2-D triangulation.
%
% [node,elem] = BISECT(node,elem,markedElem) refine the current
% triangulation by bisecting marked elements and minimal neighboring
% elements to get a conforming and shape regular triangulation. Newest
% vertex bisection is implemented. markedElem is a vector containing the
% indices of elements to be bisected. It could be a logical vector of
% length size(elem,1).
%
% [node,elem,bdFlag] = BISECT(node,elem,markedElem,bdFlag) returns the
% updated bdFlag after the bisection. It will be used for PDEs with mixed
% boundary conditions.
%
% [node,elem,bdFlag,HB,tree] = BISECT(node,elem,markedElem,bdFlag)
% returns HB and tree arrays.
%
% - HB(:,1:3) is a hierarchical basis structure for nodes, where
%   HB(:,1) is the global index of new added nodes, and HB(:,2:3) the
%   global indices of two parent nodes of new added nodes. HB is usful
%   for the interpolation between two grids; see also nodeinterpolate.
%
% - tree(:,1:3) stores the binary tree of the refinement. tree(:,1) is the
%   index of the parent element in the coarse mesh and tree(:,2:3) are two
%   children indices in the fine mesh. tree is useful for the interpolation
%   of elementwise function; see eleminterpolate.
%
% Example
%
%      node = [0,0; 1,0; 1,1; 0,1];
%      elem = [2,3,1; 4,1,3];
%      figure(1); subplot(1,3,1); showmesh(node,elem);
%      [node,elem] = bisect(node,elem,'all');
%      figure(1); subplot(1,3,2); showmesh(node,elem);
%      bdFlag = setboundary(node,elem,'Dirichlet','all','Neumann','y==1');
%      [node,elem,bdFlag] = bisect(node,elem,[1 4],bdFlag);
%      figure(1); subplot(1,3,3); showmesh(node,elem);
%
% See also bisect3, coarsen, coarsen3, nodeinterpolate, eleminterpolate.
%
% Reference page in Help browser
% ifem meshdoc
% ifem bisectdoc

% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Set up
HB = []; tree = [];
if ~exist('bdFlag','var'), bdFlag = []; end
if ~exist('markedElem','var'), markedElem = (1:size(elem,1))'; end
if isempty(markedElem), return; end
if strcmp(markedElem,'all'), markedElem = (1:size(elem,1))'; end
if islogical(markedElem), markedElem = find(markedElem); end

%% Construct auxiliary data structure
T = auxstructure(elem);
neighbor = T.neighbor; elem2edge = T.elem2edge; edge = T.edge;
clear T;
%[neighbor,elem2edge,edge] = auxstructurec(int32(elem));
N = size(node,1); NT = size(elem,1); NE = size(edge,1);

%% Add new nodes
isCutEdge = false(NE,1);
while sum(markedElem)>0
    isCutEdge(elem2edge(markedElem,1)) = true;
    refineNeighbor = neighbor(markedElem,1);
    markedElem = refineNeighbor(~isCutEdge(elem2edge(refineNeighbor,1)));
end
edge2newNode = zeros(NE,1,'uint32');
edge2newNode(isCutEdge) = N+1:N+sum(isCutEdge);
HB = zeros(sum(isCutEdge),3,'uint32');
HB(:,1) = edge2newNode(isCutEdge);
HB(:,[2 3]) = edge(isCutEdge,[1 2]);
node(HB(:,1),:) = (node(HB(:,2),:) + node(HB(:,3),:))/2;

%% Refine marked elements
Nb = 0; tree = zeros(3*NT,3,'uint32');
for k = 1:2
    t = find(edge2newNode(elem2edge(:,1))>0);
    newNT = length(t);
    if (newNT == 0), break; end
    L = t;
    R = NT+1:NT+newNT;
    p1 = elem(t,1); p2 = elem(t,2); p3 = elem(t,3);
    p4 = edge2newNode(elem2edge(t,1));
    elem(L,:) = [p4, p1, p2];
    elem(R,:) = [p4, p3, p1];
	if nargin==4 && ~isempty(bdFlag) % Refine boundary edges
   		bdFlag(R,[1 3]) = bdFlag(t,[2 1]);
   		bdFlag(L,[1 2]) = bdFlag(t,[3 1]);
        bdFlag(L,3) = 0;
    else
        bdFlag = [];
	end
    tree(Nb+1:Nb+newNT,1) = L;
    tree(Nb+1:Nb+newNT,2) = L;
    tree(Nb+1:Nb+newNT,3) = R;
    elem2edge(L,1) = elem2edge(t,3);
    elem2edge(R,1) = elem2edge(t,2);
    NT = NT + newNT;
    Nb = Nb + newNT;
end
tree = tree(1:Nb,:);

repo: ifem


text:
function [ x, y, a, b, c, ierror ] = r8poly2_ex2 ( x1, y1, x2, y2, x3, y3 )

%*****************************************************************************80
%
%% R8POLY2_EX2 finds the extremal point of a parabola determined by three points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X1, Y1, X2, Y2, X3, Y3, the coordinates of
%    three points on the parabola.  X1, X2 and X3 must be distinct.
%
%    Output, real X, Y, the X coordinate of the extremal
%    point of the parabola, and the value of the parabola at that point.
%
%    Output, real A, B, C, the coefficients that define the
%    parabola: P(X) = A * X**2 + B * X + C.
%
%    Output, integer IERROR, error flag.
%    0, no error.
%    1, two of the X values are equal.
%    2, the data lies on a straight line; there is no finite extremal
%    point.
%
  ierror = 0;
  x = x1;
  y = y1;
  a = 0.0;
  b = 0.0;
  c = 0.0;

  if ( x1 == x2 || x2 == x3 || x3 == x1 )
    ierror = 1;
    return
  end

  if ( y1 == y2 && y2 == y3 && y3 == y1 )
    x = x1
    y = y1
    return
  end
%
%  Set up the Vandermonde matrix.
%
  v(1,1) = 1.0;
  v(1,2) = x1;
  v(1,3) = x1 * x1;

  v(2,1) = 1.0;
  v(2,2) = x2;
  v(2,3) = x2 * x2;

  v(3,1) = 1.0;
  v(3,2) = x3;
  v(3,3) = x3 * x3;
%
%  Get the inverse.
%
  [ w, det ] = r8mat_inverse_3d ( v );
%
%  Compute the parabolic coefficients.
%
  c = w(1,1) * y1 + w(1,2) * y2 + w(1,3) * y3;
  b = w(2,1) * y1 + w(2,2) * y2 + w(2,3) * y3;
  a = w(3,1) * y1 + w(3,2) * y2 + w(3,3) * y3;
%
%  Determine the extremal point.
%
  if ( a == 0.0 )
    ierror = 2;
    return
  end

  x = -b / ( 2.0 * a );
  y = a * x * x + b * x + c;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%--- help for expand_optimal_policy_solution ---
%
%  INTERNAL FUNCTION
%
%
%    Other functions named expand_optimal_policy_solution
%
%       dsge/expand_optimal_policy_solution
%

author: jmaih


text:
function edge = clipLine3d(line, box)
%CLIPLINE3D Clip a line with a box and return an edge.
%
%   EDGE = clipLine3d(LINE, BOX);
%   Clips the line LINE with the bounds given in BOX, and returns the
%   corresponding edge.
%
%   If the line lies totally outside of the box, returns a 1-by-6 row array
%   containing only NaN's.
%
%   If LINE is a N-by-6 array, with one line by row, returns the clipped
%   edge coresponding to each line in a N-by-6 array.
%
%   See also:
%   lines3d, edges3d, createLine3d
%

% ---------
% author : David Legland
% e-mail: david.legland@inra.fr
% INRA - TPV URPOI - BIA IMASTE
% created the 30/10/2008 from drawLine3d

%   HISTORY
%   30/10/2008 replace intersectPlaneLine by intersectLinePlane
%   25/11/2008 improve test for bounds, and use more explicit code
%   22/06/2009 fig bug, add support for several lines
%   16/11/2010 use middle point for checking edge bounds

% get box limits
xmin = box(1); xmax = box(2);
ymin = box(3); ymax = box(4);
zmin = box(5); zmax = box(6);

% extreme corners of the box
p000 = [xmin ymin zmin];
p111 = [xmax ymax zmax];

% main vectors
ex   = [1 0 0];
ey   = [0 1 0];
ez   = [0 0 1];

% box faces parallel to Oxy
planeZ0 = [p000 ex ey];
planeZ1 = [p111 ex ey];

% box faces parallel to Oxz
planeY0 = [p000 ex ez];
planeY1 = [p111 ex ez];

% box faces parallel to Oyz
planeX0 = [p000 ey ez];
planeX1 = [p111 ey ez];

% number of lines
nLines = size(line, 1);

% allocate memory for result
edge = zeros(nLines, 6);

% iterate over lines to clip
for i = 1:nLines

    % compute intersection point with each plane
    ipZ0 = intersectLinePlane(line(i,:), planeZ0);
    ipZ1 = intersectLinePlane(line(i,:), planeZ1);
    ipY0 = intersectLinePlane(line(i,:), planeY0);
    ipY1 = intersectLinePlane(line(i,:), planeY1);
    ipX1 = intersectLinePlane(line(i,:), planeX1);
    ipX0 = intersectLinePlane(line(i,:), planeX0);

    % concatenate resulting points
    points  = [ipX0;ipX1;ipY0;ipY1;ipZ0;ipZ1];

    % compute position of each point on the line
    pos     = linePosition3d(points, line(i,:));

    % keep only defined points
    ind     = find(~isnan(pos));
    pos     = pos(ind);
    points  = points(ind,:);

    % sort points with respect to their position
    [pos, ind] = sort(pos); %#ok<ASGLU>
    points  = points(ind, :);

    % keep median points wrt to position. These points define the limit of
    % the clipped edge.
    nv      = length(ind)/2;

    % create resulting edge.
    edge(i,:)   = [points(nv, :) points(nv+1, :)];
end

% check that middle point of the edge is contained in the box
midX = mean(edge(:, [1 4]), 2);
xOk  = xmin <= midX & midX <= xmax;
midY = mean(edge(:, [2 5]), 2);
yOk  = ymin <= midY & midY <= ymax;
midZ = mean(edge(:, [3 6]), 2);
zOk  = zmin <= midZ & midZ <= zmax;

% if one of the bounding condition is not met, set edge to NaN
edge (~(xOk & yOk & zOk), :) = NaN;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function  z = embedding(x,types,positions,w,dropout)
% embedding   The BERT embeddings of encoded tokens, token types and token
% positions.
%
%   Z = embedding(X,types,positions,weights,dropoutProbability) computes
%   the embedding of encoded tokens X, token types specified by types, and
%   token positions. Inputs X, types and positions are
%   1-by-numInputTokens-by-numObs unformatted dlarray-s. The types take
%   values 1 or 2. The weights input is a struct of embedding weights such
%   as mdl.Parameters.Weights.embeddings where mdl = bert(). The
%   dropoutProbability is a scalar double between 0 and 1 corresponding to
%   the post-embedding dropout probability.

% Copyright 2021 The MathWorks, Inc.
wordEmbedding = embed(x,w.word_embeddings,'DataFormat','CTB');
typeEmbedding = embed(types,w.token_type_embeddings,'DataFormat','CTB');
positionEmbedding = embed(positions,w.position_embeddings,'DataFormat','CTB');
z = wordEmbedding+typeEmbedding+positionEmbedding;
z = transformer.layer.normalization(z,w.LayerNorm.gamma,w.LayerNorm.beta);
z = transformer.layer.dropout(z,dropout);
end

author: matlab-deep-learning

repo: transformer-models


text:
function GenerateSNIRFReaderWriter()

dirnameApp = [getAppDir(), 'DataTree/AcquiredData/'];
dirnameUtils = [getAppDir(), 'Utils/'];
dirnameInstall = fileparts(which('GenerateSNIRFReaderWriter.m'));
dirnameInstall(dirnameInstall=='\') = '/';
dirnameInstall(end+1)='/';

if exist([dirnameInstall, '/snirf_homer3'], 'dir')
    rmdir([dirnameInstall, '/snirf_homer3'], 's');
end
mkdir([dirnameInstall, '/snirf_homer3']);
mkdir([dirnameInstall, '/snirf_homer3/Utils']);

copyFiles([dirnameApp, 'DataFiles'],           [dirnameInstall, 'snirf_homer3/DataFiles']);
copyFiles([dirnameApp, 'Nirs'],                [dirnameInstall, 'snirf_homer3/Nirs']);
copyFiles([dirnameApp, 'Snirf'],               [dirnameInstall, 'snirf_homer3/Snirf']);
copyFiles([dirnameUtils, 'pretty_print_struct.m'], [dirnameInstall, 'snirf_homer3']);
copyFiles([dirnameUtils, 'isproperty.m'],      [dirnameInstall, 'snirf_homer3/Utils']);
copyFiles([dirnameUtils, 'iswholenum.m'],      [dirnameInstall, 'snirf_homer3/Utils']);
copyFiles([dirnameUtils, 'propnames.m'],       [dirnameInstall, 'snirf_homer3/Utils']);
copyFiles([dirnameUtils, 'CopyHandles.m'],     [dirnameInstall, 'snirf_homer3/Utils']);
copyFiles([dirnameUtils, 'strtrim_improve.m'], [dirnameInstall, 'snirf_homer3/Utils']);
copyFiles([dirnameUtils, 'fullpath.m'],        [dirnameInstall, 'snirf_homer3']);
copyFiles([dirnameUtils, 'filesepStandard.m'], [dirnameInstall, 'snirf_homer3/Utils']);
copyFiles([dirnameUtils, 'cell2str_new.m'], [dirnameInstall, 'snirf_homer3/Utils']);
copyFiles([dirnameApp, 'AcqDataClass.m'],      [dirnameInstall, 'snirf_homer3']);
copyFiles([dirnameApp, 'README.md'],           [dirnameInstall, 'snirf_homer3']);
copyFiles([dirnameUtils, 'cd_safe.m'],          [dirnameInstall, 'snirf_homer3']);
copyFiles([getAppDir(), 'str2cell.m'],          [dirnameInstall, 'snirf_homer3/Utils']);

% -------------------------------------------------------------------
function copyFiles(src, dst, type)

if ~exist('type', 'var')
    type = 'file';
end
if ~exist('errtype', 'var')
    errtype = 'Error';
end

try
    % If src is one of several possible filenames, then src to any one of
    % the existing files.
    if iscell(src)
        for ii=1:length(src)
            if ~isempty(dir(src{ii}))
                src = src{ii};
                break;
            end
        end
    end

    assert(logical(exist(src, type)));

    % Check if we need to untar the file
    k = findstr(src,'.tar.gz');
    if ~isempty(k)
        untar(src,fileparts(src));
        src = src(1:k-1);
    end

    % Copy file from source to destination folder
    fprintf('Copying %s to %s\n', src, dst);
    copyfile(src, dst);

catch ME

    printStack();
    if iscell(src)
        src = src{1};
    end
    MenuBox(sprintf('Error: Could not copy %s to installation folder.', src), {'OK'});
    pause(5);
    rethrow(ME);

end

author: BUNPC


text:
function prob_test126 ( )

%*****************************************************************************80
%
%% TEST126 tests POISSON_CDF, POISSON_CDF_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST126:\n' );
  fprintf ( 1, '  POISSON_CDF evaluates the cumulative distribution\n' );
  fprintf ( 1, '    function for the discrete Poisson probability\n' );
  fprintf ( 1, '    density function.\n' );
  fprintf ( 1, '  POISSON_CDF_VALUES returns some exact values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A is the expected mean number of successes per unit time;\n' );
  fprintf ( 1, '  X is the number of successes;\n' );
  fprintf ( 1, '  POISSON_CDF is the probability of having up to X\n' );
  fprintf ( 1, '  successes in unit time.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   A          X   Exact F     POISSON_CDF(A,X)\n' );
  fprintf ( 1, '\n' );
  n_data = 0;

  while ( 1 )

    [ n_data, a, x, fx ] = poisson_cdf_values ( n_data );

    if ( n_data == 0 );
      break
    end

    fx2 = poisson_cdf ( x, a );

    fprintf ( 1, '  %8f  %8d  %14f  %14f\n', a, x, fx, fx2 );

  end

  return
end


text:
function plotIntCon(prob,data)
%PLOTINTCON Plot Integer Constraints on the current figure
%   plotIntCon(prob)

%   Copyright (C) 2011 Jonathan Currie (IPL)

%Check we have integer variables in this plot
if(all(prob.int.str(data.idx) == 'C'))
    return; %nothing to do!
end

%Colours
intC = [0 114/255 189/255]; %feasible
intIC = [123/255 142/255 147/255]; %infeasible

xl = round(data.xl); yl = round(data.yl);

hold on;

%Get Problem Components
if(~isempty(prob.rl))
    [A,b,Aeq,beq] = row2gen(prob.A,prob.rl,prob.ru);
else
    A = prob.A; b = prob.b;
    Aeq = prob.Aeq; beq = prob.beq;
end
Q = prob.Q; l = prob.l; qrl = prob.qrl; qru = prob.qru;
if(isempty(prob.lb)), lb = -Inf; else lb = prob.lb; end
if(isempty(prob.ub)), ub = Inf; else ub = prob.ub; end

%Check for 1D problem
if(prob.sizes.ndec == 1)
    if(prob.int.ind(1) == -1)
        X = 0:1; Y = [];
    else
        X = (floor(xl(1))+1):(ceil(xl(2))-1); Y = [];
    end
else
    %Generate Grid and check for binary constraints
    if(prob.int.ind(1) == -1)
        x = 0:1;
    else
        try
            x = (floor(xl(1))+1):(ceil(xl(2))-1);
        catch ME
            optiwarn('opti:plotint',['Could not plot integer constraints due to plotting error: ' ME.message]);
            return;
        end
    end
    if(prob.int.ind(2) == -1)
        y = 0:1;
    else
        try
            y = (floor(yl(1))+1):(ceil(yl(2))-1);
        catch ME
            optiwarn('opti:plotint',['Could not plot integer constraints due to plotting error: ' ME.message]);
            return;
        end
    end
    [X,Y] = meshgrid(x,y);
end

%Check for row nonlinear constraints
if(~isempty(prob.cl))
    prob = nrow2mix(prob,0,false);
end

%Work out feasible integer points
idx = zeros(numel(X),1); Xv = X(:); Yv = Y(:);
for i=1:numel(X)
    if(isempty(Yv))
        xy = Xv(i);
    else
        xy = data.fixval;
        xy(data.idx) = [Xv(i), Yv(i)]';
    end
    if(~isempty(prob.A))
        idx(i) = all(b >= A*xy) && all(xy >= lb) && all(xy <= ub);
    else
        idx(i) = all(xy >= lb) && all(xy <= ub);
    end
    if(~isempty(Aeq))
        idx(i) = idx(i) && (abs(Aeq*xy-beq) < 1e-6);
    end
    if(~isempty(Q))
        for n = 1:prob.sizes.nqc
            if(iscell(Q))
                nQ = Q{n}; nl = l(:,n); nrl = qrl(n); nru = qru(n);
            else
                nQ = Q; nl = l; nrl = qrl; nru = qru;
            end
            if(~isinf(nru))
                qu = all(xy'*nQ*xy + nl'*xy <= nru);
            else
                qu = true;
            end
            if(~isinf(nrl))
                ql = all(xy'*-nQ*xy - nl'*xy <= -nrl);
            else
                ql = true;
            end
            idx(i) = idx(i) && ql && qu;
        end
    end
    if(~isempty(prob.nle))
        vals = prob.nlcon(xy);
        leI = find(prob.nle == -1);
        geI = find(prob.nle == 1);
        eqI = find(prob.nle == 0);
        nlrhs = prob.nlrhs;
        idx(i) = idx(i) && all(vals(leI) <= nlrhs(leI)) && all(vals(geI) >= nlrhs(geI)) && all(abs(vals(eqI)-nlrhs(eqI)) < 1e-6);
    end
end
%Plot points
idx = logical(idx);
if(isempty(Y))
    Y = zeros(size(X));
    for i = 1:length(Y)
        Y(i) = prob.objective(X(i));
    end
end
plot(X(idx),Y(idx),'o','color',intC,'markerface',intC,'markersize',5);
plot(X(~idx),Y(~idx),'o','color',intIC,'markersize',5);

hold off;

end


text:
%	OVERVIEW:
%       This demonstration analyzes a segment of 5-minutes 'raw' data
%       with known atrial fibrillation to show the operation of the
%       AF detection algorithm.
%   OUTPUT:
%       HRV Metrics exported to .cvs files
%
%   DEPENDENCIES & LIBRARIES:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   REFERENCE:
%       Vest et al. "An Open Source Benchmarked HRV Toolbox for Cardiovascular
%       Waveform and Interval Analysis" Physiological Measurement (In Press), 2018.
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Giulia Da Poian
%	COPYRIGHT (C) 2018
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc; close all;

run(['..' filesep 'startup.m'])

% Remove old files generated by this demo
OldFolder = [pwd,filesep, 'OutputData', filesep, 'ResultsAFData'];
if exist(OldFolder, 'dir')
    rmdir(OldFolder, 's');
    fprintf('Old Demo Folder deleted \n');
end

HRVparams = InitializeHRVparams('demoAF'); % include the project name
HRVparams.poincare.on = 0; % Poincare analysis off for this demo
HRVparams.DFA.on = 0; % DFA analysis off for this demo
HRVparams.MSE.on = 0; % MSE analysis off for this demo
HRVparams.HRT.on = 0; % HRT analysis off for this demo

[subjectIDs,filesTBA] = GenerateListOfFilesTBA(HRVparams.ext,...
                                                    HRVparams.readdata,0);
idx = find(strcmp(subjectIDs,'TestAFdata'));
i_patient = idx;

% 1. Load Raw Patient Data

load(filesTBA{i_patient});
% 2. Analyze data using HRV PhysioNet Cardiovascular Signal Toolbox
[results, resFilename] = Main_HRV_Analysis(signal(:,1),[],'ECGWaveform',...
                                          HRVparams,subjectIDs(i_patient));

% 3. Compare generated output file with the reference one

currentFile = strcat(HRVparams.writedata, filesep, resFilename.HRV, '.csv');
referenceFile = strcat('ReferenceOutput', filesep, 'AFDemo.csv');
testHRV = CompareOutput(currentFile,referenceFile);

% 3. Load QRS annotation saved by Main_HRV_Analysis
annotName = strcat(HRVparams.writedata, filesep, 'Annotation',filesep,...
                                                    subjectIDs(i_patient));
jqrs_ann = read_ann( annotName{1} , 'jqrs');
wqrs_ann = read_ann( annotName{1} , 'wqrs');

% For demo pourpose recompute bsqi
[sqijw, StartIdxSQIwindows] = bsqi(jqrs_ann,wqrs_ann,HRVparams);

HRVparams.gen_figs = 1;
% Plot detected beats
if HRVparams.gen_figs
    Plot_SignalDetection_SQI(time, signal(:,1), jqrs_ann, sqijw,'ECG')
end

if testHRV
    fprintf('** DemoRawDataAF: TEST SUCCEEDED ** \n ')
    fprintf('A file named %s.csv \n has been saved in %s \n', ...
    resFilename.HRV, HRVparams.writedata);
else
    fprintf('** DemoRawDataAF: TEST FAILED ** \n')
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function sVF = grad(sF, varargin)
% gradient of a spherical function
%
% Syntax
%   sVF = grad(sF)  % returns the gradient as a spherical vector field
%   g = grad(sF, v) % return the gradient in point v as vector3d
%
% Input
%  sF - @S2FunHarmonic
%  v - @vector3d
%
% Output
%  sVF - @sphericalVectorFieldHarmonic
%    g - @vector3d
%

if nargin > 1

  delta = 0.1*degree;

  N = varargin{1};   % normal direction
  t1 = perp(N);      % first tangential vector
  t2 = normalize(cross(N,t1));  % second tangential vector

  % the exponential map on the sphere
  v2 = rotation.byAxisAngle(t1,-delta) .* N;
  v1 = rotation.byAxisAngle(t2,delta) .* N;

  sFN = sF.eval(N);

  sVF = ((sF.eval(v1)-sFN) .* t1 + (sF.eval(v2)-sFN) .* t2) ./ delta;

else
  sVF = S2VectorFieldHarmonic.quadrature(@(v) grad(sF,v));
end

end

author: mtex-toolbox

repo: mtex

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