text:
function [HBLog2D, HBLogAv, ...
    mB2D, mBAv, ...
    names2D, namesAv, ...
    namesmB2D, namesmBAv] = ...
    extractHillBands(data, params)

nChans = size(data,2);
nBands = numel(params.HillsBands.Range);

HBLog = NaN(nBands, nChans);
for c = 1:nChans
    % T = 1/Fs;
    L = size(data,1);
    % t = (0:L-1)*T;

    Y = fft(data(:,c));

    P2 = abs(Y/L);
    P1 = P2(1:L/2+1);
    P1(2:end-1) = 2*P1(2:end-1);

    f = params.Fs*(0:(L/2))/L;

    if params.plotOn
        plot(f,P1)
        title('Single-Sided Amplitude Spectrum of X(t)')
        xlabel('f (Hz)')
        ylabel('|P1(f)|')
    end

    for h = 1:nBands
        hz = params.HillsBands.Range(h);

        HBLog(h,c) = mean(P2(round(f) == hz));

    end
end

% Log10
HBLog2D = log10(HBLog);
% Remove infs
HBLog2D(isinf(HBLog2D)) = 0;
% Mean across channels
HBLogAv = mean(HBLog2D,2)';

[~, mIdx] = max(HBLog2D);
mB2D = single(mIdx);
mBAv = mean(mB2D);

namesmB2D = (string('maxHills_c') + (1:16)')';
namesmBAv = 'maxHillsAv';

nFs = 47;
HBLog2D = reshape(HBLog2D, 1, nChans*nFs);

names2D = cellstr([repmat('hillBands2D_', nFs*16, 1), ...
    repmat(num2str((1:nFs)'), 16,1), ...
    repmat('_c', nFs*16,1), ...
    num2str(reshape(repmat((1:16),nFs,1),nFs*16,1))])';
names2D = strrep(names2D, ' ', '');

namesAv = cellstr([repmat('hillsBandAv_',47,1), num2str((1:47)')])';
namesAv = strrep(namesAv, ' ', '');

author: garethjns

repo: Kaggle-EEG


text:
classdef TestConvertPointsFromHomogeneous
    %TestConvertPointsFromHomogeneous

    methods (Static)
        function test_numeric_3d
            pts = shiftdim([1,2,1; 4,5,1], -1);
            dst = cv.convertPointsFromHomogeneous(pts);
            validateattributes(dst, {class(pts)}, {'size',[2 1 3-1]});
        end

        function test_numeric_4d
            pts = shiftdim([1,2,3,1; 4,5,6,1], -1);
            dst = cv.convertPointsFromHomogeneous(pts);
            validateattributes(dst, {class(pts)}, {'size',[2 1 4-1]});
        end

        function test_cellarray_3d
            pts = {[1,2,1], [4,5,1]};
            dst = cv.convertPointsFromHomogeneous(pts);
            validateattributes(dst, {'cell'}, {'vector', 'numel',numel(pts)});
            cellfun(@(p) validateattributes(p, {'numeric'}, ...
                {'vector', 'numel',3-1}), dst);
        end

        function test_cellarray_4d
            pts = {[1,2,3,1], [4,5,6,1]};
            dst = cv.convertPointsFromHomogeneous(pts);
            validateattributes(dst, {'cell'}, {'vector', 'numel',numel(pts)});
            cellfun(@(p) validateattributes(p, {'numeric'}, ...
                {'vector', 'numel',4-1}), dst);
        end

        function test_all_formats_and_dims
            verify_numeric = @(M,n,d) assert( ...
                ismatrix(M) && isequal(size(M),[n d-1]));
            verify_cell = @(C,n,d) assert( ...
                iscell(C) && numel(C)==n && ...
                all(cellfun(@isvector, C)) && ...
                all(cellfun(@numel, C) == (d-1)));
            % 3D/4D to 2D/3D
            n = 10;
            klass = {'double', 'single'};
            for k=1:numel(klass)
                for d=[3 4]
                    % Nxd numeric matrix
                    ptsH = rand([n,d], klass{k});
                    pts = cv.convertPointsFromHomogeneous(ptsH);
                    validateattributes(pts, {klass{k}}, {'2d', 'size',[n d-1]});

                    % Nx1xd numeric matrix
                    ptsH = rand([n,1,d], klass{k});
                    pts = cv.convertPointsFromHomogeneous(ptsH);
                    validateattributes(pts, {klass{k}}, {'3d', 'size',[n 1 d-1]});

                    % 1xNxd numeric matrix
                    ptsH = rand([1,n,d], klass{k});
                    pts = cv.convertPointsFromHomogeneous(ptsH);
                    validateattributes(pts, {klass{k}}, {'3d', 'size',[n 1 d-1]});

                    % Nx1 cell array of 1xd points
                    ptsH = num2cell(rand([n,d], klass{k}), 2);
                    pts = cv.convertPointsFromHomogeneous(ptsH);
                    validateattributes(pts, {'cell'}, {'vector', 'numel',n});
                    cellfun(@(p) validateattributes(p, {'numeric'}, ...
                        {'vector', 'numel',d-1}), pts);

                    % 1xN cell array of 1xd points
                    ptsH = num2cell(rand([n,d], klass{k}), 2).';
                    pts = cv.convertPointsFromHomogeneous(ptsH);
                    validateattributes(pts, {'cell'}, {'vector', 'numel',n});
                    cellfun(@(p) validateattributes(p, {'numeric'}, ...
                        {'vector', 'numel',d-1}), pts);

                    % 1xN cell array of dx1 points
                    ptsH = num2cell(rand([d,n], klass{k}), 1);
                    pts = cv.convertPointsFromHomogeneous(ptsH);
                    validateattributes(pts, {'cell'}, {'vector', 'numel',n});
                    cellfun(@(p) validateattributes(p, {'numeric'}, ...
                        {'vector', 'numel',d-1}), pts);

                    % Nx1 cell array of dx1 points
                    ptsH = num2cell(rand([d,n], klass{k}), 1).';
                    pts = cv.convertPointsFromHomogeneous(ptsH);
                    validateattributes(pts, {'cell'}, {'vector', 'numel',n});
                    cellfun(@(p) validateattributes(p, {'numeric'}, ...
                        {'vector', 'numel',d-1}), pts);
                end
            end
        end

        function test_last_dimension
            % xn = 0, xn = 1, and arbitrary xn
            for xn=[0 1 rand()]
                for d=[3 4]
                    ptsH = {[rand(1,d-1),xn], [rand(1,d-1),xn]};
                    pts = cv.convertPointsFromHomogeneous(ptsH);
                end
            end
        end

        function test_error_argnum
            try
                cv.convertPointsFromHomogeneous();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%########################################################################
%
%   - PPGI Toolbox -
%   A MATLAB toolbox for Photoplethysmography Imaging (PPGI)
%
% Author   : Christian S. Pilz
% Company  : The Nature of Space of Time
% Date     : 07.05.2019
%
% Contact  : cpi@partofthestars.com
% Web Page : www.partofthestars.com
%
% Version  : beta0.1
%
%########################################################################
%
%	circular_buffer.m:
%
% Description:
%
%   A  MATLAB circular buffer class to handle streaming data
%

classdef circular_buffer

   properties
        buf_;
        head_;
        tail_;
        max_size_;
        full_;
        cell_
   end

   methods

       function obj = circular_buffer(max_size,cell)

          if nargin > 0
             if isnumeric(max_size)
                obj.max_size_ = max_size;
                obj.buf_= [];
                obj.head_= 1;
                obj.tail_= 1;
                obj.full_= logical(false);
                obj.cell_= logical(cell);
             else
                error('size must be numeric')
             end
          end

       end

       function ret = capacity(obj)
           ret = obj.max_size_;
       end

       function ret =size(obj)
           ret = obj.max_size_;

            if obj.full_== 0
                if obj.head_ >= obj.tail_
                    ret = obj.head_ - obj.tail_;
                else
                    %ret = obj.max_size_ + obj.head_ - obj.tail_;
                    ret = obj.head_ - obj.tail_;
                end
            end
       end

       function obj = put(obj,data)
           if obj.cell_
               obj.buf_{obj.head_} = data;
           else
               obj.buf_(obj.head_,:) = data;
           end

%            if obj.full_
%                obj.tail_ = mod(obj.head_, obj.max_size_)+1;
%            end

           obj.head_ = mod(obj.head_, obj.max_size_)+1;
           obj.full_ = (obj.head_ == obj.tail_);
       end

       function [values, obj] = get(obj,window_size,shift)
           if obj.empty()
               %values=[];
               return;
           end
           try
               if obj.cell_
                   values=obj.buf_(1,obj.tail_:obj.tail_+window_size-1)
               else
                   if (obj.tail_+window_size-1)<=obj.max_size_
                       values=obj.buf_(obj.tail_:obj.tail_+window_size-1,:);
                   else
                       tail_delta=obj.max_size_-obj.tail_;
                       values=obj.buf_(obj.tail_:obj.tail_+tail_delta,:);
                       head_delta=window_size-tail_delta-1;
                       values=[values; obj.buf_(1:head_delta,:)];
                   end
               end
           catch exception
                throw(exception)
           end

           obj.tail_ = mod(obj.tail_ + shift, obj.max_size_);
           if obj.tail_ == 0
               obj.tail_=1;
           end
       end

       function obj = reset(obj)
           obj.head_ = obj.tail_;
           obj.full_ = logical(false);
       end

       function ret = empty(obj)
           %if head and tail are equal, we are empty
           ret = (obj.full_ == 0) & (obj.head_ == obj.tail_);
       end

       function ret = full(obj)
           %if tail is ahead the head by 1, we are full
           ret = obj.full_;
       end
   end
end

author: partofthestars

repo: PPGI-Toolbox


text:
%% computeTRNrisk.m
% This script does validation of vectors W acquired from training data
% (see test_bmrm.m)
%
% 27-04-11 Michal Uricar

clc;
close all; clearvars;

%% Timestamp

fprintf(1,'Started on %s\n\n', datestr(now));

%% Add path

addpath('./Functions/');
addpath('../../matlab_toolbox/mex/');

%% Select functions

load('./MAT/functions_gdisp.mat');

%%
% load config
load('./MAT/config.mat');
lbpconfig_fname = ['./MAT/exp03_lbpconfig_nImages-' config.nImages '_hop-' num2str(config.heightOfPyramid) '.mat'];
images_fname = ['./MAT/Images_nImages-' config.nImages '.mat'];
gt_fname = ['./MAT/GT_nImages-' config.nImages '.mat'];

fprintf('Loading VAL data...\n');
load(lbpconfig_fname);
load(images_fname);
load(gt_fname);
load('./results/exp03_options.mat');
S = options.S;

switch (config.nImages)
    case 'all'
        nImages = size(Images, 2);
    case 'half'
        nImages = size(Images, 2)/2;
    otherwise
        nImages = str2num(config.nImages);
end;

tic
data.lbp = lbpconfig; data.imSize = imSize; data.Images = Images; data.nImages = nImages;
options.S = S; data.Y = gt; data.options = options;
data.tmpData = getPsiMat(data, 1);
if (displacement)
    data.options.PsiG = getDisplacements(data.options);
end;
clear lbpconfig imSize Images nImages gt gt_fname images_fname lbpconfig_fname options S

% create mapTable for extracting components of W
data.mapTable = f_createMapTable(data.options, data.tmpData);
toc
fprintf('Done.\n');

%%

lambdas = [1e-2, 1e-1, 1];

% start matlabpool
try
    matlabpool;
catch ME1
    idSegLast = regexp(ME1.identifier, '(?<=:)\w+$', 'match');
    if (strcmp(idSegLast, 'OpenConnection'))
        fprintf('Already connected to matlab pool. Skipping...\n');
    end;
end

fprintf('Running validation...\n');

Rtrn = zeros(1, length(lambdas));
% for lambda = lambdas
for j = 1 : length(lambdas)
    lambda = lambdas(j);
    load(['./results/exp03_W_nImages-' config.nImages '_hop-' num2str(config.heightOfPyramid) '_lambda_' num2str(lambda) '.mat']);
    L = zeros(1, data.nImages);
    fprintf('TRN: Calculating Rtrn for lambda %f.\n', lambda);
%     for i = 1 : data.nImages
    parfor par_i = 1 : data.nImages
        GT = data.Y{par_i};
        Test = getPsiMat(data, par_i);
        [ q, g ] = getQG(data.options, W, Test, data.mapTable);
        P = f_argmax(data.options, q, g);
        L(par_i) = 1/data.options.M * sum( sqrt( sum( (GT - P ).^2 ) ) );
    end;
    Rtrn(j) = 1/data.nImages * sum(L);
end;

try
    matlabpool close;
catch ME2
    idSegLast = regexp(ME2.identifier, '(?<=:)\w+$', 'match');
    if (strcmp(idSegLast, 'OpenConnection'))
        fprintf('Cannot close matlabpool. Skipping...\n');
    end;
end;

fprintf('TRN: lambda = %s\n', num2str(lambdas, '%4.5f '));
fprintf('TRN: R      = %s\n', num2str(Rtrn, ' %4.5f '));

fprintf('Saving lambdas and Rtrns... ');
save(['./results/exp03_Rtrns_nImages-' config.nImages '_hop-' num2str(config.heightOfPyramid) '.mat'], 'lambdas', 'Rtrn');
fprintf('done.\n');

%% Timestamp

fprintf(1,'Finished on %s\n\n', datestr(now));

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function [K, sK] = indexKernCompute(kern, x, x2)

% INDEXKERNCOMPUTE Compute the INDEX kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the index based covariance function
% kernel given inputs associated with rows and columns.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 : the input matrix associated with the columns of the kernel.
% RETURN k : the kernel matrix computed at the given points.
%
% FORMAT
% DESC computes the kernel matrix for the index based covariance function
% kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : the kernel matrix computed at the given points.
%
% SEEALSO : indexKernParamInit, kernCompute, kernCreate, indexKernDiagCompute, indexardKernCompute
%
% COPYRIGHT : Neil D. Lawrence, 2011

% KERN
  if size(x, 2)>1
    error('Index kernel requires 1-dimensional input.')
  end

  if nargin<3
    x2 = x;
  end
  sK = zeros(size(x, 1), size(x2, 1));
  for i = 1:size(x, 1)
    sK(i, :) = round(x(i))==round(x2)';
  end
  sK = sparse(sK);
  K = sK*kern.variance;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function varargout = drawEdge(varargin)
%DRAWEDGE Draw an edge given by 2 points
%
%   drawEdge(x1, y1, x2, y2);
%   draw an edge between the points (x1 y1) and  (x2 y2).
%
%   drawEdge([x1 y1 x2 y2]) ;
%   drawEdge([x1 y1], [x2 y2]);
%   specify data either as bundled edge, or as 2 points
%
%   The function supports 3D edges:
%   drawEdge([x1 y1 z1 x2 y2 z2]);
%   drawEdge([x1 y1 z1], [x2 y2 z2]);
%   drawEdge(x1, y1, z1, x2, y2, z2);
%
%   Arguments can be single values or array of size [N*1]. In this case,
%   the function draws multiple edges.
%
%   H = drawEdge(..., OPT), with OPT being a set of pairwise options, can
%   specify color, line width and so on...
%
%   H = drawEdge(...) return handle(s) to created edges(s)
%
%   See also:
%   edges2d, drawCenteredEdge, drawLine
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 31/10/2003.
%

%   HISTORY
%   19/02/2004 add support for arrays of edges.
%   31/03/2004 change format of edges to [P1 P2] and variants.
%   28/11/2004 add support for 3D edges
%   01/08/2005 add support for drawing options
%   31/05/2007 update doc, and code makeup
%   03/08/2010 re-organize code

% separate edge and optional arguments
[ax, edge, options] = parseInputArguments(varargin{:});

% draw the edges
if size(edge, 2) == 4
    h = drawEdge_2d(ax, edge, options);
else
    h = drawEdge_3d(ax, edge, options);
end

% eventually return handle to created edges
if nargout > 0
    varargout = {h};
end

function h = drawEdge_2d(ax, edge, options)

h = -1 * ones(size(edge, 1), 1);

for i = 1:size(edge, 1)
    if isnan(edge(i,1))
        continue;
    end

    h(i) = plot(ax, edge(i, [1 3]), edge(i, [2 4]), options{:});
end

function h = drawEdge_3d(ax, edge, options)

h = -1 * ones(size(edge, 1), 1);

for i = 1:size(edge, 1)
    if isnan(edge(i,1))
        continue;
    end

    h(i) = plot3(ax, edge(i, [1 4]), edge(i, [2 5]), edge(i, [3 6]), options{:});
end

function [ax, edge, options] = parseInputArguments(varargin)

% extract handle of axis to draw on
if isAxisHandle(varargin{1})
    ax = varargin{1};
    varargin(1) = [];
else
    ax = gca;
end

% find the number of arguments defining edges
nbVal = 0;
for i = 1:length(varargin)
    if isnumeric(varargin{i})
        nbVal = nbVal+1;
    else
        % stop at the first non-numeric value
        break;
    end
end

% extract drawing options
options = varargin(nbVal+1:end);

% ensure drawing options have correct format
if length(options) == 1
    options = [{'color'}, options];
end

% extract edges characteristics
switch nbVal
    case 1
        % all parameters in a single array
        edge = varargin{1};

    case 2
        % parameters are two points, or two arrays of points, of size N*2.
        p1 = varargin{1};
        p2 = varargin{2};
        edge = [p1 p2];

    case 4
        % parameters are 4 parameters of the edge : x1 y1 x2 and y2
        edge = [varargin{1} varargin{2} varargin{3} varargin{4}];

    case 6
        % parameters are 6 parameters of the edge : x1 y1 z1 x2 y2 and z2
        edge = [varargin{1} varargin{2} varargin{3} varargin{4} varargin{5} varargin{6}];

    otherwise
        error('drawEdge:WrongNumberOfParameters', 'Wrong number of parameters');
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [ x, seed ] = anglit_sample ( seed )

%*****************************************************************************80
%
%% ANGLIT_SAMPLE samples the Anglit PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ cdf, seed ] = r8_uniform_01 ( seed );

  x = anglit_cdf_inv ( cdf );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [retX,retY] = get_xy(dat,ind)

%   get_xy(data)           returns X and Y matrices of the data_global object
%   get_xy(data,indexes)   returns X and Y matrices of the data_global object for given
%                           indexes

global X;
global Y;

if nargin==1
  if isempty(dat.myX)
      ind=[1:length(dat.index)];
  else
      ind=[1:size(dat.myX,1)];
  end;
end;
if ~isempty(dat.myX)
  retX=dat.myX(ind,:);
else
  retX=X(dat.index(ind),dat.findex);
end
if ~isempty(dat.myY)
  retY=dat.myY(ind,:);
else
  retY=Y(dat.index(ind),:);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function tap_delay = tapped_delay_fcn(input)
% The Tapped Delay function delays its input by the specified number
% of sample periods, and outputs all the delayed versions in a vector
% form. The output includes current input

% NOTE: To instruct Embedded MATLAB to compile an external function,
% add the following compilation directive or pragma to the function code
%#eml

persistent u_d;
if isempty(u_d)
    u_d = fi(zeros(1,40), numerictype(input), fimath(input));
end

u_d = [u_d(2:40), input];

tap_delay = u_d;

author: Sable

repo: mcbench-benchmarks


text:
function aedat = TrimTimeStampReset(aedat)

%{
If timestamp reset events occur, timestamps become non-monotonic.
Resolve this by searching for non-monotonicites in the timestamps for each
data type and trimming away everything before the last non-monotonicity.
A better technique would look at the packetTimestamps, at least in aedat3
imports.
%}

dbstop if error

% Special

if ~isfield(aedat, 'data')
    disp('No data found')
    return
end
if isfield(aedat.data, 'special')
    ts = int64(aedat.data.special.timeStamp);
    lastTsReset = find(ts(2 : end) - ts(1 : end - 1) < 0, 1, 'last');
    if ~isempty(lastTsReset)
        keepLogical = true(aedat.data.special.numEvents, 1);
        keepLogical(1 : lastTsReset) = false;
        aedat.data.special.timeStamp = aedat.data.special.timeStamp (keepLogical);
        aedat.data.special.address       = aedat.data.special.address  (keepLogical);
        aedat.data.special.numEvents = nnz(keepLogical);
        if isfield(aedat.data.special, 'valid')
            aedat.data.special.valid    = aedat.data.special.valid (keepLogical);
        end
    end
end
% Polarity

if isfield(aedat.data, 'polarity')
    ts = int64(aedat.data.polarity.timeStamp);
    lastTsReset = find(ts(2 : end) - ts(1 : end - 1) < 0, 1, 'last');
    if ~isempty(lastTsReset)
        keepLogical = true(aedat.data.polarity.numEvents, 1);
        keepLogical(1 : lastTsReset) = false;
        aedat.data.polarity.timeStamp   = aedat.data.polarity.timeStamp (keepLogical);
        aedat.data.polarity.x           = aedat.data.polarity.x         (keepLogical);
        aedat.data.polarity.y           = aedat.data.polarity.y         (keepLogical);
        aedat.data.polarity.polarity    = aedat.data.polarity.polarity  (keepLogical);
        aedat.data.polarity.numEvents   = nnz(keepLogical);
        if isfield(aedat.data.polarity, 'valid')
            aedat.data.polarity.valid    = aedat.data.polarity.valid (keepLogical);
        end
    end
end

% Frames

% This assumes that timestamps have been simplified to aedat2 standard, if
% they came from aedat3 file
if isfield(aedat.data, 'frame')
    ts = int64(aedat.data.frame.timeStampStart);
    lastTsReset = find(ts(2 : end) - ts(1 : end - 1) < 0, 1, 'last');
    if ~isempty(lastTsReset)
        keepLogical = true(aedat.data.frame.numEvents, 1);
        keepLogical(1 : lastTsReset) = false;
        aedat.data.frame.timeStampStart = aedat.data.frame.timeStampStart (keepLogical);
        aedat.data.frame.timeStampEnd = aedat.data.frame.timeStampEnd (keepLogical);
        aedat.data.frame.samples           = aedat.data.frame.samples         (keepLogical);
        aedat.data.frame.xLength           = aedat.data.frame.xLength         (keepLogical);
        aedat.data.frame.yLength           = aedat.data.frame.yLength         (keepLogical);
        aedat.data.frame.xPosition         = aedat.data.frame.xPosition       (keepLogical);
        aedat.data.frame.yPosition         = aedat.data.frame.yPosition       (keepLogical);
        aedat.data.frame.numEvents         = nnz(keepLogical);
        if isfield(aedat.data.frame, 'valid')
            aedat.data.frame.valid    = aedat.data.frame.valid (keepLogical);
        end
    end
end
% Imu6

if isfield(aedat.data, 'imu6')
    ts = int64(aedat.data.imu6.timeStamp);
    lastTsReset = find(ts(2 : end) - ts(1 : end - 1) < 0, 1, 'last');
    if ~isempty(lastTsReset)
        keepLogical = true(aedat.data.imu6.numEvents, 1);
        keepLogical(1 : lastTsReset) = false;
        aedat.data.imu6.timeStamp = aedat.data.imu6.timeStamp (keepLogical);
        aedat.data.imu6.accelX       = aedat.data.imu6.accelX (keepLogical);
        aedat.data.imu6.accelY       = aedat.data.imu6.accelY (keepLogical);
        aedat.data.imu6.accelZ       = aedat.data.imu6.accelZ (keepLogical);
        aedat.data.imu6.gyroX        = aedat.data.imu6.gyroX (keepLogical);
        aedat.data.imu6.gyroY        = aedat.data.imu6.gyroY (keepLogical);
        aedat.data.imu6.gyroZ        = aedat.data.imu6.gyroZ (keepLogical);
        aedat.data.imu6.temperature  = aedat.data.imu6.temperature (keepLogical);
        aedat.data.imu6.numEvents    = nnz(keepLogical);
        if isfield(aedat.data.imu6, 'valid')
            aedat.data.imu6.valid    = aedat.data.imu6.valid (keepLogical);
        end
    end
end

% To do: handle other event types

% To do - correct first and last timestamp in info, also numPackets, ...
% startTime, endTime, startPacket, endPacket, startEvent, endEvent

author: panpanfei


text:
function [BAS1,BAS2] = basIFE3DNed1coef(vert,intpt,coef)

%% USAGE: Coefficients of P1 IFE Shape Function on a Tetrahedron
%
% INPUTS:
% vert --- 4-by-3 matrix stores the
%          vertices of the tetrahedron
%          [A1;A2;A3;A4] = [x1 y1 z1; x2 y2 z2; x3 y3 z3; x4 y4 z4]
% intpt --- 3-by-3 or 4-by-3 matrix stores three intersection points
%          [D;E;F] = [xd yd zd; xe ye ze; xf yf zf];
%           the order is 1,2,3 or 2,3,4,5 of the index of edges
% coef --- [bt1,bt2] jump coefficient where bt1 is associated with A1.
%                                    bt2 is associated with A2 A3 A4
%
%     Type 1 Interface Element (3 intersection points), it means
%         A1 is on piece #1 whose coefficient is coef(1)
%         A2 A3 A4 are on piece #2 whose coefficient is coef(2)
%
%     Type 2 Interface Element (4 intersection points), it means
%         A1 A2 are on piece #1 whose coefficient is coef(1)
%         A3 A4 are on piece #2 whose coefficient is coef(2)
%
% OUTPUTS:
% BAS1 --- 6-by-6 matrix stores the coefficients of shape fun on piece #1
% BAS2 --- 6-by-6 matrix stores the coefficients of shape fun on piece #2
%          The ith basis: a * x + b on piece #j = 1 or 2
%              1st basis [a11,a12,a13,b11,b12,b13]
%              2nd basis [a21,a22,a23,b21,b22,b23]
%              3rd basis [a31,a32,a33,b31,b32,b33]
%              4th basis [a41,a42,a43,b41,b42,b43]
%              5th basis [a51,a52,a53,b51,b52,b53]
%              6th basis [a61,a62,a63,b61,b62,b63]
%
% Last Modified by Ruchi Guo on 19/11/20

%% 0. Initialization
bt1 = coef(1);  bt2 = coef(2); rb = bt2/bt1;
at1 = coef(3);  at2 = coef(4); ra = at2/at1;
 M = zeros(6,6); l = zeros(6,2); % the 2nd component is only associated with the shape fun on the piece2
if size(intpt,1) == 3
    D = intpt(1,:); E = intpt(2,:); F = intpt(3,:); C = 1/3*(D+E+F);
    M(1,1:3) = 1/2*(vert(1,:)+D); l(1,1) = norm(D-vert(1,:));
    M(1,4:6) = 1/2*(vert(2,:)+D); l(1,2) = norm(D-vert(2,:));
    M(2,1:3) = 1/2*(vert(1,:)+E); l(2,1) = norm(E-vert(1,:));
    M(2,4:6) = 1/2*(vert(3,:)+E); l(2,2) = norm(E-vert(3,:));
    M(3,1:3) = 1/2*(vert(1,:)+F); l(3,1) = norm(F-vert(1,:));
    M(3,4:6) = 1/2*(vert(4,:)+F); l(3,2) = norm(F-vert(4,:));
    M(4,4:6) = 1/2*(vert(2,:)+vert(3,:)); l(4,2) = norm(vert(2,:)-vert(3,:));
    M(5,4:6) = 1/2*(vert(2,:)+vert(4,:)); l(5,2) = norm(vert(2,:)-vert(4,:));
    M(6,4:6) = 1/2*(vert(3,:)+vert(4,:)); l(6,2) = norm(vert(3,:)-vert(4,:));
elseif size(intpt,1) == 4
    [intptU,~,~] = vert4to3(intpt);
    D = intptU(1,:); E = intptU(2,:); F = intptU(3,:); C = 1/3*(D+E+F);
    M(1,1:3) = 1/2*(vert(1,:)+vert(2,:)); l(1,1) = norm(vert(1,:)-vert(2,:));
    M(2,1:3) = 1/2*(vert(1,:)+intpt(1,:)); l(2,1) = norm(vert(1,:)-intpt(1,:));
    M(2,4:6) = 1/2*(vert(3,:)+intpt(1,:)); l(2,2) = norm(vert(3,:)-intpt(1,:));
    M(3,1:3) = 1/2*(vert(1,:)+intpt(2,:)); l(3,1) = norm(vert(1,:)-intpt(2,:));
    M(3,4:6) = 1/2*(vert(4,:)+intpt(2,:)); l(3,2) = norm(vert(4,:)-intpt(2,:));
    M(4,1:3) = 1/2*(vert(2,:)+intpt(3,:)); l(4,1) = norm(vert(2,:)-intpt(3,:));
    M(4,4:6) = 1/2*(vert(3,:)+intpt(3,:)); l(4,2) = norm(vert(3,:)-intpt(3,:));
    M(5,1:3) = 1/2*(vert(2,:)+intpt(4,:)); l(5,1) = norm(vert(2,:)-intpt(4,:));
    M(5,4:6) = 1/2*(vert(4,:)+intpt(4,:)); l(5,2) = norm(vert(4,:)-intpt(4,:));
    M(6,4:6) = 1/2*(vert(3,:)+vert(4,:)); l(6,2) = norm(vert(3,:)-vert(4,:));
end

t = zeros(6,6);
t(1,1:3) = vert(2,:)-vert(1,:); t(1,1:3) = t(1,1:3)/norm(t(1,1:3));
t(2,1:3) = vert(3,:)-vert(1,:); t(2,1:3) = t(2,1:3)/norm(t(2,1:3));
t(3,1:3) = vert(4,:)-vert(1,:); t(3,1:3) = t(3,1:3)/norm(t(3,1:3));
t(4,1:3) = vert(3,:)-vert(2,:); t(4,1:3) = t(4,1:3)/norm(t(4,1:3));
t(5,1:3) = vert(4,:)-vert(2,:); t(5,1:3) = t(5,1:3)/norm(t(5,1:3));
t(6,1:3) = vert(4,:)-vert(3,:); t(6,1:3) = t(6,1:3)/norm(t(6,1:3));

n = cross(E-D,F-D);
n = n/norm(n); % unit normal of plane DEF
% L = @(x,y,z) dot([x-C(1),y-C(2),z-C(3)],n);

%% for the piece2
bas0 = zeros(6,6);
bas0(:,1) = (M(:,2).*t(:,3) - M(:,3).*t(:,2)).*l(:,1) + (M(:,5).*t(:,3) - M(:,6).*t(:,2)).*l(:,2);
bas0(:,2) = (M(:,3).*t(:,1) - M(:,1).*t(:,3)).*l(:,1) + (M(:,6).*t(:,1) - M(:,4).*t(:,3)).*l(:,2);
bas0(:,3) = (M(:,1).*t(:,2) - M(:,2).*t(:,1)).*l(:,1) + (M(:,4).*t(:,2) - M(:,5).*t(:,1)).*l(:,2);
bas0(:,4:6) = t(:,1:3).*((l(:,1)+l(:,2))*ones(1,3));

%% for the piece1
MM = M - [ones(6,1)*C, ones(6,1)*C];
bas1 = zeros(6,6);
bas1(:,1) = ((ra-1)*(MM(:,2).*t(:,3) - MM(:,3).*t(:,2)) +...
    (rb-1)*(t(:,1:3)*n')*(C(2)*n(3)-C(3)*n(2))).*l(:,1);
bas1(:,2) = ((ra-1)*(MM(:,3).*t(:,1) - MM(:,1).*t(:,3)) +...
    (rb-1)*(t(:,1:3)*n')*(C(3)*n(1)-C(1)*n(3))).*l(:,1);
bas1(:,3) = ((ra-1)*(MM(:,1).*t(:,2) - MM(:,2).*t(:,1)) +...
    (rb-1)*(t(:,1:3)*n')*(C(1)*n(2)-C(2)*n(1))).*l(:,1);
bas1(:,1:3) = bas1(:,1:3) -...
    (ra-1)*(sum(cross(ones(6,1)*n,MM(:,1:3)).*t(:,1:3),2).*l(:,1))*n;
bas1(:,4:6) = (rb-1)*(l(:,1).*(t(:,1:3)*n'))*n;

A = bas0 + bas1;

ss = A\eye(6);
ss = ss';
BAS2 = ss;
BAS1 = ss;
BAS1(:,1:3) = ss(:,1:3) + (ra-1)*(ss(:,1:3)-(ss(:,1:3)*n')*n);
BAS1(:,4:6) = ss(:,4:6) - (ra-1)*cross(ss(:,1:3)-(ss(:,1:3)*n')*n,ones(6,1)*C) +...
    (rb-1)*((cross(ss(:,1:3),ones(6,1)*C)+ss(:,4:6))*n')*n;
%BAS1 = BAS1'; BAS2 = BAS2';

return;

repo: ifem


text:
function [dist,proj] = ref_3d_projection(ref,EAST,NORTH,h)

% SYNTAX:
%   [dist,proj] = ref_3d_projection(ref,EAST,NORTH,h)
%
% INPUT:
%   ref = reference path (X,Y,Z coordinates of the vertices)
%   EAST = estimated trajectory in UTM coordinates (EAST)
%   NORTH = estimated trajectory in UTM coordinates (NORTH)
%   h = ellipsoid height
%
% OUTPUT:
%   dist = 3D distance of each estimated point from the reference
%   proj = projected trajectory
%
% DESCRIPTION:
%   3D projection on a reference path.
%   At the moment working only for adjacency matrix as in the following
%   example:
%
%   adj_mat = [ 0 1 0 0 0 1
%               1 0 1 0 0 0
%               0 1 0 1 0 0
%               0 0 1 0 1 0
%               0 0 0 1 0 1
%               1 0 0 0 1 0 ];

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

%close the path by connecting the first and the last points
ref = [ref;ref(1,:)];

%computation of the angular coefficients
ax = ref(2:end,1) - ref(1:end-1,1);
ay = ref(2:end,2) - ref(1:end-1,2);
az = ref(2:end,3) - ref(1:end-1,3);

%normalization on the segment distance
ad = sqrt(ax.^2 + ay.^2 + az.^2);
ax = ax ./ ad;
ay = ay ./ ad;
az = az ./ ad;

%offset of the curvilinear coordinate
s0 = [0; cumsum(ad)];

for j = 1 : length(EAST)
    pos_R(1) = EAST(j);
    pos_R(2) = NORTH(j);
    pos_R(3) = h(j);

    d0 = sqrt((pos_R(1) - ref(:,1)).^2 + ...
              (pos_R(2) - ref(:,2)).^2 + ...
              (pos_R(3) - ref(:,3)).^2);

    [dmin0 i0] = min(d0);

    %projection on the reference path
    bx = pos_R(1) - ref(1:end-1,1) + ax.*s0(1:end-1);
    by = pos_R(2) - ref(1:end-1,2) + ay.*s0(1:end-1);
    bz = pos_R(3) - ref(1:end-1,3) + az.*s0(1:end-1);

    s_R = (ax.*bx + ay.*by + az.*bz) ./ (ax.^2 + ay.^2 + az.^2);

    pos_R_proj(:,1) = ref(1:end-1,1) + ax .* (s_R - s0(1:end-1));
    pos_R_proj(:,2) = ref(1:end-1,2) + ay .* (s_R - s0(1:end-1));
    pos_R_proj(:,3) = ref(1:end-1,3) + az .* (s_R - s0(1:end-1));

    %computation of the minimum distance
    d = sqrt((pos_R(1) - pos_R_proj(:,1)).^2 + ...
        (pos_R(2) - pos_R_proj(:,2)).^2 + ...
        (pos_R(3) - pos_R_proj(:,3)).^2);

    [dmin i] = min(d);

    %position in cartesian coordinates
    while (dmin < dmin0) & ((pos_R_proj(i,1) < min(ref(i,1),ref(i+1,1))) | (pos_R_proj(i,1) > max(ref(i,1),ref(i+1,1))) | ...
            (pos_R_proj(i,2) < min(ref(i,2),ref(i+1,2))) | (pos_R_proj(i,2) > max(ref(i,2),ref(i+1,2))) | ...
            (pos_R_proj(i,3) < min(ref(i,3),ref(i+1,3))) | (pos_R_proj(i,3) > max(ref(i,3),ref(i+1,3))))

        d(i) = 9e99;
        [dmin i] = min(d);

    end

    if dmin0 < dmin
        dist(j,1) = dmin0;

        proj(j,1) = ref(i0,1);
        proj(j,2) = ref(i0,2);
        proj(j,3) = ref(i0,3);

    else
        dist(j,1) = dmin;

        proj(j,1) = pos_R_proj(i,1);
        proj(j,2) = pos_R_proj(i,2);
        proj(j,3) = pos_R_proj(i,3);
    end
end

author: goGPS-Project


text:
classdef testTopOptDesignAndShapes < testTopOptPrinting

    properties (Access = protected)
        fileOutputName = 'testTopOptDesignAndShapes';
        printMode = 'DesignAndShapes';
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function fd_to_tec ( prefix )

%*****************************************************************************80
%
%% FD_TO_TEC converts an FD model into a TECPLOT ASCII file.
%
%  Usage:
%
%    fd_to_tec ( 'prefix' )
%
%    where 'prefix' is the common prefix for the FD files:
%
%    * 'prefix'_nodes.txt,    the node coordinates.
%    * 'prefix'_values.txt,   the values defined at each node.
%
%    reads the data from these files,
%    creates a Delaunay triangulation or tetrahedralization of the data,
%    creates an equivalent TECPLOT ASCII file named 'prefix'.dat
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 March 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD_TO_TEC:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read an FD model consisting of:\n' );
  fprintf ( 1, '  * a nodes file (2D or 3D),\n' );
  fprintf ( 1, '  * a values file (1 or more values at each node);\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Write an equivalent TECPLOT ASCII "dat" file.\n' );
%
%  Get the file prefix.
%
  if ( nargin < 1 )

    fprintf ( 1, '\n' );
    prefix = input ( 'Enter the file prefix:' );

  end
%
%  Create the filenames.
%
  node_filename = strcat ( prefix, '_nodes.txt' );
  element_filename = strcat ( prefix, '_elements.txt' );
  value_filename = strcat ( prefix, '_values.txt' );
  tec_filename = strcat ( prefix, '.dat' );
%
%  Read the FD data (node coordinates, node values).
%
  [ dim_num, node_num ] = r8mat_header_read ( node_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The spatial dimension is %d\n', dim_num );
  fprintf ( 1, '  The number of nodes is %d\n', node_num );

  node_coord = r8mat_data_read ( node_filename, dim_num, node_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read node data from "%s".\n', node_filename );

  [ value_num, node_num2 ] = r8mat_header_read ( value_filename );

  if ( node_num2 ~= node_num )
    fprintf ( 1, '\n' );
    fprintf ( 1, ...
      '  The number of nodes in the node coordinate file is %d.\n', ...
      node_num );
    fprintf ( 1, ...
      '  but the number in the node data file is %d\n', node_num2 );
    fprintf ( 1, '  Because of this, no node data will be stored.\n' );
    value_num = 0;
  else
    fprintf ( 1, '  The number of values per node is %d\n', value_num );
    values = r8mat_data_read ( value_filename, value_num, node_num );
    fprintf ( 1, '  Read node value data from "%s".\n', value_filename );
  end
%
%  Triangulate/Tetrahedralize the data.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Calling DELAUNAYN to create elements.\n' );

  element_node = delaunayn ( node_coord' );
  element_node = element_node';
  [ element_order, element_num ] = size ( element_node );

  fprintf ( 1, '  Number of elements is %d\n', element_num );

  i4mat_write ( element_filename, element_order, element_num, element_node );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wrote element data to "%s"\n', element_filename );
%
%  Write the TEC data.
%
  tec_write ( tec_filename, dim_num, node_num, element_num, ...
    element_order, value_num, node_coord, element_node, values );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wrote the TEC data to "%s".\n', tec_filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD_TO_TEC\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

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

  return
end
function value = ch_is_alpha ( c )

%*****************************************************************************80
%
%% CH_IS_ALPHA returns TRUE if C is an alphabetic character.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, a character to check.
%
%    Output, logical VALUE is TRUE if C is an alphabetic character.
%
  FALSE = 0;
  TRUE = 1;

  if ( ( 'a' <= c & c <= 'z' ) | ...
       ( 'A' <= c & c <= 'Z' ) )
    value = TRUE;
  else
    value = FALSE;
  end

  return
end
function [ dim_num, node_num, element_num, element_order, value_num, ...
  node_coord, element_node, values ] = fem_read ( node_filename, ...
  element_filename, value_filename )

%*****************************************************************************80
%
%% FEM_READ reads data files associated with a finite element solution.
%
%  Discussion:
%
%    This program reads the node, element and data files that define
%    a finite element geometry and data based on that geometry:
%    * a set of nodes,
%    * a set of elements based on those nodes,
%    * a set of data values associated with each node.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string NODE_FILENAME, the name of the node coordinate
%    file.  If this argument is not supplied, it will be requested.  If the interactive
%    response is blank, or otherwise defective, then the program terminates.
%
%    Input, string ELEMENT_FILENAME, the name of the element file.  If
%    this argument is not supplied, it will be requested.  If the interactive
%    response is blank, then the program will assume that no element information
%    is to be supplied.  (But the node coordinates must be available and may be plotted.
%    And if a node data file is supplied, then the data can be plotted against
%    the node coordinates without using any finite element structure.)
%
%    Input, string VALUES_FILENAME, the name of the node data file.  If
%    this argument is not supplied, it will be requested.  If the interactive
%    response is blank, then the program will assume that no node data information
%    is to be supplied.  (But the node coordinates will be available and may be plotted.
%    And if an element file is supplied, then the elements can also be displayed.)
%
%    Output, integer DIM_NUM, the spatial dimension, inferred from the
%    "shape" of the data in the node file.
%
%    Output, integer NODE_NUM, the number of nodes, inferred from the
%    number of lines of data in the node coordinate file.
%
%    Output, integer ELEMENT_NUM, the number of elements, inferred from the
%    number of lines of data in the element file.
%
%    Output, integer ELEMENT_ORDER, the order of the elements, inferred from
%    the number of items in the first line of the element file.
%
%    Output, integer VALUE_NUM, the number of data items per node,
%    inferred from the number of items in the first line of the node data file.
%
%    Output, real NODE_COORD(DIM_NUM,NODE_NUM), the coordinates of nodes.
%
%    Output, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM);
%    the global index of local node I in element J.
%
%    Output, real VALUES(VALUE_NUM,NODE_NUM), the data values associated
%    with each node.
%

%
%  Read the node coordinate file.
%
  [ dim_num, node_num ] = r8mat_header_read ( node_filename );

  [ element_order, element_num ] = i4mat_header_read ( element_filename );

  [ value_num, node_num2 ] = r8mat_header_read ( value_filename );

  if ( node_num2 ~= node_num )
    fprintf ( 1, '\n' );
    fprintf ( 1, ...
      '  The number of nodes in the node coordinate file is %d.\n', ...
      node_num );
    fprintf ( 1, ...
      '  but the number in the node data file is %d\n', node_num2 );
    fprintf ( 1, '  Because of this, no node data will be stored.\n' );
    value_num = 0;
  end

  node_coord = r8mat_data_read ( node_filename, dim_num, node_num );

  element_node = i4mat_data_read ( element_filename, element_order, ...
    element_num );

  values = r8mat_data_read ( value_filename, value_num, node_num );

  return
end
function column_num = file_column_count ( input_filename )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_filename );
    error ( 'FILE_COLUMN_COUNT - Error!' );
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return;
  end

  column_num = s_word_count ( line );

  return
end
function filename = file_name_inc ( filename )

%*****************************************************************************80
%
% FILE_NAME_INC generates the next filename in a series.
%
%  Discussion:
%
%    It is assumed that the digits in the name, whether scattered or
%    connected, represent a number that is to be increased by 1 on
%    each call.  If this number is all 9's on input, the output number
%    is all 0's.  Non-numeric letters of the name are unaffected..
%
%    If the name is empty, then the routine stops.
%
%    If the name contains no digits, the empty string is returned.
%
%  Example:
%
%      Input            Output
%      -----            ------
%      'a7to11.txt'     'a7to12.txt'  (typical case.  Last digit incremented)
%      'a7to99.txt'     'a8to00.txt'  (last digit incremented, with carry.)
%      'a9to99.txt'     'a0to00.txt'  (wrap around)
%      'cat.txt'        ' '           (no digits in input name.)
%      ' '              STOP!         (error.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILENAME, the string to be incremented.
%
%    Output, string FILENAME, the incremented string.
%
  lens = length ( filename );

  if ( lens <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_NAME_INC - Fatal error!\n' );
    fprintf ( 1, '  The input filename is empty.\n' );
    error ( 'FILE_NAME_INC - Fatal error!' );
  end

  change = 0;

  for i = lens : -1 : 1

    c = filename(i);

    if ( '0' <= c & c <= '8' )

      change = change + 1;

      c = c + 1;

      filename(i) = c;

      return

    elseif ( c == '9' )

      change = change + 1;

      c = '0';

      filename(i) = c;

    end

  end

  if ( change == 0 )
    filename = ' ';
  end

  return
end
function row_num = file_row_count ( input_filename )

%*****************************************************************************80
%
%% FILE_ROW_COUNT counts the number of row records in a file.
%
%  Discussion:
%
%    Each input line is a "RECORD".
%
%    The records are divided into three groups:
%
%    * BLANK LINES (nothing but blanks)
%    * COMMENT LINES (begin with a '#')
%    * DATA RECORDS (anything else)
%
%    The value returned by the function is the number of data records.
%
%    By the way, if the MATLAB routine FGETS is used, instead of
%    FGETL, then the variable LINE will include line termination
%    characters, which means that a blank line would not actually
%    have zero characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer ROW_NUM, the number of rows found.
%
  input_unit = fopen ( input_filename );

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

  blank_num = 0;
  comment_num = 0;
  row_num = 0;

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    record_num = record_num + 1;
    record_length = s_len_trim ( line );

    if ( record_length <= 0 )
      blank_num = blank_num + 1;
    elseif ( line(1) == '#' )
      comment_num = comment_num + 1;
    else
      row_num = row_num + 1;
    end

  end

  fclose ( input_unit );

  return
end
function i4mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%
%    Input, logical HEADER, is TRUE if the header is to be included.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function table = r8mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %f' );
  end

  input_unit = fopen ( input_filename );

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

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = r8mat_header_read ( input_filename )

%*****************************************************************************80
%
%% R8MAT_HEADER_READ reads the header from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function iloc = s_alpha_last ( s )

%*****************************************************************************80
%
%% S_ALPHA_LAST returns the location of the last alphabetic character.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be searched.
%
%    Output, integer ILOC, the location of the last alphabetic
%    character in the string.  If there are no alphabetic
%    characters, ILOC is returned as -1.
%
  s_len = s_len_trim ( s );

  for i = s_len : -1 : 1
    if ( ch_is_alpha ( s(i) ) )
      iloc = i;
      return
    end
  end

  iloc = -1;

  return
end
function [ s, ierror ] = s_inc ( s )

%*****************************************************************************80
%
%% S_INC "increments" a string.
%
%  Discussion:
%
%    The routine tries to produce the next string, in dictionary order,
%    following the input value of a string.  Digits, spaces, and other
%    nonalphabetic characters are ignored.  Case is respected; in other
%    words, the case of every alphabetic character on input will be the
%    same on output.
%
%    The following error conditions can occur:
%
%      There are no alphabetic characters in the string.  No
%      incrementing is possible.
%
%      All alphabetic characters are equal to 'Z' or 'z'.  In this
%      case, an error value is returned, but the string is also "wrapped
%      around" so that all alphabetic characters are "A" or "a".
%
%    If the word "Tax" were input, the successive outputs would be
%    "Tay", "Taz", "Tba", "Tbb", ...  If the input word "January 4, 1989"
%    were input, the output would be "Januarz 4, 1989".
%
%    This routine could be useful when trying to create a unique file
%    name or variable name at run time.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character ( len = * ) S, the string whose
%    alphabetic successor is desired.
%
%    Output, character ( len = * ) S, if IERROR = 0,
%    S has been replaced by its successor.  If IERROR = 2,
%    S will be "wrapped around" so that all alphabetic
%    characters equal "A" or "a".
%
%    Output, integer IERROR, an error flag.
%    0, no error.
%    1, no alphabetic characters occur in the string.
%    2, all alphabetic characters are "Z" or "z".  S is wrapped around so
%       that all alphabetic characters are "A" or "a".
%
  ierror = 0;
  ilo = 1;
  ihi = s_len_trim ( s );
%
%  Find the last alphabetic character in the string.
%
  while ( 1 )

    iloc = s_alpha_last ( s(ilo:ihi) );
%
%  If there is no alphabetic character, we can't help.
%
    if ( iloc == 0 )
      ierror = 1;
      break;
    end

    if ( s(iloc) == 'Z' )

      s(iloc) = 'A';
      ihi = iloc - 1;

      if ( ihi <= 0 )
        ierror = 2;
        break;
      end

    elseif ( s(iloc) == 'z' )

      s(iloc) = 'a';
      ihi = iloc - 1;

      if ( ihi <= 0 )
        ierror = 2;
        break;
      end

    else

      s(iloc) = s(iloc) + 1;
      break;

    end

  end

  return
end
function len = s_len_trim ( s )

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

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

  return
end
function word_num = s_word_count ( s )

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

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

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

  end

  return
end
function tec_write ( tec_filename, dim_num, node_num, element_num, ...
  element_order, value_num, node_coord, element_node, values )

%*****************************************************************************80
%
%% TEC_WRITE writes finite element data to a TEC file.
%
%  Discussion:
%
%    This program writes the node, element and data files that define
%    a finite element geometry and data based on that geometry:
%    * a set of nodes,
%    * a set of elements based on those nodes,
%    * a set of data values associated with each node.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string NODE_FILENAME, the name of the node coordinate
%    file.  If this argument is empty, no node coordinate file will be written.
%
%    Input, string ELEMENT_FILENAME, the name of the element file.  If
%    this argument is empty, no element file will be written.
%
%    Input, string VALUES_FILENAME, the name of the node data file.  If
%    this argument is empty, no node data file will be written.
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer VALUE_NUM, the number of data items per node.
%
%    Input, real NODE_COORD(DIM_NUM,NODE_NUM), the coordinates of nodes.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM);
%    the global index of local node I in element J.
%
%    Input, real VALUES(VALUE_NUM,NODE_NUM), the data values associated
%    with each node.
%
  tec_file_unit = fopen ( tec_filename, 'wt' );
%
%  Write the title.
%
  fprintf ( tec_file_unit, 'TITLE = "%s"\n', tec_filename );
%
%  Write the variable names.
%
  fprintf ( tec_file_unit, 'VARIABLES = ' );

  name = 'X';
  for dim = 1 : dim_num
    if ( dim == 1 )
      fprintf ( tec_file_unit, '"%s"', name );
    else
      fprintf ( tec_file_unit, ', "%s"', name );
    end
    name = s_inc ( name );
  end

  name = 'data_001';
  for dim = 1 : value_num
    fprintf ( tec_file_unit, ', "%s"', name );
    name = file_name_inc ( name );
  end

  fprintf ( tec_file_unit, '\n' );
%
%  Write the ZONE record.
%
  if ( dim_num == 2 & element_order == 3 )
    zonetype = 'FETRIANGLE';
  elseif ( dim_num == 2 & element_order == 4 )
    zonetype = 'FEQUADRILATERAL';
  elseif ( dim_num == 3 & element_order == 4 )
    zonetype = 'FETETRAHEDRON';
  elseif ( dim_num == 3 & element_order == 8 )
    zonetype = 'FEBRICK';
  else
    zonetype = 'FEUNKNOWN';
  end

  fprintf ( tec_file_unit, 'ZONE N = %d, E = %d, ', node_num, element_num );
  fprintf ( tec_file_unit, 'DATAPACKING = POINT, ZONETYPE = %s\n', ...
    zonetype );
%
%  Write the node coordinates and node data.
%
  for node = 1 : node_num
    for dim = 1 : dim_num
      fprintf ( tec_file_unit, '  %10f', node_coord(dim,node) );
    end
    for data = 1 : value_num
      fprintf ( tec_file_unit, '  %10f', values(data,node) );
    end
    fprintf ( tec_file_unit, '\n' );
  end
%
%  Write the element-node connectivity.
%
  for element = 1 : element_num
    for order = 1 : element_order
      fprintf ( tec_file_unit, '  %6d', element_node(order,element) );
    end
    fprintf ( tec_file_unit, '\n' );
  end

  fclose ( tec_file_unit );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEC_WRITE wrote all data to "%s".\n', ...
    tec_filename );

  return
end
function timestamp ( )

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

  return
end


text:
function M = getmassmatvec3(elem2dof,volume,Dlambda,elemType,K)
%% GETMASSMATVEC3 get a mass matrix of finite element spaces in 3D
%
% M = GETMASSMATVEC3(elem2dof,volume,Dlambda,elemType,K) get mass matrix of
% the finite element space specified by elemType. The coefficient K is
% piecewise constant.
%
% The elemType can be:
%
% - "RT0": The lowest order Raviart-Thomas element
% - "ND0": The lowest order Nedelec element
% - "ND1": The full linear Nedelec element
% - "ND2m": The incomplete quadratic Nedelec element (P\Lambda_2^{-})
% - "ND2": The full quadratic Nedelec element (P\Lambda_2)
%
% Note that for H(div) element, the coefficient is Mij/K while for H(curl)
% element, it is Mij*K.
%
% See also Maxwell, Maxwell1, Maxwell2, Poisson3RT0
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('K','var'), K = []; end
NT = size(elem2dof,1);

if contains(elemType, 'ND')
    DiDj = zeros(NT,4,4);
    for i = 1:4
        for j = i:4
            DiDj(:,i,j) = dot(Dlambda(:,:,i),Dlambda(:,:,j),2);
            DiDj(:,j,i) = DiDj(:,i,j);
        end
    end
end

%% ND0: the lowest order edge element
if strcmp(elemType,'ND0')
    NE = double(max(elem2dof(:)));
    locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
    ii = zeros(21*NT,1); jj = zeros(21*NT,1); sM = zeros(21*NT,1);
    index = 0;
    for i = 1:6
        for j = i:6
            ii(index+1:index+NT) = double(elem2dof(:,i));
            jj(index+1:index+NT) = double(elem2dof(:,j));
            % mass matrix
            % locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
            i1 = locEdge(i,1); i2 = locEdge(i,2);
            j1 = locEdge(j,1); j2 = locEdge(j,2);
            Mij = 1/20*volume.*( (1+(i1==j1))*DiDj(:,i2,j2) ...
                               - (1+(i1==j2))*DiDj(:,i2,j1) ...
                               - (1+(i2==j1))*DiDj(:,i1,j2) ...
                               + (1+(i2==j2))*DiDj(:,i1,j1));
            if ~isempty(K)
                Mij = Mij.*K;
            end
            sM(index+1:index+NT) = Mij;
            index = index + NT;
        end
    end
    diagIdx = (ii == jj);   upperIdx = ~diagIdx;
    M = sparse(ii(diagIdx),jj(diagIdx),sM(diagIdx),NE,NE);
    MU = sparse(ii(upperIdx),jj(upperIdx),sM(upperIdx),NE,NE);
    M = M + MU + MU';
end

%% ND1: The full linear Nedelec element
if strcmp(elemType,'ND1')
    NE = double(max(elem2dof(:)));
    Ndof = 2*NE;
    locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
    ii = zeros(21*NT,1); jj = zeros(21*NT,1);
    sMphi = zeros(21*NT,1); sMpsi = zeros(21*NT,1);
    index = 0;
    % two block diagonal matrix for phi and psi
    for i = 1:6
        for j = i:6
            % (phi_i,phi_j)
            ii(index+1:index+NT) = double(elem2dof(:,i));
            jj(index+1:index+NT) = double(elem2dof(:,j));
            % mass matrix
            % locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
            i1 = locEdge(i,1); i2 = locEdge(i,2);
            j1 = locEdge(j,1); j2 = locEdge(j,2);
            Mij = 1/20*volume.*( (1+(i1==j1))*DiDj(:,i2,j2) ...
                               - (1+(i1==j2))*DiDj(:,i2,j1) ...
                               - (1+(i2==j1))*DiDj(:,i1,j2) ...
                               + (1+(i2==j2))*DiDj(:,i1,j1));
            if ~isempty(K)
                Mij = Mij.*K;
            end
            sMphi(index+1:index+NT) = Mij;
            % (psi_i,psi_j)
            Mij = 1/20*volume.*( (1+(i1==j1))*DiDj(:,i2,j2) ...
                               + (1+(i1==j2))*DiDj(:,i2,j1) ...
                               + (1+(i2==j1))*DiDj(:,i1,j2) ...
                               + (1+(i2==j2))*DiDj(:,i1,j1));
            if ~isempty(K)
                Mij = Mij.*K;
            end
            sMpsi(index+1:index+NT) = Mij;
            index = index + NT;
        end
    end
    diagIdx = (ii == jj);   upperIdx = ~diagIdx;
    M = sparse([ii(diagIdx); ii(diagIdx)+NE], [jj(diagIdx); jj(diagIdx)+NE],...
               [sMphi(diagIdx); sMpsi(diagIdx)],Ndof,Ndof);
    MU = sparse([ii(upperIdx); ii(upperIdx)+NE], [jj(upperIdx); jj(upperIdx)+NE],...
               [sMphi(upperIdx); sMpsi(upperIdx)],Ndof,Ndof);
    M = M + MU + MU';
    % off-diagonal matrix (psi, phi)
    ii = zeros(36*NT,1); jj = zeros(36*NT,1); ss = zeros(36*NT,1);
    index = 0;
    for i = 1:6
        for j = 1:6
            % local to global index map and its sign
            i1 = locEdge(i,1); i2 = locEdge(i,2);
            j1 = locEdge(j,1); j2 = locEdge(j,2);
            Mij = 1/20*volume.*( (1+(i1==j1))*DiDj(:,i2,j2) ...
                       - (1+(i1==j2))*DiDj(:,i2,j1) ...
                       + (1+(i2==j1))*DiDj(:,i1,j2) ...
                       - (1+(i2==j2))*DiDj(:,i1,j1));
            if ~isempty(K)
                Mij = Mij.*K;
            end
            ii(index+1:index+NT) = double(elem2dof(:,i))+NE;
            jj(index+1:index+NT) = double(elem2dof(:,j));
            ss(index+1:index+NT) = Mij;
            index = index + NT;
        end
    end
    ML = sparse(ii,jj,ss,Ndof,Ndof);
    M  = M  + ML + ML';
end

%% ND2m: the incomplete quadratic Nedelec element
% see Maxwell2
% elem2dof = [elem2edge elem2edge+NE elem2face+2*NE elem2face+2*NE+NF];

if strcmp(elemType,'ND2m')
    NE = double(max(elem2dof(:, 1:6), [], "all"));
    NF = double(max(elem2dof(:, 12:16)-2*NE, [], "all"));
    Ndof = 2*(NE+NF);
    locBasesIdx = [1 2 0; 1 3 0; 1 4 0; 2 3 0; 2 4 0; 3 4 0; ... % Nd0
               1 2 0; 1 3 0; 1 4 0; 2 3 0; 2 4 0; 3 4 0; ... % grad(P2)
               3 2 4; 3 1 4; 2 1 4; 2 1 3; ...
               4 2 3; 4 1 3; 4 1 2; 3 1 2]; % face bubbles
    ii = zeros(210*NT,1); jj = zeros(210*NT,1);
    sM = zeros(210*NT,1);
    index = 0;
    for i = 1:20
        for j = i:20
            ii(index+1:index+NT) = double(elem2dof(:,i));
            jj(index+1:index+NT) = double(elem2dof(:,j));
            i1 = locBasesIdx(i,1); i2 = locBasesIdx(i,2); i3 = locBasesIdx(i,3);
            j1 = locBasesIdx(j,1); j2 = locBasesIdx(j,2); j3 = locBasesIdx(j,3);
            Mij = zeros(NT,1);
            if (i<=6) && (j<=6)
                % block 1: (phi_i,phi_j)
                Mij = 1/20*((1+(i1==j1))*DiDj(:,i2,j2) ...
                    - (1+(i1==j2))*DiDj(:,i2,j1) ...
                    - (1+(i2==j1))*DiDj(:,i1,j2) ...
                    + (1+(i2==j2))*DiDj(:,i1,j1));
            end
            if (i<=6) && (7<=j) && (j<=12)
                % block 2: (psi_j,phi_i)
                Mij = 1/20*( (1+(j1==i1))*DiDj(:,j2,i2) ...
                    - (1+(j1==i2))*DiDj(:,j2,i1) ...
                    + (1+(j2==i1))*DiDj(:,j1,i2) ...
                    - (1+(j2==i2))*DiDj(:,j1,i1));

            end
            if (7<=i) && (i<=12) && (7<=j) && (j<=12)
                % block 3: (psi_j,psi_i)
                Mij = 1/20*((1+(i1==j1))*DiDj(:,i2,j2) ...
                    + (1+(i1==j2))*DiDj(:,i2,j1) ...
                    + (1+(i2==j1))*DiDj(:,i1,j2) ...
                    + (1+(i2==j2))*DiDj(:,i1,j1));
            end
            if (i<=6) && (j>12)
                % block 4: (chi_j,phi_i)
                Mij = intlambda([j1,i1,j2],3)*DiDj(:,i2,j3) ...
                    -intlambda([j1,i1,j3],3)*DiDj(:,i2,j2) ...
                    -intlambda([j1,i2,j2],3)*DiDj(:,i1,j3) ...
                    +intlambda([j1,i2,j3],3)*DiDj(:,i1,j2);
            end
            if (7<=i) && (i<=12) && (j>12)
                % block 5: (chi_j,psi_i)
                Mij = intlambda([j1,i1,j2],3)*DiDj(:,i2,j3) ...
                    -intlambda([j1,i1,j3],3)*DiDj(:,i2,j2) ...
                    +intlambda([j1,i2,j2],3)*DiDj(:,i1,j3) ...
                    -intlambda([j1,i2,j3],3)*DiDj(:,i1,j2);
            end
            if (i>12) && (j>12)
                % block 6: (chi_j,chi_i)
                Mij = intlambda([i1,j1,i2,j2],3)*DiDj(:,i3,j3) ...
                    -intlambda([i1,j1,i2,j3],3)*DiDj(:,i3,j2) ...
                    -intlambda([i1,j1,i3,j2],3)*DiDj(:,i2,j3) ...
                    +intlambda([i1,j1,i3,j3],3)*DiDj(:,i2,j2);
            end
            if ~isempty(K)
                Mij = Mij.*K;
            end
            Mij = Mij.*volume;
            sM(index+1:index+NT) = Mij;
            index = index + NT;
        end
    end

    diagIdx = (ii == jj);   upperIdx = ~diagIdx;
    M = sparse(ii(diagIdx),jj(diagIdx),sM(diagIdx),Ndof,Ndof);
    MU = sparse(ii(upperIdx),jj(upperIdx),sM(upperIdx),Ndof,Ndof);
    M = M + MU + MU';
end

%% ND2: the complete quadratic Nedelec element
% elem2dof = [elem2edge elem2edge+NE elem2edge+2*NE ...
%             elem2face+3*NE elem2face+3*NE+NF elem2face+3*NE+2*NF]
% index 1 to 18: edge basis
% index 19 to 30: face basis
% phi_1 = lambda_i^2 \nabla lambda_j
% phi_2 = lambda_j^2 \nabla lambda_i
% phi_3 = lambda_i lambda_j \nabla (lambda_i - lambda_j)
% psi_{F_{ijk}; 1,2,3} = ijk \in cyc {\lambda_i \lambda_j \nabla \lambda_k}

if strcmp(elemType,'ND2')
    NE = double(max(elem2dof(:, 1:6), [], "all"));
    NF = double(max(elem2dof(:, 19:22)-3*NE, [], "all"));
    Ndof = 3*(NE+NF);
    locBasesIdx = [1 2 0; 1 3 0; 1 4 0; 2 3 0; 2 4 0; 3 4 0; ... % phi_1
               2 1 0; 3 1 0; 4 1 0; 3 2 0; 4 2 0; 4 3 0; ... % psi_2
               1 2 0; 1 3 0; 1 4 0; 2 3 0; 2 4 0; 3 4 0; ... % psi_3
               2 3 4; 3 4 2; 4 2 3; ... % face 1
               1 3 4; 3 4 1; 4 1 3; ... % face 2
               1 2 4; 2 4 1; 4 1 2; ... % face 3
               1 2 3; 2 3 1; 3 1 2]; % face 4
    ii = zeros(465*NT,1); jj = zeros(465*NT,1);
    sM = zeros(465*NT,1);
    index = 0;
    for i = 1:30
        for j = i:30
            ii(index+1:index+NT) = double(elem2dof(:,i));
            jj(index+1:index+NT) = double(elem2dof(:,j));
            i1 = locBasesIdx(i,1); i2 = locBasesIdx(i,2); i3 = locBasesIdx(i,3);
            j1 = locBasesIdx(j,1); j2 = locBasesIdx(j,2); j3 = locBasesIdx(j,3);
            Mij = zeros(NT,1);
            if (i<=12) && (j<=12)
                % block 1: (phi_{1,2},phi_{1,2})
                Mij = intlambda([i1,i1,j1,j1],3)*DiDj(:,i2,j2);
            end
            if (i<=12) && (12<j) && (j<=18)
                % block 2: (phi_{1,2},phi_{3})
                Mij = intlambda([i1,i1,j1,j2],3)*...
                    (DiDj(:,i2,j1) - DiDj(:,i2,j2));
            end
            if (12<i) && (i<=18) && (12<j) && (j<=18)
                % block 3: (phi_{3},phi_{3})
                Mij = intlambda([i1,i2,j1,j2],3)* ...
                    (DiDj(:,i1,j1) + DiDj(:,i2,j2) ...
                    -DiDj(:,i1,j2) - DiDj(:,i2,j1));
            end
            if (i<=12) && (j>18)
                % block 4: (phi_{1,2}, psi_F)
                Mij = intlambda([i1,i1,j1,j2],3)*DiDj(:,i2,j3);
            end
            if (12<i) && (i<=18) && (j>18)
                % block 5: (phi_{3}, psi_F)
                Mij = intlambda([i1,i2,j1,j2],3)*...
                    (DiDj(:,i1,j3) - DiDj(:,i2,j3));
            end
            if (i>18) && (j>18)
                % block 6: (psi_F,psi_F)
                Mij = intlambda([i1,i2,j1,j2],3)*DiDj(:,i3,j3);
            end
            if ~isempty(K)
                Mij = Mij.*K;
            end
            Mij = Mij.*volume;
            sM(index+1:index+NT) = Mij;
            index = index + NT;
        end
    end

    diagIdx = (ii == jj);   upperIdx = ~diagIdx;
    M = sparse(ii(diagIdx),jj(diagIdx),sM(diagIdx),Ndof,Ndof);
    MU = sparse(ii(upperIdx),jj(upperIdx),sM(upperIdx),Ndof,Ndof);
    M = M + MU + MU';
end

%% RT0: the lowest order face element
if strcmp(elemType,'RT0')
    NF = double(max(elem2dof(:)));
    localFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3]; % ascend ordering
    M = sparse(NF,NF);
    for i = 1:4
        for j = i:4
            % local to global index map
            ii = double(elem2dof(:,i));
            jj = double(elem2dof(:,j));
            i1 = localFace(i,1); i2 = localFace(i,2); i3 = localFace(i,3);
            j1 = localFace(j,1); j2 = localFace(j,2); j3 = localFace(j,3);
            % computation of mass matrix --- (phi_i, phi_j)
            Mij = 1/5*volume.*( ...
                  (1+(i1==j1))*dot(mycross(Dlambda(:,:,i2),Dlambda(:,:,i3),2), ...
                                   mycross(Dlambda(:,:,j2),Dlambda(:,:,j3),2),2)...
                 +(1+(i1==j2))*dot(mycross(Dlambda(:,:,i2),Dlambda(:,:,i3),2), ...
                                   mycross(Dlambda(:,:,j3),Dlambda(:,:,j1),2),2)...
                 +(1+(i1==j3))*dot(mycross(Dlambda(:,:,i2),Dlambda(:,:,i3),2), ...
                                   mycross(Dlambda(:,:,j1),Dlambda(:,:,j2),2),2)...
                 +(1+(i2==j1))*dot(mycross(Dlambda(:,:,i3),Dlambda(:,:,i1),2), ...
                                   mycross(Dlambda(:,:,j2),Dlambda(:,:,j3),2),2)...
                 +(1+(i2==j2))*dot(mycross(Dlambda(:,:,i3),Dlambda(:,:,i1),2), ...
                                   mycross(Dlambda(:,:,j3),Dlambda(:,:,j1),2),2)...
                 +(1+(i2==j3))*dot(mycross(Dlambda(:,:,i3),Dlambda(:,:,i1),2), ...
                                   mycross(Dlambda(:,:,j1),Dlambda(:,:,j2),2),2)...
                 +(1+(i3==j1))*dot(mycross(Dlambda(:,:,i1),Dlambda(:,:,i2),2), ...
                                   mycross(Dlambda(:,:,j2),Dlambda(:,:,j3),2),2)...
                 +(1+(i3==j2))*dot(mycross(Dlambda(:,:,i1),Dlambda(:,:,i2),2), ...
                                   mycross(Dlambda(:,:,j3),Dlambda(:,:,j1),2),2)...
                 +(1+(i3==j3))*dot(mycross(Dlambda(:,:,i1),Dlambda(:,:,i2),2), ...
                                   mycross(Dlambda(:,:,j1),Dlambda(:,:,j2),2),2));
            if ~isempty(K)
                Mij = Mij./K;
            end
            if (j==i)
                M = M + sparse(ii,jj,Mij,NF,NF);
            else
                M = M + sparse([ii;jj],[jj;ii],[Mij; Mij],NF,NF);
            end
        end
    end
end


text:
% OP_UXN_VXN_2D: assemble the matrix M = [m(i,j)], m(i,j) = (mu u_j x n, v_i x n), with n the exterior normal vector.
%
%   mat = op_uxn_vxn_2d (spu, spv, msh, coeff);
%   [rows, cols, values] = op_uxn_vxn_2d (spu, spv, msh, coeff);
%
% INPUT:
%
%  spu:   structure representing the space of trial functions (see sp_vector/sp_eval_boundary_side)
%  spv:   structure representing the space of test functions  (see sp_vector/sp_eval_boundary_side)
%  msh:   structure containing the domain partition and the quadrature rule (see msh_cartesian/msh_eval_boundary_side)
%  coeff: physical parameter
%
% OUTPUT:
%
%   mat:    assembled matrix
%   rows:   row indices of the nonzero entries
%   cols:   column indices of the nonzero entries
%   values: values of the nonzero entries
%
% Copyright (C) 2009, 2010 Carlo de Falco, Rafael Vazquez
% Copyright (C) 2011, 2017 Rafael Vazquez
%
%    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/>.

function varargout = op_uxn_vxn_2d (spu, spv, msh, coeff)

  rows = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  cols = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  values = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);

  jacdet_weights = msh.jacdet .* msh.quad_weights .* coeff;

  ncounter = 0;
  for iel = 1:msh.nel
    if (all (msh.jacdet(:, iel)))
      shpu_iel = reshape (spu.shape_functions(:, :, :, iel), spu.ncomp, msh.nqn, spu.nsh_max);
      shpv_iel = reshape (spv.shape_functions(:, :, :, iel), spv.ncomp, msh.nqn, spv.nsh_max);

      normal_iel = reshape (msh.normal(:,:,iel), 2, msh.nqn);
      shpu_x_n = bsxfun (@times, shpu_iel(1,:,:), normal_iel(2,:)) - ...
                 bsxfun (@times, shpu_iel(2,:,:), normal_iel(1,:));
      shpv_x_n = bsxfun (@times, shpv_iel(1,:,:), normal_iel(2,:)) - ...
                 bsxfun (@times, shpv_iel(2,:,:), normal_iel(1,:));
      shpu_x_n = reshape (shpu_x_n, msh.nqn, 1, spu.nsh_max);
      shpv_x_n = reshape (shpv_x_n, msh.nqn, spv.nsh_max, 1);

      jacdet_iel = jacdet_weights(:,iel);
      jacdet_shpu = bsxfun (@times, jacdet_iel, shpu_x_n);
      tmp1 = bsxfun (@times, jacdet_shpu, shpv_x_n);
      elementary_values = reshape (sum (tmp1, 1), spv.nsh_max, spu.nsh_max);

      [rows_loc, cols_loc] = ndgrid (spv.connectivity(:,iel), spu.connectivity(:,iel));
      indices = rows_loc & cols_loc;
      rows(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = rows_loc(indices);
      cols(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = cols_loc(indices);
      values(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = elementary_values(indices);
      ncounter = ncounter + spu.nsh(iel)*spv.nsh(iel);

    else
      warning ('geopdes:jacdet_zero_at_quad_node', 'op_uxn_vxn_2d: singular map in element number %d', iel)
    end
  end

  if (nargout == 1 || nargout == 0)
    varargout{1} = sparse (rows(1:ncounter), cols(1:ncounter), ...
                           values(1:ncounter), spv.ndof, spu.ndof);
  elseif (nargout == 3)
    varargout{1} = rows(1:ncounter);
    varargout{2} = cols(1:ncounter);
    varargout{3} = values(1:ncounter);
  else
    error ('op_uxn_vxn_2d: wrong number of output arguments')
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function quad = gauss_legendre_integrate_fast ( f, n )

%*****************************************************************************80
%
%% GAUSS_LEGENDRE_INTEGRATE_FAST applies a Gauss-Legendre quadrature rule.
%
%  Discussion:
%
%    The integration is carried out on the standard interval [-1,1].
%
%    The computation should be very efficient in MATLAB.
%
%  Modified:
%
%    03 March 2007
%
%  Author:
%
%    Lloyd Trefethen
%
%  Reference:
%
%    Lloyd Trefethen,
%    Is Gauss Quadrature Better than Clenshaw-Curtis?
%    SIAM Review,
%    Volume 50, Number 1, 2008, pages 67-87.
%
%  Parameters:
%
%    Input, function F ( x ), the function to be integrated.
%
%    Input, integer N, the order of the quadrature rule.
%
%    Output, real QUAD, the estimate for the integral of F(X).
%
  beta = 0.5 ./ sqrt ( 1 - ( 2 * ( 1:n ) ) .^ (-2) );

  tridiag = diag ( beta, 1 ) + diag ( beta, -1 );

  [ eigval, eigvec ] = eig ( tridiag );

  x = diag ( eigvec );
  [ x, i ] = sort ( x );
  w = 2 * eigval(1,i) .^ 2;

  quad = w * feval ( f, x );

  return
end


text:
function test_pull769

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_read_event

% the data was generated following the example at https://mne.tools/0.15/auto_tutorials/plot_creating_data_structures.html#tut-creating-data-structures
% and written to disk following https://mne.tools/0.15/auto_tutorials/plot_object_epochs.html
filename = dccnpath('/home/common/matlab/fieldtrip/data/test/pull769-epo.fif');

%%

hdr = ft_read_header(filename);
dat = ft_read_data(filename);
evt = ft_read_event(filename);

%%

cfg = [];
cfg.dataset = filename;
cfg.trialdef.eventvalue = 'smiling';
cfg = ft_definetrial(cfg);
data_smiling = ft_preprocessing(cfg);

assert(numel(data_smiling.time)==5);
assert(length(data_smiling.time{1})==200);

cfg = [];
cfg.dataset = filename;
cfg.trialdef.eventvalue = 'frowning';
cfg = ft_definetrial(cfg);
data_frowning = ft_preprocessing(cfg);

assert(numel(data_frowning.time)==5);
assert(length(data_frowning.time{1})==200);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Short test script for running algorithms related to the glottal source
%
% Please find the path management in the startup.m script in the root directory
% of this repository. Note that by starting matlab in the root directory, this
% script should automatically run. If it is not the case, you can also move to the
% root directory and run this script manually.
%
% License
%  This file is under the LGPL license,  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. This file 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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%

clear all;

% Settings
F0min = 80; % Minimum F0 set to 80 Hz
F0max = 500; % Maximum F0 set to 80 Hz
frame_shift = 10; % Frame shift in ms

% Load soundfile
[x,fs] = audioread(['howtos' filesep 'arctic_a0007.wav']);

% Check the speech signal polarity
polarity = polarity_reskew(x,fs);
x=polarity*x;

% Extract the pitch and voicing information
[srh_f0,srh_vuv,srh_vuvc,srh_time] = pitch_srh(x,fs,F0min,F0max,frame_shift);

% Extract the maximum voiced frequency
[max_voiced_freq] = maximum_voiced_frequency(x,fs,srh_f0.*srh_vuv,srh_time);

% Creaky probability estimation
warning off
try
    [creak_pp,creak_bin] = detect_creaky_voice(x,fs); % Detect creaky voice
    creak=interp1(creak_bin(:,2),creak_bin(:,1),1:length(x));
    creak(creak<0.5)=0; creak(creak>=0.5)=1;
    do_creak=1;
catch
    disp('Version or toolboxes do not support neural network object used in creaky voice detection. Creaky detection skipped.')
    creak=zeros(length(x),1);
    creak_pp=zeros(length(x),2);
    creak_pp(:,2)=1:length(x);
    do_creak=0;
end
warning on

% GCI estimation
sd_gci = gci_sedreams(x,fs,median(srh_f0),1);        % SEDREAMS
se_gci = se_vq(x,fs,median(srh_f0),creak);           % SE-VQ

%res = lpcresidual(x,25/1000*fs,5/1000*fs,fs/1000+2); % LP residual
res = maxlpresidual(x,fs,round(fs/1000)+2);

m = mdq(res,fs,se_gci); % Maxima dispersion quotient measurement

ps = peakslope(x,fs);   % peakSlope extraction

[gf_iaif,gfd_iaif] = iaif_ola(x,fs);    % Glottal flow (and derivative) by the IAIF method

dgf_cc = complex_cepstrum(x,fs,sd_gci,srh_f0,srh_vuv); % Glottal flow derivative by the complex cesptrum-based method

[NAQ,QOQ,H1H2,HRF,PSP] = get_vq_params(gf_iaif,gfd_iaif,fs,se_gci); % Estimate conventional glottal parameters

% Estimate the Rd parameter of the Liljencrants-Fant (LF) model
srh_f0(find(srh_f0==0)) = 100;
opt = sin_analysis();
opt.fharmonic  = true;
opt.use_ls     = false;
opt.debug = 0;
frames = sin_analysis(x, fs, [srh_time(:),srh_f0(:)], opt);
rds = rd_msp(frames, fs);

% Extract Cepstral Peak Prominence
CPPv = cpp( x, fs, 1, 'mean' );

% Plots
t=(0:length(x)-1)/fs;

fig(1) = subplot(311);
    plot(t,x, 'b');
    hold on
    plot(srh_time, srh_vuv, 'g');
    stem(sd_gci,ones(1,length(sd_gci))*-.1,'m');
    stem(se_gci,ones(1,length(se_gci))*-.1,'r');
    legend('Speech signal','Voicing (SRH)', 'GCI (SEDREAMS)','GCI (SE-VQ)');
    xlabel('Time [s]');
    ylabel('Amplitude');

fig(2) = subplot(312);
    plot(srh_time, srh_f0, 'r');
    legend('f0 (SRH)');
    xlabel('Time [s]');
    ylabel('Hz');

fig(3) = subplot(313);
    plot(srh_time, max_voiced_freq, 'g');
    legend('Maximum voiced frequency');
    xlabel('Time [s]');
    ylabel('Hz');

linkaxes(fig, 'x');
xlim([0 srh_time(end)]);
figure

fig(4) = subplot(311);
    plot(t,x, 'b');
    hold on
    plot(ps(:,1), ps(:,2),'--b');
    plot(m(:,1),m(:,2),'m');
    plot(rds(:,1), rds(:,2), 'r');
    legend('Speech signal','PeakSlope','MDQ','Rd');
    xlabel('Time [s]');

fig(5) = subplot(312);
    plot(t,x, 'b');
    hold on
    plot(t,norm(x)*gf_iaif./norm(gf_iaif), 'g');
    plot(creak_pp(:,2)/fs,creak_pp(:,1),'r')
    xlabel('Time [s]');
    ylabel('Amplitude');
    legend('Speech signal','Glottal flow (IAIF)','Creaky voice probability','Location','NorthWest');

fig(6) = subplot(313);
    plot(t,x, 'b');
    hold on
    plot(t,norm(x)*dgf_cc./norm(dgf_cc), 'r');
    xlabel('Time [s]');
    ylabel('Amplitude');
    legend('Speech signal','Glottal flow derivative (CC)');

linkaxes(fig, 'x');
xlim([0 srh_time(end)]);

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [Q, R] = stiefelgeneralized_retraction_MGS(B, X, U, t)
% Retraction for generalized Stiefel based on Modified Gram-Schmidt.
% When used just as a retraction, only the output Q is relevant.
% NB, Dec. 16, 2018.
    if ~exist('t', 'var') || isempty(t)
        A = X + U;   % t = 1 by default
    else
        A = X + t*U;
    end
    [n, p] = size(X);
    Q = zeros(n, p);
    R = zeros(p, p);
    for j = 1 : p
        v = A(:, j);
        R(j, j) = sqrt(v'*B*v);
        Q(:, j) = v / R(j, j);
        R(j, (j+1):p) = Q(:, j)' * B * A(:, (j+1):p);
        A(:, (j+1):p) = A(:, (j+1):p) - Q(:, j) * R(j, (j+1):p);
    end
end

author: NicolasBoumal

repo: manopt


text:
function result = GetGoogleSpreadsheet(DOCID)
% result = GetGoogleSpreadsheet(DOCID)
% Download a google spreadsheet as csv and import into a Matlab cell array.
%
% [DOCID] see the value after 'key=' in your spreadsheet's url
%           e.g. '0AmQ013fj5234gSXFAWLK1REgwRW02hsd3c'
%
% [result] cell array of the the values in the spreadsheet
%
% IMPORTANT: The spreadsheet must be shared with the "anyone with the link" option
%
% This has no error handling and has not been extensively tested.
% Please report issues on Matlab FX.
%
% DM, Jan 2013
%

loginURL = 'https://www.google.com';
csvURL = ['https://docs.google.com/spreadsheet/ccc?key=' DOCID '&output=csv&pref=2'];

%Step 1: go to google.com to collect some cookies
cookieManager = java.net.CookieManager([], java.net.CookiePolicy.ACCEPT_ALL);
java.net.CookieHandler.setDefault(cookieManager);
handler = sun.net.www.protocol.https.Handler;
connection = java.net.URL([],loginURL,handler).openConnection();
connection.getInputStream();

%Step 2: go to the spreadsheet export url and download the csv
connection2 = java.net.URL([],csvURL,handler).openConnection();
result = connection2.getInputStream();
result = char(readstream(result));

%Step 3: convert the csv to a cell array
result = parseCsv(result);

end

function data = parseCsv(data)
% splits data into individual lines
data = textscan(data,'%s','whitespace','\n');
data = data{1};
for ii=1:length(data)
   %for each line, split the string into its comma-delimited units
   %the '%q' format deals with the "quoting" convention appropriately.
   tmp = textscan(data{ii},'%q','delimiter',',');
   data(ii,1:length(tmp{1})) = tmp{1};
end

end

function out = readstream(inStream)
%READSTREAM Read all bytes from stream to uint8
%From: http://stackoverflow.com/a/1323535

import com.mathworks.mlwidgets.io.InterruptibleStreamCopier;
byteStream = java.io.ByteArrayOutputStream();
isc = InterruptibleStreamCopier.getInterruptibleStreamCopier();
isc.copyStream(inStream, byteStream);
inStream.close();
byteStream.close();
out = typecast(byteStream.toByteArray', 'uint8');

end


text:
function [err_p,elerr_p] = diffpost_bc(aez,fez,elerror,xy,ev,ebound);
%diffpost_bc postprocesses local Poisson error estimator
%   [err_p,elerr_p] = diffpost_bc(aez,fez,elerror,xy,ev,ebound);
%   input
%          aez       elementwise Poisson problem matrices
%          fez       elementwise rhs vectors
%          elerror   elementwise error estimate (without BC imposition)
%          xy        vertex coordinate vector
%          ev        element mapping matrix
%          ebound    element edge boundary matrix
%   output
%          err_p     global error estimate
%          elerr_p   elementwise error estimate
%   IFISS function: DJS; 4 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
      x=xy(:,1); y=xy(:,2);
      nel=length(ev(:,1));
      lev=[ev,ev(:,1)]; elerr_p=elerror;
%
% recompute contributions from elements with Dirichlet boundaries
      nbde=length(ebound(:,1));
      ebdy = zeros(nel,1);
      edge = zeros(nel,1);
% isolate boundary elements
      for el = 1:nbde
      ee = ebound(el,1);
      ebdy(ee) = ebdy(ee)+1; edge(ee)=ebound(el,2);
      end
%
% two edge elements
      k2=find(ebdy==2);
      nel2b=length(k2);
% loop over two edge elements
      for el = 1:nel2b
      el2e=k2(el);
      kk=find(ebound(:,1) == el2e);
      edges=ebound(kk,2);
% set up original matrix and RHS vector
	  ae=squeeze(aez(el2e,1:5,1:5));
      fe=fez(el2e,:)';
% set up local coordinates and impose interpolated error as Dirichlet bc
      xl=x(lev(el2e,:)); yl=y(lev(el2e,:));
      [bae,fe] = localbc_p(ae,fe,edges,xl,yl);
% solve local problem
      err=bae\fe;
      elerr_p(el2e,1) = err'*fe;
      end
% end of element loop
%
% one edge elements
      k1=find(ebdy==1);
      nel1b=length(k1);
% loop over one edge elements
      for el = 1:nel1b
      el1e=k1(el);
      kk=find(ebound(:,1) == el1e);
      edges=ebound(kk,2);
% set up original matrix and RHS vector
      fe=fez(el1e,:)';
	  ae=squeeze(aez(el1e,1:5,1:5));
% set up local coordinates and impose interpolated error as Dirichlet bc
      xl=x(lev(el1e,:)); yl=y(lev(el1e,:));
      [bae,fe] = localbc_p(ae,fe,edges,xl,yl);
% solve local problem
      err=bae\fe;
      elerr_p(el1e,1) = err'*fe;
      end
% end of element loop
%
      err_p = sqrt(sum(elerr_p));
      elerr_p = sqrt(elerr_p);
      fprintf('done\n')
      fprintf('estimated global error (in energy):  %10.6e\n',err_p)
 return


text:
function set_debug(debuglevel)
    % SET_DEBUG set the debug level.
    %
    %    SET_DEBUG is used in conjunction with PRINT_DEBUG to control the verboseness of MATLAB output.
    %    You use it to specify the debug level (an integer) at which to display a certain message.
    %
    %    EXAMPLE:
    %    Set the debug level to 1:
    %        debug.set_debug(1)
    %
    %    See also debug.get_debug, debug.print_debug, setpref.

    % AUTHOR: Ronni Grapenthin and Glenn Thompson, UAF-GI
    % $Date$
    % $Revision$
    % Modified by Glenn Thompson: 2010-03-10: To use multiple levels of debug (0 = none, 1 = some, 2 = higher)
    % GT 20150916: the setpref doesn't seem to work now print_debug is
    % using persistent variable. ANd now print_debug(level) without a text
    % message sets level. So replace this:

        %setpref('runmode', 'debug', debuglevel);

    % with
        debug.print_debug(debuglevel);

end

author: geoscience-community-codes

repo: GISMO


text:
function [permuteMatrix, nNonRepeats, nUniqueDirs, nUniqueMeasurements] = dtiBootGetPermMatrix(bvecs, bvals)
%
% [permuteMatrix, nNonRepeats, nUniqueDirs, nUniqueMeasurements] = dtiBootGetPermMatrix(bvecs, bvals)
%
%
% HISTORY:
% 2007.06.07 RFD wrote it

n = min(size(bvals,2),size(bvecs,2));
bv = [bvecs(:,1:n).*repmat(bvals(:,1:n),[3 1])];
nMeasurements = size(bv,2);
permuteMatrix = cell(nMeasurements,1);
for(ii=1:nMeasurements)
    dist1 = sqrt((bv(1,:)-bv(1,ii)).^2+(bv(2,:)-bv(2,ii)).^2+(bv(3,:)-bv(3,ii)).^2);
    dist2 = sqrt((bv(1,:)+bv(1,ii)).^2+(bv(2,:)+bv(2,ii)).^2+(bv(3,:)+bv(3,ii)).^2);
    permuteMatrix{ii} = unique([find(dist1<1e-3) find(dist2<1e-3)]);
end
numPerms = cellfun('length',permuteMatrix);
nNonRepeats = sum(numPerms==1);
m = zeros(nMeasurements,max(numPerms));
for(ii=1:nMeasurements)
    m(ii,1:length(permuteMatrix{ii})) = permuteMatrix{ii};
end
nUniqueMeasurements = size(unique(m,'rows'),1);
if(any(bvals<max(bvals)*0.01))
    nUniqueDirs = nUniqueMeasurements-1;
else
    nUniqueDirs = nUniqueMeasurements;
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ishappy, cutoff] = classicCheck(f, values, data, pref)
%CLASSICCHECK   Attempt to trim trailing Fourier coefficients in a TRIGTECH.
%   [ISHAPPY, CUTOFF] = CLASSICCHECK(F, VALUES, DATA) returns an estimated
%   location, the CUTOFF, at which the TRIGTECH F could be truncated to
%   maintain an accuracy of EPSLEVEL (see documentation below) relative to
%   DATA.VSCALE and DATA.HSCALE. ISHAPPY is TRUE if CUTOFF <
%   MIN(LENGTH(VALUES),2) or VSCALE(F) = 0 and FALSE otherwise.
%
%   [ISHAPPY, CUTOFF] = CLASSICCHECK(F, VALUES, DATA, PREF) allows additional
%   preferences to be passed. In particular, one can adjust the target accuracy
%   with PREF.CHEBFUNEPS.
%
%   CLASSICCHECK first queries HAPPINESSREQUIREMENTS to obtain TESTLENGTH and
%   EPSLEVEL (see documentation below). If |F.COEFFS(1:TESTLENGTH)|/VSCALE <
%   EPSLEVEL, then the representation defined by F.COEFFS is deemed happy. The
%   value returned in CUTOFF is essentially that from TESTLENGTH (although it
%   can be reduced if there are further COEFFS which fall below EPSLEVEL).
%
%   HAPPINESSREQUIREMENTS defines what it means for a TRIGTECH to be happy.
%   [TESTLENGTH, EPSLEVEL] = HAPPINESSREQUIREMENTS(VALUES, COEFFS, POINTS,
%   DATA, EPS) returns two scalars TESTLENGTH and EPSLEVEL.  POINTS
%   is the vector of points at which F was sampled to get the values in
%   VALUES.  EPS is the desired accuracy.  A TRIGTECH is deemed to be
%   'happy' if the coefficients COEFFS(END-TESTLENGTH+1:END) (recall that
%   COEFFS are stored in ascending order) are all below EPSLEVEL.
%   The default choice of the test length is:
%       TESTLENGTH = n,             for n = 1:2
%       TESTLENGTH = 3,             for n = 3:44
%       TESTLENGTH = round((n-1)/8) for n > 44
%
%   EPSLEVEL is essentially the maximum of:
%       * pref.chebfuneps
%       * eps*TESTLENGTH
%       * eps*condEst (where condEst is an estimate of the condition number
%                      based upon a finite difference approximation to the
%                      gradient of the function from F.VALUES.).
%   However, the final two estimated values can be no larger than 1e-4.
%
%   Note that the accuracy check implemented in this function is the (roughly)
%   same as that employed in Chebfun v4.x.
%
% See also STRICTCHECK, LOOSECHECK.

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

% Deal with special cases ------------------------------------------------------

% Determine n (the length of the input).
n = length(f);

% Assume we're not happy. (N'aww! :( )
ishappy = false;

% Grab some preferences:
if ( nargin == 1 )
    pref = f.techPref();
    epslevel = pref.chebfuneps;
elseif ( isnumeric(pref) )
    epslevel = pref;
else
    epslevel = pref.chebfuneps;
end

% Convert scalar epslevel/tolerance inputs into vectors.
if ( isscalar(epslevel) )
    epslevel = repmat(epslevel, 1, size(f.coeffs, 2));
end

% Deal with the trivial case:
if ( n < 2 ) % (Can't be simpler than a constant!)
    cutoff = n;
    return
end

% Check the vertical scale:
if ( max(data.vscale) == 0 )
    % This is the zero function, so we must be happy!
    ishappy = true;
    cutoff = 1;
    return
elseif ( any(isinf(data.vscale)) )
    % Inf located. No cutoff.
    cutoff = n;
    return
end

% If one column of f is the zero function, we will get into trouble further
% down when we take the absolute value of the coefficients relative to vscale
% and compute the relative condition number estimate in happinessCheck.  We
% replace zero vscales by eps to avoid division by zero.
data.vscale(data.vscale == 0) = eps;

% NaNs are not allowed.
if ( any(isnan(f.coeffs(:))) )
    error('CHEBFUN:FUN:classicCheck:NaNeval', ...
        'Function returned NaN when evaluated.')
end

% We require values. Get these before we alter the coeffs:
if ( isempty(values) )
    values = f.coeffs2vals(f.coeffs);
end

% Do the test on the vector formed by the sum of the absolute value of the
% positive and negative mode coefficients.
% [TODO] This reversal of coeffs can be removed but then classicCheck will need
% to be written carefully:
absCoeffs = abs(f.coeffs(end:-1:1,:));

% Need to handle odd/even cases separately.
isEven = ~mod(n, 2);
if ( isEven )
    % In this case the negative cofficients have an additional term
    % corresponding to the cos(N/2*x) coefficient.
    f.coeffs = [absCoeffs(n,:);absCoeffs(n-1:-1:n/2+1,:)+absCoeffs(1:n/2-1,:);absCoeffs(n/2,:)];
else
    f.coeffs = [absCoeffs(n:-1:(n+1)/2+1,:)+absCoeffs(1:(n+1)/2-1,:);absCoeffs((n+1)/2,:)];
end

n = size(f.coeffs, 1);

% Check for convergence and chop location --------------------------------------

% Absolute value of coefficients, relative to vscale: (max across columns)
ac = bsxfun(@rdivide, abs(f.coeffs), data.vscale);

% Happiness requirements:
[testLength, epslevel] = ...
    happinessRequirements(values, f.coeffs, f.points(), data, epslevel);

if ( all(max(ac(1:testLength, :)) < epslevel) ) % We have converged! Chop tail:
    % We must be happy.
    ishappy = true;

    % Find first row of coeffs with entry above epslevel:
    rowsWithLargeCoeffs = any(bsxfun(@ge, ac, epslevel), 2);
    Tloc = find(rowsWithLargeCoeffs, 1, 'first') - 1;

    % Check for the zero function!
    if ( isempty(Tloc) )
        cutoff = 1;
        return
    end

    % Compute the cumulative max of eps/4 and the tail entries:
    t = .25*eps*ones(1, size(ac, 2));
    ac = ac(1:Tloc, :);             % Restrict to coefficients of interest.
    for k = 1:size(ac, 1)           % Cumulative maximum.
        ind = ac(k,:) < t;
        ac(k, ind) = t(ind);
        ind = ac(k,:) >= t;
        t(ind) = ac(k,ind);
    end

    % Obtain an estimate for much accuracy we'd gain compared to reducing
    % length ("bang for buck"):
    bang = log(1e3*bsxfun(@rdivide, epslevel, ac));
    buck = n - (1:Tloc).';
    Tbpb = bsxfun(@rdivide, bang, buck);

    % Compute position at which to chop.  Keep greatest number of coefficients
    % demanded by any of the columns.
    [ignored, perColTchop] = max(Tbpb(3:Tloc, :));
    Tchop = min(perColTchop);

    % Cutoff value.
    cutoff = n - Tchop - 2;
    % We want to keep [c(-cutoff), ...,c(-1), c(0), c(1), ..., c(cutoff)],
    % so the number of coefficients that will be thrown away is actually 2*cutoff-1.
    cutoff = 2*cutoff-1;

else

    % We're unhappy. :(
    cutoff = n;

end

end

function [testLength, epslevel] = ...
    happinessRequirements(values, coeffs, x, data, epslevel) %#ok<INUSL>
%HAPPINESSREQUIREMENTS   Define what it means for a TRIGTECH to be happy.
%   See documentation above.

% Grab the size:
n = size(coeffs, 1);

% We will not allow the estimated rounding errors to be cruder than this value:
minPrec = 1e-4; % Worst case precision!

% Length of tail to test.
testLength = min(n, max(3, round((n-1)/8)));

% Look at length of tail to loosen tolerance:
tailErr = eps*testLength;
tailErr = min(tailErr, minPrec);

% Estimate the condition number of the input function by
% ||f(x+eps(x)) - f(x)||_inf / ||f||_inf ~~ (eps(hscale)/vscale)*f'.
dy = diff(values);
dx = diff(x)*ones(1, size(values, 2));
gradEst = max(abs(dy./dx));                       % Finite difference approx.
condEst = eps(data.hscale)./data.vscale.*gradEst; % Condition number estimate.
condEst = min(condEst, minPrec);

% Choose maximum between prescribed tolerance and estimated rounding errors:
epslevel = max(max(epslevel, condEst), tailErr);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun