text:
classdef AMICO_VERDICTPROSTATE

properties
    id, name                % id and name of the model
    max_dirs                % maximum number of directions to fit
    dIC                     %
    Rs                      %
    dEES                    %
    P                       %
    OUTPUT_names            % suffix of the output maps
    OUTPUT_descriptions     % description of the output maps
end

methods

    % =================================
    % Setup the parameters of the model
    % =================================
    function obj = AMICO_VERDICTPROSTATE()
        global CONFIG

        % set the parameters of the model
        obj.id        = 'VerdictProstate';
        obj.name      = 'VERDICT prostate';
        obj.max_dirs  = 0; % no need to estimate directions, it's an isotropic model
        obj.dIC       = 2.0 * 1E-3;
        obj.Rs        = linspace(0.01,20.1,20);
        obj.dEES      = 2.0 * 1E-3;
        obj.P         = 8.0 * 1E-3;
        obj.OUTPUT_names        = {'R', 'fIC', 'fEES', 'fVASC', 'Fobj'};
        obj.OUTPUT_descriptions = {'R', 'fIC', 'fEES', 'fVASC', 'Fobj'};

        % set the parameters to fit it
        CONFIG.OPTIMIZATION.SPAMS_param.mode    = 2;
        CONFIG.OPTIMIZATION.SPAMS_param.pos     = true;
        CONFIG.OPTIMIZATION.SPAMS_param.lambda  = 0;    % l1 regularization
        CONFIG.OPTIMIZATION.SPAMS_param.lambda2 = 1e-3; % l2 regularization
    end

    % ==================================================================
    % Generate high-resolution kernels and rotate them in harmonic space
    % ==================================================================
    function GenerateKernels( obj, ATOMS_path, schemeHR, AUX, idx_IN, idx_OUT )
        global CONFIG AMICO_data_path CAMINO_path

        % check if high-resolution scheme has been created
        schemeHrFilename = fullfile(ATOMS_path,'protocol_HR.scheme');
        if ~exist( schemeHrFilename, 'file' )
            error( '[AMICO_VERDICTPROSTATE.GenerateKernels] File "protocol_HR.scheme" not found in folder "%s"', ATOMS_path )
        end

        filenameHr = [tempname '.Bfloat'];
        progress = ProgressBar( numel(obj.Rs) + 2 );

        % IC compartment
        % ==============
        for R = obj.Rs
            % generate
            if exist( filenameHr, 'file' ), delete( filenameHr ); end
            CMD = sprintf( '%s/datasynth -synthmodel compartment 1 SPHEREGPD %E %E -schemefile %s -voxels 1 -outputfile %s 2> /dev/null', CAMINO_path, obj.dIC*1e-6, R*1e-6, schemeHrFilename, filenameHr );
            [status result] = system( CMD );
            if status>0
                disp(result)
                error( '[AMICO_VERDICTPROSTATE.GenerateKernels] Problems generating the signal with datasynth' );
            end

            % rotate and save
            fid = fopen( filenameHr, 'r', 'b' );
            signal = fread(fid,'float');
            fclose(fid);
            delete( filenameHr );
            lm = AMICO_RotateKernel( signal, AUX, idx_IN, idx_OUT, true );
            save( fullfile( ATOMS_path, sprintf('A_%03d.mat',progress.i) ), '-v6', 'lm' )
            progress.update();
        end

        % EES compartment
        % ===============
        % generate
        if exist( filenameHr, 'file' ), delete( filenameHr ); end
        CMD = sprintf( '%s/datasynth -synthmodel compartment 1 BALL %E -schemefile %s -voxels 1 -outputfile %s 2> /dev/null', CAMINO_path, obj.dEES*1e-6, schemeHrFilename, filenameHr );
        [status result] = system( CMD );
        if status>0
            disp(result)
            error( '[AMICO_VERDICTPROSTATE.GenerateKernels] problems generating the signal' );
        end

        % resample and save
        fid = fopen( filenameHr, 'r', 'b' );
        signal = fread(fid,'float');
        fclose(fid);
        delete( filenameHr );
        lm = AMICO_RotateKernel( signal, AUX, idx_IN, idx_OUT, true );
        save( fullfile( ATOMS_path, sprintf('A_%03d.mat',progress.i) ), '-v6', 'lm' )
        progress.update();

        % VASC compartment
        % ================
        % generate
        if exist( filenameHr, 'file' ), delete( filenameHr ); end
        CMD = sprintf( '%s/datasynth -synthmodel compartment 1 ASTROSTICKS %E -schemefile %s -voxels 1 -outputfile %s 2> /dev/null', CAMINO_path, obj.P*1e-6, schemeHrFilename, filenameHr );
        [status result] = system( CMD );
        if status>0
            disp(result)
            error( '[AMICO_VERDICTPROSTATE.GenerateKernels] problems generating the signal' );
        end

        % resample and save
        fid = fopen( filenameHr, 'r', 'b' );
        signal = fread(fid,'float');
        fclose(fid);
        delete( filenameHr );
        lm = AMICO_RotateKernel( signal, AUX, idx_IN, idx_OUT, true );
        save( fullfile( ATOMS_path, sprintf('A_%03d.mat',progress.i) ), '-v6', 'lm' )
        progress.update();

        progress.close();
end

    % ==============================================
    % Project kernels from harmonic to subject space
    % ==============================================
    function ResampleKernels( obj, ATOMS_path, idx_OUT, Ylm_OUT )
        global CONFIG AMICO_data_path KERNELS

        % Setup the KERNELS structure
        % ===========================
        nIC = numel(obj.Rs);

        KERNELS = {};
        KERNELS.model    = 'VERDICTPROSTATE';
        KERNELS.nS       = CONFIG.scheme.nS;
        KERNELS.nA       = nIC + 2; % number of atoms

        KERNELS.Aic      = zeros( [KERNELS.nS nIC], 'single' );
        KERNELS.Aic_R    = zeros( 1, nIC, 'single' );

        KERNELS.Aees      = zeros( [KERNELS.nS 1], 'single' );
        KERNELS.Aees_d    = NaN;

        KERNELS.Avasc     = zeros( [KERNELS.nS 1], 'single' );
        KERNELS.Avasc_P   = NaN;

        progress = ProgressBar( KERNELS.nA );

        % IC compartment
        % ==============
        for i = 1:nIC
            load( fullfile( ATOMS_path, sprintf('A_%03d.mat',progress.i) ), 'lm' );
            KERNELS.Aic(:,i) = AMICO_ResampleKernel( lm, idx_OUT, Ylm_OUT, true );
            KERNELS.Aic_R(i) = obj.Rs(i);
            progress.update();
        end

        % Precompute norms of coupled atoms (for the l1 minimization)
        A = double( KERNELS.Aic(CONFIG.scheme.dwi_idx,:) );
        KERNELS.Aic_norm = repmat( 1./sqrt( sum(A.^2) ), [size(A,1),1] );
        clear A

        % EES compartment
        % ===============
        load( fullfile( ATOMS_path, sprintf('A_%03d.mat',progress.i) ), 'lm' );
        KERNELS.Aees   = AMICO_ResampleKernel( lm, idx_OUT, Ylm_OUT, true );
        KERNELS.Aees_d = obj.dEES;
        progress.update();

        % VASC compartment
        % ================
        load( fullfile( ATOMS_path, sprintf('A_%03d.mat',progress.i) ), 'lm' );
        KERNELS.Avasc   = AMICO_ResampleKernel( lm, idx_OUT, Ylm_OUT, true );
        KERNELS.Avasc_P = obj.P;
        progress.update();

        progress.close();
    end

    % ===========================
    % Fit the model to each voxel
    % ===========================
    function [ MAPs ] = Fit( obj, y, i1, i2 )
        global CONFIG KERNELS

%         A = double( [ KERNELS.Aic(CONFIG.scheme.dwi_idx,:) KERNELS.Aees(CONFIG.scheme.dwi_idx) KERNELS.Avasc(CONFIG.scheme.dwi_idx) ] );
%         AA = [ ones(1,KERNELS.nA) ; A ];
%         yy = [ 1 ; y(CONFIG.scheme.dwi_idx) ];
        A = double( [ KERNELS.Aic KERNELS.Aees KERNELS.Avasc ] );
        AA = A;%[ ones(1,KERNELS.nA) ; A ];
        yy = y;%[ 1 ; y(CONFIG.scheme.dwi_idx) ];

        switch( 1 )

            case 1
            % [ normal fitting ]
            norms = repmat( 1./sqrt(sum(AA.^2)), [size(AA,1),1] );
            AA = AA .* norms;
            x = full( mexLasso( yy, AA, CONFIG.OPTIMIZATION.SPAMS_param ) );
            x = x .* norms(1,:)';

            case 2
            params = CONFIG.OPTIMIZATION.SPAMS_param;
            params.loss = 'square';
            params.regul = 'group-lasso-l2';
            params.groups = int32([ repmat(1,1,numel(KERNELS.Aic_R)) 2 3 ]);  % all the groups are of size 2
            params.intercept = false;
            x = full( mexFistaFlat( yy, AA, zeros(size(A,2),1), params ) );

            case 3
            % [NODDI style] estimate and remove EES and VASC contributions
            y = y(CONFIG.scheme.dwi_idx);
            x = lsqnonneg( AA, yy, CONFIG.OPTIMIZATION.LS_param );
            y = y - x(end)*A(:,end);
            y = y - x(end-1)*A(:,end-1);

            % find sparse support for remaining signal
            An = A(:,1:size(KERNELS.Aic,2)) .* KERNELS.Aic_norm;
            x = full( mexLasso( y, An, CONFIG.OPTIMIZATION.SPAMS_param ) );

            % debias coefficients
            idx = [ x>0 ; true ; true ];
            x(idx) = lsqnonneg( AA(:,idx), yy, CONFIG.OPTIMIZATION.LS_param );
        end

        % compute MAPS
        xIC = x( 1:end-2 );
        fIC = sum( xIC );
        MAPs(1) = KERNELS.Aic_R * xIC / ( fIC + eps );           % cell radius
        MAPs(2) = fIC;                                           % fIC
        MAPs(3) = x( end-1 );                                    % fEES
        MAPs(4) = x( end );                                      % fVASC

        y_predicted = A*x;
        MAPs(5) = norm(y-y_predicted) / norm(y);
%         sigma = 0;
%         MAPs(5) = sum( (y - sqrt((y_predicted).^2+sigma^2) ) ).^2;
    end

end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% Build and initialize the computational graph for regression based speech
% enhancement/dereverberation
%
function [layer, para] = Build_EnhanceNet_Masking(Data_tr, para)
para.output = 'tmp';

layer = genNetworkDereverb_Masking(para.topology);     % generate the network graph
para.preprocessing{1} = {};                     % optional preprocessing for each data stream
para.preprocessing{2} = {};
para.cost_func.layer_idx = length(layer);       % specify which layers are cost function layers
para.cost_func.layer_weight = [1];              % set the weights of each cost function layer
para.IO.inputFeature = [1 1];
para.IO.isTensor = [1 1];
para = ParseOptions2(para);

% generating the scaling factor for the input, as we will need to use a
% small constant in the logarithm. We need to make sure that the power of
% speech are larger than this constant most of the time.
scale = 1e4;        % we hard code the scale to be a constant so that all network will use the same number
scale = scale/2^16; % note that we are using int16 to store waveform samples, so need to scale down
if para.topology.useWav
    layer = InitWavScaleLayer(layer, scale);
end

if strcmpi(para.topology.RegressionNetType, 'DNN')    % if use DNN, splice the frames
    idx = ReturnLayerIdxByName(layer, 'splice');
else
    idx = ReturnLayerIdxByName(layer, 'delta'); % if use LSTM, use dynamic features
end
fft_net_length = idx(1);
fft_net = layer(1:fft_net_length);
paraTmp = para;
paraTmp.out_layer_idx = fft_net_length;
paraTmp.IO = RemoveIOStream(paraTmp.IO, [2]);
paraTmp.IO.nStream = 1;
paraTmp.IO.inputFeature = paraTmp.IO.DynamicDistortion.inputFeature(1);
if isfield(paraTmp.IO, 'fileReader'); paraTmp.IO = rmfield(paraTmp.IO, 'fileReader'); end
paraTmp.IO.fileReader(1)= paraTmp.IO.DynamicDistortion.fileReader(1);
fprintf('Generate global MVN weights for mask subnet - %s\n', datestr(now));
[layer{fft_net_length+1}.W, layer{fft_net_length+1}.b] = computeGlobalCMVN(Data_tr(1), 100, paraTmp, fft_net);
VerifyPreprocessingTree(layer(1:fft_net_length+1), Data_tr(1), paraTmp, 100);

% set weight of static, velocity, and accelration features in the MSE cost
% function.
delta_idx = ReturnLayerIdxByName(layer, 'delta');
weight_idx = delta_idx+1;
for i=1:length(weight_idx)
    if isfield(layer{weight_idx(i)}, 'W') && numel(layer{weight_idx(i)}.W) == prod(layer{weight_idx(i)}.dim)
        continue;
    end
    layer{weight_idx(i)}.W = diag([ones(para.topology.nFreqBin,1)*para.topology.MSECostWeightSDA(1); ...
        ones(para.topology.nFreqBin,1)*para.topology.MSECostWeightSDA(2); ...
        ones(para.topology.nFreqBin,1)*para.topology.MSECostWeightSDA(3)]);
    layer{weight_idx(i)}.b = zeros(para.topology.nFreqBin*3,1);
end

end

author: singaxiong

repo: SignalGraph


text:
function gX =  whiteXwhiteKernGradX(whiteKern1, whiteKern2, X, X2)

% WHITEXWHITEKERNGRADX

% KERN

if nargin < 3,
    X2 = X;
else
    U = X;
    X = X2;
    X2 = U;
end

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

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function p = SquareSubplot(n,i)

%SquareSubplot - Layout subplots in a square arrangement.
%
%
%  USAGE
%
%    p = SquareSubplot(n,i)
%
%    n              total number of subplots
%    i              current subplot

% Copyright (C) 2008-2011 by Michaël Zugaro
%
% 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.

if nargin < 2,
  error('Incorrect number of parameters (type ''help SquareSubplot'' for details).');
end

M = round(sqrt(n));
N = ceil(n/M);
p = subplot(M,N,i);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [mwcn,wpixc, globmax,edgedist]=localminW(w,nx,ny,psf)
% [mwcn,wpixc, globmax]=analyzeWconfPSF(w,nx,ny,psf)
% Finds local maxima in the results w (nx*ny,ncomp) convoluted with psf
% wpixc - w convolved with estimated PSF.
% mwcn = image of local maxima (non zero pixels) with values of hte convoluted image, nornalised to the maximum of each column.
% globmax = maximum of each column of W convolved with PSF. - PSF and W are L1 normalised - teh maximum can be used as a "quality measure."
%
% To get number of local maxima with relative strength < .5: sum(mwcn>.5)
ncomp=size(w,2);
autoconvmax=max(max(conv2(psf,psf,'same'))); % maximum of the "autoconvolution" - use to normalise the convolution with w.
wpix=reshape(w,nx, ny, ncomp); % each frame normalized to 1
wpixc=convstack(wpix,psf,'same');
edgedist=maxtoedgedist(wpixc);
wc=reshape(wpixc,nx*ny,ncomp);
mpix=maximastack(wpixc); % binary image of local maxima locations in each frame
m=reshape(mpix,nx*ny,ncomp);
mwc=m.*wc; % pixels indicates locations and value of hte local maxima
globmax=max(mwc)/autoconvmax; % The value of the global maximum for each column of W normalised with respect to the 'autoconvolution' of the PSF.
mwcn=normcMax(mwc); % normalised to the maximum of each column.

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function TR=SubdivideSphericalMesh(TR,k)
% Subdivide triangular mesh representing the surface of the unit sphere
% k times using triangular quadrisection (see function TriQuad for more
% info). The newly inserted vertices are re-projected onto the unit sphere
% after every iteration.
%
% INPUT ARGUMENTS:
%   - TR   : input mesh. TR can must be specified as a TriRep object.
%   - k    : desired number of subdivisions. k=0 is default.
%
% OUTPUT:
%   - TR  : subdivided mesh. Same format as input.
%
% AUTHOR: Anton Semechko (a.semechko@gmail.com)
% DATE: June.2012
%

if nargin<2 || isempty(k), k=0; end

if k==0, return; end
for i=1:k

    % Subdivide the mesh
    TR=TriQuad(TR);

    % Project the points onto the surface of the unit sphere
    x=TR.X;
    x_L2=sqrt(sum(x.^2,2));
    x=bsxfun(@rdivide,x,x_L2);
    TR=TriRep(TR.Triangulation,x);

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function tanh_prime = tanh_prime(a)
    %a = tanh(input);
    tanh_prime = (1-a.^2);
end

author: jacoxu


text:
function value = p31_exact ( dim_num )

%*****************************************************************************80
%
%% P31_EXACT returns the exact integral for problem 31.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Kenneth Hanson,
%    Quasi-Monte Carlo: halftoning in high dimensions?
%    in Computatinal Imaging,
%    Edited by CA Bouman and RL Stevenson,
%    Proceedings SPIE,
%    Volume 5016, 2003, pages 161-172.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Output, real VALUE, the exact value of the integral.
%

%
%  Get the limits of integration.
%
  [ a, b ] = p31_lim ( dim_num );
%
%  Get the coefficient vector C.
%
  c = [];
  c = p31_r8vec ( 'G', 'C', dim_num, c );
%
%  Get the location of Z.
%
  z = [];
  z = p31_r8vec ( 'G', 'Z', dim_num, z );
%
%  The value of the DIM_NUM dimensional integral is separable
%  into the product of integrals over each dimension.
%
%  Each of these 1 dimensional integrals, in turn, is
%  easily computed, depending on where Z(I) lies with
%  respect to the limits of integration A(I) and B(I).
%
  value = 1.0;

  for i = 1 : dim_num
%
%  Z < A < B
%
%  | X - Z | = X - Z from A to B.
%
    if ( z(i) < a(i) )

      value = value * ...
      ( exp ( - c(i) * ( a(i) - z(i) ) ) ...
      - exp ( - c(i) * ( b(i) - z(i) ) ) ) / c(i);
%
%  A < Z < B
%
%  | X - Z | = Z - X from B to Z,
%            = X - Z from      Z to A.
%
    elseif ( z(i) < b(i) )

      value = value * ( 2.0 ...
          - exp ( - c(i) * ( z(i) - a(i) ) ) ...
          - exp ( - c(i) * ( b(i) - z(i) ) ) ) / c(i);
%
%  A < B < Z
%
%  | X - Z | = Z - X from A to B.
%
    else

      value = value * ...
      ( exp ( - c(i) * ( z(i) - b(i) ) ) ...
      - exp ( - c(i) * ( z(i) - a(i) ) ) ) / c(i);

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = pyramid_grid_size ( n )

%*****************************************************************************80
%
%% PYRAMID_GRID_SIZE sizes a pyramid grid.
%
%  Discussion:
%
%    0:  x
%
%    1:  x  x
%        x  x
%
%    2:  x  x  x
%        x  x  x
%        x  x  x
%
%    3:  x  x  x  x
%        x  x  x  x
%        x  x  x  x
%        x  x  x  x
%
%    N  Size
%
%    0     1
%    1     5 = 1 + 4
%    2    14 = 1 + 4 + 9
%    3    30 = 1 + 4 + 9 + 16
%    4    55 = 1 + 4 + 9 + 16 + 25
%    5    91 = 1 + 4 + 9 + 16 + 25 + 36
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of subintervals.
%
%    Output, integer VALUE, the number of
%    nodes in the grid of size N.
%
  np1 = n + 1;

  value = ( np1 * ( np1 + 1 ) * ( 2 * np1 + 1 ) ) / 6;

  return
end


text:
classdef NDK
    % NDK works with NDK files can convert it to a ZmapCatalog
    %
    % from : https://www.ldeo.columbia.edu/~gcmt/projects/CMT/catalog/allorder.ndk_explained
    %
    % CMT files for download at: http://www.globalcmt.org/CMTfiles.html
    %
    % CMT reference:
    %   Dziewonski, A. M., T.-A. Chou and J. H. Woodhouse, Determination of earthquake source parameters from waveform data for studies of global and regional seismicity, J. Geophys. Res., 86, 2825-2852, 1981. doi:10.1029/JB086iB04p02825
    %   Ekström, G., M. Nettles, and A. M. Dziewonski, The global CMT project 2004-2010: Centroid-moment tensors for 13,017 earthquakes, Phys. Earth Planet. Inter., 200-201, 1-9, 2012. doi:10.1016/j.pepi.2012.04.002
    properties
        allNDKs table

        % TypeOfSourceInvertedFor is
        %"CMT: 0" - general moment tensor;
        %"CMT: 1" - moment tensor with constraint of zero trace (standard);
        %"CMT: 2" - double-couple source.

        % from line 3 of NDK format : CMT info (2)
        % Centroid_Time is  wrt reference time

        % CentroidDepthType is one of FREE, FIX, BODY

        % from line 4 of NDK format : CMT info (3)
        % ExponentForAllMomentValues most all moment measurements need to be multipleid by exponent
        %Mrr % r is up
        %Mtt % t is south
        %Mpp % p is east
    end
    properties(Constant)
        base_directory = 'https://www.ldeo.columbia.edu/~gcmt/projects/CMT/catalog/';
        % different dates are stored in differing locations.
        filelocations = {... [startdate>=, <enddate],location relative to base.
            datetime([1962 1976], 1, 1), {'PRE1976/deep_1962-1976.ndk','PRE1976/intdep_1962-1976.ndk'};
            datetime([1976 2017], 1, 1), 'jan76_dec17.ndk';...
            datetime([2017 year(datetime)], 1, 1), 'NEW_MONTHLY/{yyyy}';...
            [datetime(year(datetime), 1, 1), datetime], 'NEW_QUICK/E{yyyymmddHHMM}.ndk'...
            }

    end

    methods
        function obj = NDK(ndktext)
            %
            %
            % parsing based on documentation found at :
            % https://www.ldeo.columbia.edu/~gcmt/projects/CMT/catalog/allorder.ndk_explained

            if nargin==0 || isempty(ndktext)
                return
            end
            if ischar(ndktext) && isvector(ndktext)
                ndktext=splitlines(strip(ndktext));
            end
            obj=obj.fillLine1(strjoin(ndktext(1:5:end),newline));
            obj=obj.fillLine2(strjoin(ndktext(2:5:end),newline));
            obj=obj.fillLine3(strjoin(ndktext(3:5:end),newline));
            obj=obj.fillLine4(strjoin(ndktext(4:5:end),newline));
            obj=obj.fillLine5(strjoin(ndktext(5:5:end),newline));

        end
        function obj = fillLine1(obj,s)
            % make sure reference catalog is not empty, otherwise textscan will choke
            noCat = strfind(s,[newline '    ']);
            s(noCat+1)='-';
            s(noCat+2)='-';
            s(noCat+3)='-';

            % make sure location is not empyt, otherwise textscan will choke
            sidexes = 57:81:numel(s);
            no_location= s(sidexes)==' ';
            s(sidexes(no_location))='-';
            % NDK's exist with 60 seconds. datetime chokes on this. Fix it.
            sidexes = 23:81:numel(s);
            add_a_minute = s(sidexes)=='6';
            s(sidexes(add_a_minute))='0';

            [c,POS]=textscan(s,'%4C %{uuuu/MM/DD}D %{HH:mm:ss.s}D %f %f %f %f %f %24c');
            obj.allNDKs.HypocenterReferenceCatalog=c{1};
            obj.allNDKs.ReferenceEventDateTime =c{2} + timeofday(c{3}) + minutes(add_a_minute(:));
            obj.allNDKs.Latitude = c{4};
            obj.allNDKs.Longitude=c{5};
            obj.allNDKs.Depth=c{6};
            obj.allNDKs.ReportedMagnitudes = [c{7},c{8}];
            obj.allNDKs.GeographicalLocation=categorical(string(c{9}));
        end

        function obj = fillLine2(obj,s)

            c = textscan(s, '%16s B:%d %d %d S:%d %d %d M:%d %d %d CMT:%C %5C:%f');

            % c is now a 1x13 cell, each cell is a field, and contains an array with one item per NDK
            obj.allNDKs.CMTEventName=c{1};
            obj.allNDKs.LPBodyWaves_NumStationsUsed=c{2};
            obj.allNDKs.LPBodyWaves_NumComponentsUsed=c{3};
            obj.allNDKs.LPBodyWaves_ShortestPeriodUsed=c{4};
            obj.allNDKs.IntermedPSurfaceWaves_NumStationsUsed=c{5};
            obj.allNDKs.IntermedPSurfaceWaves__NumComponentsUsed=c{6};
            obj.allNDKs.IntermedPSurfaceWaves_ShortestPeriodUsed=c{7};
            obj.allNDKs.LPMantleWaves_NumStationsUsed=c{8};
            obj.allNDKs.LPMantleWaves__NumComponentsUsed=c{9};
            obj.allNDKs.LPMantleWaves_ShortestPeriodUsed=c{10};
            obj.allNDKs.TypeOfSourceInvertedFor=c{11};
            obj.allNDKs.MomentRateFunctionType=c{12};
            obj.allNDKs.MomentRateFunctionDuration=c{13};
        end

        function obj = fillLine3(obj,s)

            [c, POS] = textscan(s, 'CENTROID: %f %f %f %f %f %f %f %f %C %c-%14c');
            obj.allNDKs.Centroid_Time =c{1};% wrt reference time
            obj.allNDKs.Centroid_TimeStderr = c{2};
            obj.allNDKs.Centroid_Lat = c{3};
            obj.allNDKs.Centroid_LatStderr = c{4};
            obj.allNDKs.Centroid_Lon = c{5};
            obj.allNDKs.Centroid_LonStderr = c{6};
            obj.allNDKs.Centroid_Depth = c{7};
            obj.allNDKs.Centroid_DepthStderr = c{8};
            obj.allNDKs.CentroidDepthType = c{9}; % one of FREE, FIX, BODY
            obj.allNDKs.TimeStampPrefix = c{10};
            obj.allNDKs.TimeStamp = c{11};
        end

        %from line 3 of NDK format : CMT info (2)

        function obj = fillLine4(obj,s)

            [c, POS] = textscan(s, '%d %f %f %f %f %f %f %f %f %f %f %f %f');
            % from line 4 of NDK format : CMT info (3)
            obj.allNDKs.ExponentForAllMomentValues=c{1};
            obj.allNDKs.Mrr=c{2}; % r is up
            obj.allNDKs.MrrStderr=c{3};
            obj.allNDKs.Mtt=c{4}; % t is south
            obj.allNDKs.MttStderr=c{5};
            obj.allNDKs.Mpp=c{6}; % p is east
            obj.allNDKs.MppStderr=c{7};
            obj.allNDKs.Mrt=c{8};
            obj.allNDKs.MrtStderr=c{9};
            obj.allNDKs.Mrp=c{10};
            obj.allNDKs.MrpStderr=c{11};
            obj.allNDKs.Mtp=c{12};
            obj.allNDKs.MtpStderr=c{13};
        end

        function obj = fillLine5(obj,s)
            [c,POS] = textscan(s, '%3C %f %d %d %f %d %d %f %d %d %f %d %d %d %d %d %d');

            obj.allNDKs.VersionCode=c{1};
            obj.allNDKs.MTEigenvalue=[c{2} c{3} c{4}];
            obj.allNDKs.MTPlunge=[c{5} c{6} c{7}];
            obj.allNDKs.MTAzimuth=[c{8} c{9} c{10}];
            obj.allNDKs.ScalarMoment=c{11}; % needs to be multiplied by exponent
            obj.allNDKs.Strike_NodalPlane1=c{12};
            obj.allNDKs.Dip_NodalPlane1=c{13};
            obj.allNDKs.Rake_NodalPlane1=c{14};
            obj.allNDKs.Strike_NodalPlane2=c{15};
            obj.allNDKs.Dip_NodalPlane2=c{16};
            obj.allNDKs.Rake_NodalPlane2=c{17};
        end

        function c = toZmapCatalog(obj)
            c = ZmapCatalog;
            c.Name          = 'ndk';
            c.Longitude     = obj.allNDKs.Longitude;
            c.Latitude      = obj.allNDKs.Latitude;
            c.Depth         = obj.allNDKs.Depth;
            c.Date          = obj.allNDKs.ReferenceEventDateTime;
            c.Magnitude     = obj.allNDKs.ReportedMagnitudes(:,1);
            c.MagnitudeType = repmat(categorical({''}),size(c.Magnitude));
            mt      = obj.allNDKs(:,{'Mrr','Mtt','Mpp','Mrt','Mrp','Mtp'});
            mt{:,:} = mt{:,:}.* 10.^double(obj.allNDKs.ExponentForAllMomentValues);
            mt.Properties.VariableNames = {'mrr', 'mtt', 'mff', 'mrt', 'mrf', 'mtf'};
            mtAddon                 = MomentTensorAddon(mt);
            mtAddon.Dip             = double(obj.allNDKs.Dip_NodalPlane1);
            mtAddon.DipDirection    = mod(double(obj.allNDKs.Strike_NodalPlane1)+ 90, 360 );
            mtAddon.Rake            = double(obj.allNDKs.Rake_NodalPlane1);
            c.setAddon(mtAddon);
        end
    end

    methods(Static)

        function NDKS = read(f)
            tx=fileread(f);
            NDKS = NDK(splitlines(strip(tx)));
        end

        function NDKS = download(varargin)
            p=inputParser();
            p.addRequired('starttime');
            p.addRequired('endtime');
            p.addParameter('minmagnitude',-inf);
            p.addParameter('maxmagnitude',inf);
            p.addParameter('minlatitude',-90);
            p.addParameter('maxlatitude',90);
            p.addParameter('minlongitude',-180);
            p.addParameter('maxlongitude',180);
            p.addParameter('mindepth',inf);
            p.addParameter('maxdepth',inf);
            % TODO finish this

        end
        function tb = getFileLocations(starttime, endtime)
            % unpack theoretical file locations into real ones.
            floc_cells = NDK.filelocations;
            fls = NDK.filelocations{:,1}(1);
            fle = NDK.filelocations{:,1}(2);
            ignoreMe =   starttime > fle || endtime <= fls ;
            floc_cells(ignoreMe,:) = [];
            tb = table('Size', [1,3],...
                'VariableTypes', {'datetime',   'datetime', 'string'},...
                'VariableNames', {'starttime',  'endtime',  'url'});
            % TODO finish this

        end

    end
end

author: CelsoReyes


text:
function representations = batch_flipped_chol_omp(data_matrix, nk, threshold, quiet)
    if nargin < 3
        threshold = 1e-3;
    end
    if nargin < 4
        quiet = false;
    end
    % Number of data vectors
    ns = size(data_matrix, 2);
    % Computes sparse representations of the data vectors
    data_matrix = spx.norm.normalize_l2(data_matrix);
    % support set for each vector [one row for each vector]
    support_sets = ones(ns, nk);
    % termination vector
    % contains the number of iterations in each a particular vector was solved.
    % initially it is set to K
    % the moment residual norm reduces below the threshold
    % we mark it as finished and store the number of iterations
    % in this vector
    termination_vector = nk * ones(ns, 1);
    % gram matrix
    gram_matrix = data_matrix' * data_matrix;
    % initialize the correlation matrix
    correlation_matrix = abs(gram_matrix);
    prev_delta = zeros(1, ns);
    % Initialization of residual norms
    res_norm_sqr = ones(1, ns);
    % the Cholesky factors
    L = zeros(ns, nk, nk);
    % Initialize it
    L(:, 1,1) = 1;
    for iter=1:nk
        if ~quiet
            fprintf('.');
        end
        k2 = iter-1;
        % set all the diagonal entries (inner product with self) to zero.
        correlation_matrix(1:ns+1:end) = 0;
        % the inner product of a residual with each atom is stored along a column
        % we need to take maximum of each column to identity the best matching atom
        [~, best_match_indices] = max(correlation_matrix, [], 1);
        % we fill in the support set
        support_sets(:, iter) = best_match_indices;
        if (iter ~= nk)
            % we need to compute residuals except for the last iteration.
            % iterate for each data point
            for s=1:ns
                if termination_vector(s) == nk
                    % pick support for this vector
                    support_set = support_sets(s, 1:iter);
                    lambda = support_set(end);
                    if iter > 1
                        b = gram_matrix(support_set(1:k2), lambda);
                        opts = struct();
                        opts.LT= true;
                        LL = reshape(L(s, 1:k2, 1:k2), [k2, k2]);
                        w = linsolve(LL, b, opts);
                        L(s, iter, 1:iter) = [w' sqrt(1 - w' * w)];
                    end
                    % now solve for the coefficients using LL' z = p(:,Gamma)
                    h0 = gram_matrix(:,s);
                    h0lamba = h0(support_set);
                    LL = reshape(L(s, 1:iter, 1:iter), [iter, iter]);
                    opts = struct();
                    opts.LT = true;
                    z1 = linsolve(LL, h0lamba, opts);
                    opts = struct();
                    opts.LT = true;
                    opts.TRANSA = true;
                    z = linsolve(LL, z1, opts);
                    % pick gram submatrix
                    submatrix = gram_matrix(:, support_set);
                    beta = submatrix * z;
                    h = h0 - beta;
                    % put the new correlations into correlation matrix
                    correlation_matrix(:, s) = abs(h);
                    delta = z' * beta(support_set);
                    res_norm_sqr(s) = res_norm_sqr(s) - delta + prev_delta(s);
                    prev_delta(s) = delta;
                    if res_norm_sqr(s) < threshold
                        % The processing of this vector is complete
                        termination_vector(s) = iter;
                    end
                end
            end
        end
    end
    % disp(termination_vector');
    % final computation of coefficients
    % the values of coefficients corresponding to support sets
    coefficients_matrix = zeros(ns, nk);
    % Creation of sparse representation matrix
    for s=1:ns
        % number of iterations for this vector
        k = termination_vector(s);
        support_set  = support_sets(s, 1:k);
        x = data_matrix(:, s);
        % pick atoms from the unnormalized data matrix
        submatrix = data_matrix(:, support_set);
        % solve the least squares problem
        coeff = submatrix \ x;
        % if (s == 1)
        %     disp(submatrix(1:10, :));
        %     disp(x(1:10));
        %     disp(coeff);
        % end
        % put the coefficients back into coefficients matrix
        coefficients_matrix(s, 1:k) = coeff';
    end

    % each column varies from 1 to ns
    % each row is just a repetition
    column_indices  = repmat((1:ns)', 1, nk);

    % the support set, column_indices and coefficients_matrix are combined to form
    % the representation matrix at the end of function as follows
    % the support set contains row number
    % the column_indices contains the column number
    % the values contains the value of the non-zero entry
    % ns, ns is the size of the sparse matrix.
    representations = sparse( support_sets(:), column_indices(:), coefficients_matrix(:), ns, ns);
    iterations = termination_vector;
    if ~quiet
        fprintf('\n');
    end
end

author: indigits

repo: sparse-plex


data = [ 10 37 5 6 6];
explode = [ 0 1 0 0 0];
pie(data,explode);
title('\bfExample of a Pie Plot');
legend('One','Two','Three','Four','Five');


text:
%  Figure 6.65      Feedback Control of Dynamic Systems, 5e
%                   Franklin, Powell, Emami
%

clear all
close all;

k=10;
num=k;
den=[1 1 0];
w=logspace(-3,1,100);
[m,p]=bode(num,den,w);
num=conv(num,[10 1]);
den=conv(den,[100 1]);
[mc,pc]=bode(num,den,w);
subplot(2,1,1)
w2=[.001 10];
m2=[1 1];
loglog(w,m,w,mc,w2,m2,'r');
grid;
%xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 6.65 Frequency response of lag-compensation design: magnitude.');
subplot(2,1,2)
semilogx(w,p,w,pc);
grid;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
title('Fig. 6.65 Frequency response of lag compensation design: phase.');
%return;
%numcl=[1 0.1];
%dencl=[1 1.01 1.01 0.1];
%t=0:.01:20;
%y=step(numcl,dencl,t);
%plot(t,y);
%grid;
%xlabel('sec');
%ylabel('y');

author: Sable

repo: mcbench-benchmarks


text:
function [meanMap,stdMap,TcircSq,p]=computeStatsOnComplexVectors(cVectors)
% [meanMap,stdMap]=computeStatsOnComplexVectors(cVectors)
% PURPOSE : Computes the mean and std along the 3rd dimension
%   on a 3D or 4D set of complex vectors
%   produced by getSingleCycleVectors or getSingleCycleVectorsMultipleScans
%   It's almost redundant except that if you pass in a 4D array, it
%   reshapes to a 3D array before computing the std and mean along the 3rd
%   dimension.

s=size(cVectors);
if (length(s)==4)
    cVectors=reshape(cVectors,[s(1),s(2),s(3)*s(4)]);
end
s=size(cVectors);
if (length(s)~=3)
    error('S must be a 3 or 4D array');
end

meanMap=mean(cVectors,3);
stdMap=std(cVectors,[],3);

zeroMean=cVectors-repmat(meanMap,[1,1,s(3)]);

rStd=sum(real(zeroMean).^2/s(3),3);
iStd=sum(imag(zeroMean).^2/s(3),3);

Vindiv=(rStd+iStd)/(2*(s(3)-1));

Vgroup=(s(3)/2)*(abs(meanMap).^2);
TcircSq=(1/s(3))*(Vgroup./Vindiv);
p=FTest(2,(2*s(3)-2),s(3)*TcircSq);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% tutorial2_solutions
%   This script contains the solutions for Tutorial 2, see Tutorial 2 in
%   "RAVEN tutorials.docx" for more details.

%Import the Excel model into a RAVEN model structure
smallModel=importExcelModel('small.xlsx');

%This solves the linear programming problem.
%NOTE: if sol.f is equal to zero then the problem is not solvable. Ensure
%that uptake and excretion of all necessary stuff are added to the model
%and run the solveLP again.
sol=solveLP(smallModel);

%Print the resulting exchange fluxes
printFluxes(smallModel,sol.x,true);

%Print all fluxes
printFluxes(smallModel,sol.x,false,10^-5,[],'%rxnID (%rxnName):\n\t%eqn\n\t%flux\n');

author: SysBioChalmers

repo: RAVEN


text:
 function out = subsref(ob, args)
%function out = subsref(ob, args)
% handle subscript references like ob.ref or ob(ref,:)
% Copyright 2002-2-20, Jeff Fessler, University of Michigan

out = fatrix2_subsref1(ob, args(1)); % handle first subscript

if numel(args) > 1 % handle multiple subscripts, e.g., ob.field()
	out = subsref(out, args(2:end));
end

% fatrix2_subsref1()
function out = fatrix2_subsref1(ob, arg)

switch(arg.type)

case '.' % ob.?
	out = fatrix2_subsref_dot(struct(ob), arg);
	return

case '()' % ob(?)
	subs = arg.subs;

	% A(:,:) or A(:,j)
	if length(subs) == 2 && streq(subs{1}, ':')
		out = fatrix2_subsref_colon(ob, subs{2});

	% A(i,:)
	elseif length(subs) == 2 && streq(subs{2}, ':')
		out = fatrix2_subsref_colon(ob', subs{1})'; % trick: (A'(:,i))'

	% A(i,j)
	elseif length(subs) == 2
		out = fatrix2_subsref_colon(ob, subs{2}); % A(:,j) first
		if numel(subs{2}) > 1
			out = full(out); % fatrix2 to array
		end
		out = out(subs{1},:); % A(i,j) using normal array indexing

	% A(:)
	elseif length(subs) == 1 && streq(subs{1}, ':')
		out = col(fatrix2_subsref_colon(ob, ':')); % col(A(:,:))

	% A(k)
	elseif length(subs) == 1 && isnumeric(subs{1})
		[ii, jj] = ind2sub(size(ob), subs{1});
		out = zeros(numel(jj),1);
		for ll=1:numel(jj)
			tmp = fatrix2_subsref_colon(ob, jj(ll)); % A(:,j)
			out(ll) = tmp(ii(ll)); % A(i,j)
		end

	else
		printm('A subsref of type () called with these args:')
		pr subs
		error('That subsref type is unknown.  Do you mean {1}?')
	end

case '{}' % ob{?} (for ordered-subsets / block methods)
	if length(arg.subs) ~= 1, fail('"{?}" usage'), end

	out = fatrix2_subsref_brace(ob, arg.subs{1});

otherwise
	fail('subscript type %s unknown', arg.type)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function T = multiply_pots(T1, T2)
% MULTIPLY_POTS Multiply a pair of dpots together pointwise.
% T = multiply_pots(pots)

dom = myunion(T1.domain, T2.domain);
%ns = sparse(1, max(dom)); % causes problems in myreshape on NT
ns = zeros(1, max(dom));
ns(T1.domain) = T1.sizes;
ns(T2.domain) = T2.sizes;
T = dpot(dom, ns(dom));
T = multiply_by_pot(T, T1);
T = multiply_by_pot(T, T2);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function hhh=vline(x,in1,in2)
% function h=vline(x, linetype, label)
%
% Draws a vertical line on the current axes at the location specified by 'x'.  Optional arguments are
% 'linetype' (default is 'r:') and 'label', which applies a text label to the graph near the line.  The
% label appears in the same color as the line.
%
% The line is held on the current axes, and after plotting the line, the function returns the axes to
% its prior hold state.
%
% The HandleVisibility property of the line object is set to "off", so not only does it not appear on
% legends, but it is not findable by using findobj.  Specifying an output argument causes the function to
% return a handle to the line, so it can be manipulated or deleted.  Also, the HandleVisibility can be
% overridden by setting the root's ShowHiddenHandles property to on.
%
% h = vline(42,'g','The Answer')
%
% returns a handle to a green vertical line on the current axes at x=42, and creates a text object on
% the current axes, close to the line, which reads "The Answer".
%
% vline also supports vector inputs to draw multiple lines at once.  For example,
%
% vline([4 8 12],{'g','r','b'},{'l1','lab2','LABELC'})
%
% draws three lines with the appropriate labels and colors.
%
% By Brandon Kuczenski for Kensington Labs.
% brandon_kuczenski@kensingtonlabs.com
% 8 November 2001

if length(x)>1  % vector input
    for I=1:length(x)
        switch nargin
        case 1
            linetype='r:';
            label='';
        case 2
            if ~iscell(in1)
                in1={in1};
            end
            if I>length(in1)
                linetype=in1{end};
            else
                linetype=in1{I};
            end
            label='';
        case 3
            if ~iscell(in1)
                in1={in1};
            end
            if ~iscell(in2)
                in2={in2};
            end
            if I>length(in1)
                linetype=in1{end};
            else
                linetype=in1{I};
            end
            if I>length(in2)
                label=in2{end};
            else
                label=in2{I};
            end
        end
        h(I)=vline(x(I),linetype,label);
    end
else
    switch nargin
    case 1
        linetype='r:';
        label='';
    case 2
        linetype=in1;
        label='';
    case 3
        linetype=in1;
        label=in2;
    end

    g=ishold(gca);
    hold on

    y=get(gca,'ylim');
    h=plot([x x],y,linetype);
    if length(label)
        xx=get(gca,'xlim');
        xrange=xx(2)-xx(1);
        xunit=(x-xx(1))/xrange;
        if xunit<0.8
            text(x+0.01*xrange,y(1)+0.1*(y(2)-y(1)),label,'color',get(h,'color'))
        else
            text(x-.05*xrange,y(1)+0.1*(y(2)-y(1)),label,'color',get(h,'color'))
        end
    end

    if g==0
    hold off
    end
    set(h,'tag','vline','handlevisibility','off')
end % else

if nargout
    hhh=h;
end


text:
function b = isnull(x,eps)
% ckeck double == 0

if nargin == 1, eps = 1e-14; end

b = abs(x)<=eps;

author: mtex-toolbox

repo: mtex

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


text:
classdef ObjMeta
    %META Summary of this class goes here
    %   Detailed explanation goes here

    properties
        Name % file name
        Sensor % L_TM, L_ETM_PLUS, L_OLI, L_OLI_TIRS, S_MSI
        Resolution % Spatail Resolution.
        Dim
        UL
        Zen
        Azi
        ZC
        BBox % boundery.
        Output % output directory
    end

    methods
    end

end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
function [H_o_j, A_j, H_x_j] = calcHoj(p_f_G, msckfState, camStateIndices)
%CALCHOJ Calculates H_o_j according to Mourikis 2007
% Inputs: p_f_G: feature location in the Global frame
%         msckfState: the current window of states
%         camStateIndex: i, with camState being the ith camera pose in the window
% Outputs: H_o_j, A

N = length(msckfState.camStates);
M = length(camStateIndices);
H_f_j = zeros(2*M, 3);
H_x_j = zeros(2*M, 12 + 6*N);

c_i = 1;
for camStateIndex = camStateIndices
    camState = msckfState.camStates{camStateIndex};

    C_CG = quatToRotMat(camState.q_CG);
    %The feature position in the camera frame
    p_f_C = C_CG*(p_f_G - camState.p_C_G);

    X = p_f_C(1);
    Y = p_f_C(2);
    Z = p_f_C(3);

    J_i = (1/Z)*[1 0 -X/Z; 0 1 -Y/Z];

    H_f_j((2*c_i - 1):2*c_i, :) = J_i*C_CG;

    H_x_j((2*c_i - 1):2*c_i,12+6*(camStateIndex-1) + 1:12+6*(camStateIndex-1) + 3) = J_i*crossMat(p_f_C);
    H_x_j((2*c_i - 1):2*c_i,(12+6*(camStateIndex-1) + 4):(12+6*(camStateIndex-1) + 6)) = -J_i*C_CG;

    c_i = c_i + 1;
end

A_j = null(H_f_j');
H_o_j = A_j'*H_x_j;

end

author: utiasSTARS

repo: msckf-swf-comparison


text:
function [outStruct, replacedFields] = bst_jsondecode(inString, forceBstVersion)
% BST_JSONDECODE: Decodes a JSON string as a Matlab structure
%
% USAGE: outStruct = bst_jsondecode(inString)
%        outStruct = bst_jsondecode(filename)
%
% Limitations:
%   - requires the "root type" to be an object: a list of key:value pairs
%     enclosed in curly brackets.
%   - only supports arrays of number or string.
%   - escaped special characters \r\n\t are decoded as unicode arrow symbols as
%     MATLAB does. Other escaped characters (other than \\ and \") are kept as is
%     and therefore would not be correctly re-encoded.

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

if nargin < 2 || isempty(forceBstVersion)
    forceBstVersion = 0;
end

% Json allows any string (unicode characters) as a key, unlike Matlab var/field
% names. So some must be renamed.
iReplaced = 0;
replacedFields = cell(0, 2);

% If the input is an existing filename: read it
if exist(inString, 'file')
    jsonFile = inString;
    % Open file
    fid = fopen(jsonFile, 'r');
    if (fid < 0)
        error(['Cannot open JSON file: ' jsonFile]);
    end
    % Read file
    inString = fread(fid, [1, Inf], '*char');
    % Close file
    fclose(fid);
    % Check that something was read
    if isempty(inString)
        error(['File is empty: ' jsonFile]);
    end
else
    jsonFile = [];
end

% If possible, call built-in function
if ~forceBstVersion && exist('jsondecode', 'builtin') == 5
    outStruct = jsondecode(inString);
    return;
end

inString = strtrim(inString);
outStruct = struct();

% States:
STATE.START           = 0;
STATE.START_FIELD     = 1;
%STATE.READ_NEXT_FIELD = 2; % unused
STATE.READ_FIELD      = 3;
STATE.END_FIELD       = 4;
STATE.START_VALUE     = 5; % can also be the start of a string in an existing list
%STATE.READ_NEXT_VALUE = 6; % unused
STATE.READ_VALUE      = 7;
STATE.END_VALUE       = 8;
STATE.START_LIST      = 9;
STATE.END_LIST        = 10;
STATE.END_SUBLEVEL    = 11;

% Supported value types:
% Note: only lists of numbers or strings are supported
VAL.NONE = 0; % to simplify resetting without having to check if empty.
VAL.CHAR = 1;
VAL.NUM  = 2;
VAL.BOOL = 3;
VAL.LIST_NUM = 4;
VAL.LIST_STR = 5;

state = STATE.START;
field = ''; % json key / matlab struct field
value = [];
token = ''; % json string (key, value, or element of array value)
path  = {};
valType = VAL.NONE;
escape = false;
lineNum  = 1;
lineChar = 0;

for iChar = 1:length(inString)
    c = inString(iChar);
    lineChar = lineChar + 1;
    err = 0;

    % Process character depending on state

    % Reading a string
    if state == STATE.READ_FIELD || ...
            (state == STATE.READ_VALUE && (valType == VAL.CHAR || valType == VAL.LIST_STR))
        if escape
            if c == '"' || c == '\'
                token(end + 1) = c; %#ok<*AGROW>
                % Copy Matlab's hack to decode/encode escaped special characters.
            elseif c == 'n' % newline / linefeed
                token(end + 1) = char(8629); % ↵
            elseif c == 'r' % carriage return
                token(end + 1) = char(8592); % ←
            elseif c == 't' % tab
                token(end + 1) = char(8594); % →
            else
                % Possibly other special character. Keep as is, but won't be properly re-encoded.
                token = [token '\' c];
            end
            escape = false;
        elseif c == '\'
            % Escape character for itself and double quote
            escape = true;
        elseif c == '"'
            if state == STATE.READ_FIELD
                state = STATE.END_FIELD;
            else
                state = STATE.END_VALUE;
            end
        else
            % Any other unicode character is valid in a json string.
            token(end + 1) = c;
        end

        % Ignore all white space outside strings.
    elseif c == char(10)
        % Increment line number if we have a line break
        lineNum = lineNum + 1;
        lineChar = 0;
    elseif isspace(c)
        continue;

    elseif state == STATE.START
        if c == '{'
            state = STATE.START_FIELD;
        elseif c == '"'
            state = STATE.READ_FIELD;
        else
            err = 1;
        end
    elseif state == STATE.START_FIELD
        %field = '';
        token = '';
        if c == '"'
            state = STATE.READ_FIELD;
        elseif c == '}'
            % This was an empty "sublevel"/object, i.e. {}.
            state = STATE.END_SUBLEVEL;
        else
            err = 1;
        end
    elseif state == STATE.END_FIELD
        % Validate field name is appropriate for matlab.  Need to do it now in
        % case there are sub-fields.
        if ~isvarname(token)
            iReplaced = iReplaced + 1;
            %         field = ['REPLACED_FIELD_', str2double(iReplaced)];
            % Rename like Matlab, though 2 keys could have the same modified name.
            field = matlab.lang.makeValidName(token);
            % Append _ if the field already exists.
            tempPath = 'outStruct';
            for iP = 1:length(path)
                tempPath = [tempPath '.' path{iP}];
            end
            if isfield(eval(tempPath), field)
                field = [field '_'];
            end
            replacedFields(iReplaced, :) = {field, token};
        else
            field = token;
        end
        if c ~= ':'
            err = 1;
        end
        valType = VAL.NONE;
        value = [];
        state = STATE.START_VALUE;
    elseif state == STATE.START_VALUE
        % value was emptied in END_FIELD since we might continue a list of strings.
        if c == '{'
            % Save sub-field. Needed to later check existance of renamed fields.
            outStruct = saveField(outStruct, path, field, [], valType);
            path{end + 1} = field;
            state = STATE.START_FIELD;
        elseif c == '['
            state = STATE.START_LIST;
        elseif c == '"'
            % Check if continuing an existing list of strings
            if valType ~= VAL.LIST_STR
                valType = VAL.CHAR;
            end
            token = '';
            state = STATE.READ_VALUE;
        elseif c >= '0' && c <= '9'
            valType = VAL.NUM;
            token = c;
            state = STATE.READ_VALUE;
            % Detect special values for boolean and null
        elseif findAhead({'true', 'false', 'null'}, inString(iChar:end))
            valType = VAL.BOOL;
            token = c;
            state = STATE.READ_VALUE;
        else
            err = 1;
        end
    elseif state == STATE.READ_VALUE % NUM, BOOL or LIST_NUM
        if valType == VAL.LIST_NUM && (c == ',' || c == ']')
            n = str2double(token);
            token = '';
            if ~isempty(n)
                value(end + 1) = n;
            else
                err = 1;
            end
            if c == ']'
                state = STATE.END_LIST;
            end
        elseif c == ',' || c == '}'
            value = token;
            % Save
            outStruct = saveField(outStruct, path, field, value, valType);
            if c == ','
                state = STATE.START_FIELD;
            elseif c == '}'
                state = STATE.END_SUBLEVEL;
            end
        else % numerical, boolean or null values validated later.
            token(end + 1) = c;
        end
    elseif state == STATE.END_VALUE && valType == VAL.LIST_STR
        value{end + 1} = token;
        if c == ','
            state = STATE.START_VALUE;
        elseif c == ']'
            state = STATE.END_LIST;
        else
            err = 1;
        end
    elseif (state == STATE.END_VALUE && valType == VAL.CHAR) || ...
            state == STATE.END_LIST
        if valType == VAL.CHAR
            value = token;
        end
        % Save
        outStruct = saveField(outStruct, path, field, value, valType);
        if c == ','
            state = STATE.START_FIELD;
        elseif c == '}'
            state = STATE.END_SUBLEVEL;
        else
            err = 1;
        end
    elseif state == STATE.END_SUBLEVEL
        if ~isempty(path) % else this should be the end of the file
            path = path(1:end-1);
        end
        if c == ','
            state = STATE.START_FIELD;
        elseif c == '}'
            state = STATE.END_SUBLEVEL;
        else
            err = 1;
        end
    elseif state == STATE.START_LIST
        if c == '"'
            valType = VAL.LIST_STR;
            token = '';
        else
            valType = VAL.LIST_NUM;
            token = c;
        end
        state = STATE.READ_VALUE;
    else
        err = 1;
    end

    % Error management
    if err
        errorMsg = sprintf(...
            ['JSON syntax error. Unexpected character ''%c'' ' ...
            'at line %d, column %d (character %d)'], ...
            c, lineNum, lineChar, iChar);
        if ~isempty(jsonFile)
            errorMsg = [errorMsg ' of file ' jsonFile];
        end
        error(errorMsg);
    end

end % character loop

% Warn if there were invalid field names.
if iReplaced > 0
    if nargout > 1
        warning('One or more field names were not valid for MATLAB and were replaced.');
    else
        warning('One or more field names were not valid for MATLAB and were replaced. To see the original names, use bst_jsondecode with a second output argument.');
    end
end

    function outStruct = saveField(outStruct, path, field, value, valType)
        sPath = 'outStruct';
        for iNode = 1:length(path)
            sPath = [sPath '.' path{iNode}];
        end

        if isempty(value)
            value = 'struct';
        elseif valType == VAL.CHAR
            % Due to eval later, replace single quotes by 2 single quotes.
            value = ['''' strrep(value, '''', '''''') ''''];
        elseif valType == VAL.LIST_NUM
            elems = value;
            value = '[';
            for iElem = 1:length(elems)
                if iElem > 1
                    value = [value ','];
                end
                value = [value num2str(elems(iElem))];
            end
            value = [value ']'];
        elseif valType == VAL.LIST_STR
            elems = value;
            value = '{';
            for iElem = 1:length(elems)
                if iElem > 1
                    value = [value ','];
                end
                value = [value '''' strrep(elems{iElem}, '''', '''''') ''''];
            end
            value = [value '}'];
        elseif valType == VAL.BOOL
            value = lower(value);
            if strcmp(value, 'null')
                value = '[]';
            elseif ~strcmp(value, 'true') && ~strcmp(value, 'false')
                error(['JSON syntax error. Invalid boolean string: ' value]);
            end
        elseif valType == VAL.NUM && isempty(str2double(value))
            error(['JSON syntax error. Invalid numerical string: ' value]);
        end

        try
            eval([sPath '.' field ' = ' value ';']);
        catch ME
            sPath %#ok<*NOPRT>
            field
            value
            rethrow(ME);
        end

    end
end

function found = findAhead(needles, haystack)
found = 0;
for iNeedle = 1:length(needles)
    needle = needles{iNeedle};
    if strncmpi(needle, haystack, length(needle))
        found = 1;
        return;
    end
end
end

author: brainstorm-tools


text:
function [y, dzdx2] = vl_nncosineloss(x1, x2, c, varargin)
%VL_NNCOSINELOSS Compute cosine embedding loss
%   Y = VL_NNCOSINELOSS(X1, X2, C) computes the contrastive loss incurred by
%   the similar and dissimilar pairs in X1 and X2 labelled by C, making use
%   of the cosine similarity function to assess the distance between embeddings
%   (rather than euclidean metric used in the original contrastive loss
%   formulation of [1]).
%
%   The variables X1 nd X2 are of size H x W x D x N. The distance between
%   two pairs is computed between vectorised chunks of size HWD x N,
%   keeping the spatial and channel arrangement.
%
%   C has dimension 1 x 1 x 1 x N and specifies dissimilar pairs when
%   equal 0 and similar pairs otherwise.
%
%   The loss between two vectors X1 and X2 with cosine similarity
%   D = 1 - X1'X2 / (||X1||_2 *||X2||_2) and label C is computed as:
%   L(D, L) = sum(C * (1 - COS_SIM(X1,X2)) + (1-C) * max(COS_SIM(X1,X2) - M, 0)).
%
%   The term "cosine distance" is often used to denote the quantity
%   D = 1 - cos_sim(X1,X2), where cos_sim refers to the cosine similarity
%   between X1 and X2.  However, it is worth noting that this is not, strictly
%   speaking a distance function (it does not obey Cauchy-Schwarz). The
%   VL_NNCOSINELOSS seeks to minimise the cosine distance between similar pairs
%   and minimise cosine similarity between dissimilar pairs, up to a margin.
%
%   [DZDX1, DZDX2] = VL_NNCOSINELOSS(X1, X2, C, DZDY) computes the
%   derivative of the block projected onto the output derivative DZDY.
%   DZDX1, DZDX2 and DZDY have the same dimensions as X1, X2 and Y
%   respectively.
%
%   VL_NNCOSINELOSS(.., 'option', value, ...) accepts the following options:
%
%   `margin`:: 0.5
%    The maximum margin to be enforced between dissimilar pairs.
%
%    See also: VL_NNLOSS().
%
%    Based on the VL_NNCONTRLOSS() by Karel Lenc.
%
%    [1] Hadsell, Raia, Sumit Chopra, and Yann LeCun. "Dimensionality
%    reduction by learning an invariant mapping." CVPR 2006

% Copyright (C) 2018 Samuel Albanie
% Licensed under The MIT License [see LICENSE.md for details]

	opts.margin = 0.5 ;
  [opts, dzdy] = vl_argparsepos(opts, varargin) ;

	sx1 = size(x1) ; sx2 = size(x2) ; bsize = size(x1, 4) ;
	assert(numel(sx1) == numel(sx2), 'input sizes must match') ;
	assert(all(sx1 == sx2), 'Invalid input sizes.') ;
	assert(numel(c) == bsize, 'Invalid number of labels.') ;

  % allow element-specific margins
	if numel(opts.margin) > 1
		assert(numel(opts.margin) == bsize, 'Invalid margin.');
		opts.margin = reshape(opts.margin, [], bsize);
	end
	c = reshape(c, [], bsize);

  sims = vl_nncosinesim(x1, x2) ;
  diff = 1 - sims ;
	mdist = sims - opts.margin ;

	if isempty(dzdy)
		sims(c == 0) = max(mdist(c == 0), 0) ;
		y = sum(sims);
	else
    % outcome is invariant to order of projection
		[dcos1, dcos2] = vl_nncosinesim(x1, x2, dzdy{1}) ;
		keyboard
		one = ones(1, 'like', x1);
		mdist = squeeze(mdist);
		y1 = diff * (dzdy * 2);
		nf = mdist ./ (dist + 1e-4*one);
		neg_sel = mdist >  0 & c == 0;
		y1(:, neg_sel) = bsxfun(@times, -y1(:, neg_sel), nf(neg_sel));
		y1(:, mdist <= 0 & c == 0) = 0;
		y2 = -y1;
		y1 = reshape(y1, sx1);
		y2 = reshape(y2, sx2);
	end

repo: MFT


text:
function ecc = eccentricity(phi, params);
% From corAnal estimated phase (phi), and eccentricity params,
% compute the estimated eccentricity represented by a voxel or
% voxels, in visual degrees.
%
% ecc = eccentricity(phi, params);
%
% See retinoSetParams for info on the params struct.
%
% ASSUMPTIONS: This code makes the following assumptions:
% (1) Phi represents an accurate estimate of the zero-crossing point of
%     the best-fitting sinusoid for a given voxel's time series. This
%     is what computeCorAnal computes, but the co threshold should be
%     reasonable.
% (2) Phi corresponds to the time when the leading edge of a ring stimulus
%     -- if the ring is expanding, the outer edge; if contracting, the
%     inner edge -- just entered the receptive field of neurons within a
%     voxel (and so the neurons just started firing, and the hemodynamic
%     response started rising above the mean). This means the rise time
%     for the hemodynamic response should be about the same time as the
%     rise time of the sinusoid.
%
% ras, 01/06
if nargin < 2, error('Not enough input arguments.'); end

if isempty(phi), ecc = []; return; end

% check that all necessary fields are assigned
fields = fieldnames(params);
requiredFields = {'startAngle' 'endAngle' 'width' 'blankPeriod'};
if ~all(ismember(requiredFields, fields))
    error('Not all eccentricity params are set. Use retinoSetParams.')
end

% init output eccentricities -- will be same size as phi
ecc = zeros(size(phi));

% if specified, rotate phi mod 2*pi
if checkfields(params, 'startPhase') & params.startPhase~=0
	phi = mod(phi + 2*pi*(1-params.startPhase), 2*pi);
end

% figure out where the leading edge of the stimulus is w.r.t.
% the center angles specified by startAngle and endAngle. The
% idea here is that the rise time of the hemodynamic response
% -- which corresponds to the zero-crossing point of the fitted
% sinsuoid that the corAnal saves, and should be phi -- is
% determined by the leading edge of the stimulus, rather than
% the center.
if params.startAngle > params.endAngle
    dirFlag = -1;     % inward-moving stimulus
else
    dirFlag = +1;     % outward-moving stimulus
end

% the dynamic range of the cycle used for mapping eccentricity depends
% on the experimental design -- specifically, if there was a blank period
% the start or end of each cycle, we use the stimulus duty cycle to find
% this range, but if there was no blank period, or frequency tagging was
% used, we use the full input range of [0 2pi]. If the duty cycle is used,
% input phi values outside the range will be mapped to 0.
if isequal(lower(params.blankPeriod), 'start of cycle')
    rng = [2*pi*params.dutyCycle 2*pi];
elseif isequal(lower(params.blankPeriod), 'end of cycle')
    rng = [0 2*pi*(1-params.dutyCycle)];
else
    rng = [0 2*pi];
end
inRange = (phi>=rng(1) & phi<=rng(2));

% mark values not in range as -1
ecc(~inRange) = -1;

% the 'width' field of params can be a single, constant value in
% degrees, or a 2xN array specifying different widths at different
% time points within a cycle. In the latter case, the first row
% represents time within a cycle (0-2pi), and the second row represents
% width in degrees at that point.
% Compute the eccentricity separately for each instance:
if length(params.width)==1
    newRange = [params.startAngle params.endAngle]; % + dirFlag*params.width;
    ecc(inRange) = normalize(phi(inRange), newRange(1), newRange(2));

else
    xi = params.width(1,:);
    width = params.width(2,:);

    % interpolate to find the stimulus center angle at each point in xi
    delta = (params.startAngle - params.endAngle);
    xx = linspace(rng(1), rng(2), length(xi));
    yy = params.startAngle + xx .* delta ./ (2*pi);
    startAngle = interp1(xx, yy, xi);

    % subtract out the width at each point
    leadingEdge = startAngle + dirFlag*width;

    % now we need to interpolate back into the range
    % represented by phi:
    phiVals = unique(phi);
    eccVals = interp1(xi, leadingEdge, phiVals);
    for ii = 1:length(eccVals)
        I = find(phiVals==phiVals(i));
        ecc(I) = eccVals(ii);
    end
    ecc(~inRange) = 0;

end

return


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2010, John T. Ramshur, jramshur@gmail.com
%
% This file is part of HRVAS
%
% HRVAS 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.
%
% HRVAS 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 HRVAS.  If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function batchHRV(opt)
%batchHRV: function that calculates hrv on all ibi files in selected
%           directory
%
%Inputs:
%
%Outputs:
%

%check inputs
if (nargin < 1) || ~isstruct(opt)
    flagNoGui=true;
else
    flagNoGui=false;
end

if ~flagNoGui
    bcolor=[0.702 0.7216 0.8235]; %background color

    h.batchFigure=figure('Name','Batch Prcccess HRV', ...
        'Position',[20 50 500 500 ],...
        'Toolbar','none','Menubar','none');

    h.panel = uipanel('Parent',h.batchFigure,...
        'Units', 'normalized', 'Position',[0 0 1 1]);
    h.lbl=uicontrol(h.panel,'Style','text', ...
        'String','Choose Directory to perform batch HRV analysis:',...
        'HorizontalAlignment','left',...
        'FontWeight', 'bold', 'FontSize',12,...
        'Units', 'normalized', 'Position',[.02 .9 .8 .08]);
    h.lblDir=uicontrol(h.panel,'Style','text', 'String','',...
        'Units', 'normalized', 'Position',[.02 .85 .8 .05],...
        'BackgroundColor','white',...
        'HorizontalAlignment','left');
    h.btnChooseDir=uicontrol(h.panel,'Style','pushbutton', ...
        'String','Choose Dir',...
        'Units', 'normalized', 'Position',[.83 .85 .14 .05],...
        'Callback', @btnChooseDir_Callback);

    h.list = uicontrol(h.panel,'Style','list',...
        'Units', 'normalized','Position',[.02 .34 .8 .5],'max',3,'min',1,...
        'Callback',@list_Callback);
    h.lblProcessing=uicontrol(h.panel,'Style','edit', ...
        'String','!!! PROCESSING !!!',...
        'Units', 'normalized', 'Position',[.15 .55 .5 .1],...
        'BackgroundColor','white', 'FontSize',12, ...
        'ForegroundColor', 'red',...
        'HorizontalAlignment','center','visible','off');
    h.lblType = uicontrol(h.panel,'Style','text',...
                'String','File Type:','HorizontalAlignment','right',...
                'Units', 'normalized','Position',[.45 .25 .2 .05],...
                'Callback',@updateList);
    h.popupType = uicontrol(h.panel,'Style','popupmenu',...
                'String',{'Any (*.*)','ibi (*.ibi)','text (*.txt)'},...
                'Units', 'normalized','Position',[.67 .25 .15 .05],...
                'Value',2,'BackgroundColor','white','Callback',@updateList);

    h.lblHeader = uicontrol(h.panel,'Style','text',...
                'String','# of Header Rows:', ...
                'HorizontalAlignment','right',...
                'Units', 'normalized','Position',[.45 .21 .2 .04]);
    h.txtHeader = uicontrol(h.panel,'Style','edit',...
                'String','0',...
                'Units', 'normalized','Position',[.67 .21 .15 .04],...
                'BackgroundColor','white');

    h.btnProc = uicontrol(h.panel,'Style','pushbutton',...
                'String','Run Batch',...
                'Units', 'normalized','Position',[.02 .06 .8 .11],...
                'fontsize',10,...
                'Callback',@btnProc_Callback);
    h.chkAll = uicontrol(h.panel,'Style','checkbox',...
                'String','Include all files.',...
                'Units', 'normalized','Position',[.02 .02 .8 .03],...
                'fontsize',10,'value',1);

else %procede without GUI
    %choose file containing HRV analysis options
   reply = input('Enter full file path of HRV analysis option file: ', 's');
    if isempty(reply) || exist(reply,'file')~=2
        error('Please choose valid file path!')
        return;
    end
    opt=load(reply);
    opt=opt.settings;

    %choose dir containing input files
    reply = input('Enter dir path containing IBI files: ', 's');
    if isempty(reply) || exist(reply,'dir')~=7
        error('Please choose valid input dir path!')
        return;
    end

    %check input directory
    fileList = dir(fullfile(reply, '*.ibi')); %get list of files
    fileList(any([fileList.isdir],1))=[]; %remove folders/dir from the list
    fnames = {fileList.name}; %get file names only from fileList structure
    %build array of full file paths
    fpaths=cell(length(fnames),1);
    for ff=1:length(fnames)
        fpaths{ff}=fullfile(reply,fnames{ff});
    end

    %display results of input dir
    disp([reply ' contains ' num2str(length(fpaths)) ' .ibi files.'])

    %choose output file name to save hrv reslts
    outPath = input('Enter file name for exported HRV data: ', 's');
    if isempty(outPath)
        error('Please choose valid input dir path!')
        return;
    end

    %Batch process HRV
    hrv=batchGetHRV(fpaths,opt,flagNoGui);

    %Export/save HRV
    saveHRV(hrv,opt,fpaths,outPath)
    disp('batchHRV Done!')
end

    function btnChooseDir_Callback(hObject, eventdata)
        % Callback function executed when btnChooseDir is pressed

        %get directory path
      path = uigetdir('Select directory containg subject files to export:');
        if path~=0
            set(h.lblDir,'String',path)

            extVal=get(h.popupType,'value');
            switch extVal
                case 1
                    ext='*.*';
                case 2
                    ext='*.ibi';
                case 3
                    ext='*.txt';
            end
            fileList = dir(fullfile(path, ext)); %get list of files
            %remove folders/dir from the list
            fileList(any([fileList.isdir],1))=[];
            %get file names only from fileList structure
            fnames = {fileList.name};

            set(h.list,'String',fnames)
            set(h.chkAll,'value',1)
            str=['Run Batch (' num2str(size(fnames,2)) ' files)'];
            set(h.btnProc,'string',str)

        end
    end

    function updateList(hObject, eventdata)
        path = get(h.lblDir,'String');
        if ~isempty(path)
            extVal=get(h.popupType,'value');
            switch extVal
                case 1
                    ext='*.*';
                case 2
                    ext='*.ibi';
                case 3
                    ext='*.txt';
            end
            fileList = dir(fullfile(path, ext)); %get list of files
            %remove any folder/dir from the list
            fileList(any([fileList.isdir],1))=[];
            %get file names only from fileList structure
            fnames = {fileList.name};

            set(h.list,'String',fnames)
            set(h.chkAll,'value',1)
            str=['Run Batch (' num2str(size(fnames,2)) ' files)'];
            set(h.btnProc,'string',str)
        end
    end

    function btnProc_Callback(hObject, eventdata)
        [fname, pathname] = uiputfile(...
        {'*.csv','Comma Delimited (*.csv)';...
        '*.xlsx','Excel (*.xlsx)'},...
        'Export HRV as',...
        'hrv.csv');
        if isequal(fname,0) || isequal(pathname,0)
            error('Please choose a export file name.')
            return
        else
            outPath=fullfile(pathname,fname);
        end

        set(h.lblProcessing,'visible','on');
        drawnow;

        %get list of files
        fnames=get(h.list,'String');
        if get(h.chkAll,'value')==0 %
            iselected = get(h.list,'Value');
            fnames = fnames(iselected);
        end
        if isempty(fnames)
            error('No files to process.')
            return
        end
        %build array of full file paths
        path=get(h.lblDir,'String');
        fpaths=cell(length(fnames),1);
        for f=1:length(fnames)
            fpaths{f}=fullfile(path,fnames{f});
        end

        %get HRV
        hrv=batchGetHRV(fpaths,opt,flagNoGui);

        %Export HRV
        saveHRV(hrv,opt,fpaths,outPath)
        set(h.lblProcessing,'visible','off');
        drawnow;

        clear hrv
    end

    function list_Callback(hObject, eventdata)
        %get count of selected
        fnames=get(h.list,'String');
        iselected = get(h.list,'Value');
        fnames = fnames(iselected);

        if ~isempty(fnames)
            set(h.chkAll,'value',0)
        else
            set(h.chkAll,'value',1)
        end

        if get(h.chkAll,'value')
            str=['Run Batch (' num2str(sum(iselected)) ' files)'];
        else
            str=['Run Batch (' num2str(size(fnames,1)) ' files)'];
        end
        set(h.btnProc,'string',str)

    end
end

function saveHRV(hrv,opt,fList,outPath)
       exportHRV(outPath,fList,hrv,opt)
end

function hrv=batchGetHRV(fList,opt,flagNoGui)
    %Preallocate hrv array
    %hrv=repmat(struct('b',0,'a',0),1,length(fList));

    %create waitbar to display time remaining
    if ~flagNoGui
        abort=false;
        h=waitbar(0,'Processing: 0% complete.','CreateCancelBtn',@cancelcb);
        pause(0.05); drawnow;
    else
        disp('Processing: 0% complete.')
    end
    tic; avgElapsed=0;

    %get hrv for each file
    nFiles=length(fList);
    for f=1:nFiles
        hrv(f)=getHRV(fList{f},opt);

        %calculate remaining time for waitbar
        elapsedTime=toc; %total time elapsed since start
        %average time elapsed to process a single file
        avgElapsed=elapsedTime/f;
        %time remaining in processing (min)
        timeRem=ceil(avgElapsed*(nFiles-f)/60);
        %Update waitbar
         bar=round(f/nFiles*100);
        if ~flagNoGui
            waitbar(bar/100,h,['Processing: ' num2str(bar)  ...
                '% complete. (~' num2str(timeRem) ' min)'])
            drawnow;
            if abort
                warning('Batch canceled.');
                delete(h)
                return;
            end
        else
            disp(['Processing: ' num2str(bar) ...
                '% complete. (~' num2str(timeRem) ' min)'])
        end
    end
    if ~flagNoGui; delete(h); end;

    function cancelcb(a1, a2)
        abort = true;
    end
end

function output=getHRV(f,settings)

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % LOAD IBI
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    nIBI=[]; dIBI=[];
    IBI=loadIBI(f,settings);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Preprocess Data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %build cell array of locate artifacts methods
    methods={}; methInput=[];
    if settings.ArtLocatePer
        methods=[methods,'percent'];
        methInput=[methInput,settings.ArtLocatePerVal];
    end
    if settings.ArtLocateSD
        methods=[methods,'sd'];
        methInput=[methInput,settings.ArtLocateSDVal];
    end
    if settings.ArtLocateMed
        methods=[methods,'median'];
        methInput=[methInput,settings.ArtLocateMedVal];
    end
    %determine which window/span to use
    if strcmpi(settings.ArtReplace,'mean')
        replaceWin=settings.ArtReplaceMeanVal;
    elseif strcmpi(settings.ArtReplace,'median')
        replaceWin=settings.ArtReplaceMedVal;
    else
        replaceWin=0;
    end

    %Note: We don't need to use all the input arguments,but we will let the
    %function handle all inputs
    [dIBI,nIBI,trend,art] = preProcessIBI(IBI, ...
        'locateMethod', methods, 'locateInput', methInput, ...
        'replaceMethod', settings.ArtReplace, ...
        'replaceInput',replaceWin, ...
        'detrendMethod', settings.Detrend, ...
        'smoothMethod', settings.SmoothMethod, ...
        'smoothSpan', settings.SmoothSpan, ...
        'smoothDegree', settings.SmoothDegree, ...
        'polyOrder', settings.PolyOrder, ...
        'waveletType', ...
        [settings.WaveletType num2str(settings.WaveletType2)], ...
        'waveletLevels', settings.WaveletLevels, ...
        'lambda', settings.PriorsLambda,...
        'resampleRate',settings.Interp,...
        'meanCorrection',true);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Calculate HRV
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    output.ibiinfo.count=size(IBI,1); % total # of ibi
    output.ibiinfo.outliers=sum(art); % number of outliers

    %Time-Domain (using non-detrented ibi)
    output.time=timeDomainHRV(nIBI,settings.SDNNi*60,settings.pNNx);
    %output.time.mean=round(mean(nIBI(:,2).*1000)*10)/10;
    %output.time.meanHR=round(mean(60./nIBI(:,2))*10)/10;

    %Freq-Domain
    output.freq=freqDomainHRV(dIBI,settings.VLF,settings.LF, ...
        settings.HF,settings.AROrder,settings.WinWidth, ...
        settings.WinOverlap,settings.Points,settings.Interp);

    %Nonlinear (using non-detrented ibi)
    output.nl=nonlinearHRV(nIBI,settings.m,settings.r, ...
        settings.n1,settings.n2,settings.breakpoint);
    %Poincare
    output.poincare=poincareHRV(nIBI);

    %Time-Freq
    output.tf=timeFreqHRV(dIBI,nIBI,settings.VLF,settings.LF, ...
        settings.HF,settings.AROrder, settings.tfWinSize, ...
        settings.tfOverlap,settings.Points,settings.Interp, ...
        {'ar','lomb','wavelet'});
    %%%%%%%%%%%%%%%%%%%%%%%

    clear  nIBI dIBI IBI

end

function ibi=loadIBI(f,opt)
    if ~exist(f,'file')
        error(['Error opening file: ' f])
        return
    end

    ibi=[];
    DELIMITER = ',';
    HEADERLINES = opt.headerSize;

    %read ibi
    tmpData = importdata(f, DELIMITER, HEADERLINES);
    if HEADERLINES>0
        tmpData=tmpData.data;
    end

    %check ibi dimentions
    [rows cols] = size(tmpData);
    if rows==1
        tmpData=tmpData';
        ibi=zeros(cols,2);
        tmp=cumsum(tmpData);
        ibi(2:end,1)=tmp(1:end-1);
        ibi(:,2)=tmpData;
    elseif cols==1
        ibi=zeros(rows,2);
        tmp=cumsum(tmpData);
        ibi(2:end,1)=tmp(1:end-1);
        ibi(:,2)=tmpData;
    elseif rows<cols
        ibi=tmpData';
    else
        ibi=tmpData;
    end
    clear tmpData
end

author: jramshur

repo: HRVAS

save_path: github-repos/MATLAB/jramshur-HRVAS


text:
function vw = viewSet(vw,param,val,varargin)
%Organize methods for setting view parameters.
%
%   vw = viewSet(vw,param,val,val,varargin)
%
% Example:
%   vw = viewSet(vw, 'cothresh', 0.1);
%
% Author: Wandell
%
% See also viewSet, viewMapParameterField
%
% These are the cases that can be set (as of July 5, 2011):
%
% %%%%% Session-related properties; selected scan, slice, data type
%      'homedir'
%      'sessionname'
%      'subject'
%      'name'
%      'viewtype'
%      'subdir'
%      'curdt'
%      'curslice'
%      'curscan'
%      'datavalindex'
%
% %%%%% Traveling-Wave / Coherence Analysis properties
%      'coherence'
%      'scanco'
%      'phase'
%      'scanph'
%      'amplitude'
%      'scanamp'
%      'phwin'
%      'cothresh'
%      'refph'
%      'ampclip'
%      'framestouse'
%      'ncycles'
%
% %%%%% Map properties
%      'map'
%      'mapname'
%      'mapunits'
%      'mapclip'
%      'mapwin'
%      'scanmap'
%
% %%%%% Anatomy / Underlay-related properties
%      'anatomy'
%      'brightness'
%      'contrast'
%      'inplanepath'
%      'anatinitialize'
%      'anatomynifti'
%      'inplaneorientation'
%
% %%%%% ROI-related properties
%      'roi'
%      'rois'
%      'selectedroi'
%      'roioptions'
%      'filledperimeter'
%      'maskrois'
%      'roivertinds'
%      'showrois'
%      'hidevolumerois'
%      'roidrawmethod'
%      'roiname'
%      'roicoords'
%      'roimodified'
%      'roicomments'
%
% %%%%% Time-series related properties
%      'tseries'
%      'tseriesslice'
%      'tseriesscan'
%
% %%%%% Retinotopy/pRF Model related properties
%      'rmfile'
%      'rmmodel'
%      'rmparams'
%      'rmstimparams'
%      'rmmodelnum'
%
% %%%% Mesh-related properties
%      'leftclassfile'
%      'rightclassfile'
%      'leftgrayfile'
%      'rightgrayfile'
%
% %%%%% these params interface with the mrMesh functions
%      'mesh'
%      'currentmesh'
%      'allmeshes'
%      'addmesh'
%      'meshdata'
%      'meshn'
%      'deletemesh'
%      'recomputev2gmap'
%
% %%%%% Volume/Gray-related properties
%      'nodes'
%      'edges'
%      'allleftnodes'
%      'allleftedges'
%      'allrightnodes'
%      'allrightedges'
%      'coords'
%      'mmpervox'
%
% %%%%% UI properties
%      'initdisplaymodes'
%      'ui'
%      'anatomymode'
%      'coherencemode'
%      'correlationmode'
%      'phasemode'
%      'amplitudemode'
%      'projectedamplitudemode'
%      'mapmode'
%      'displaymode'
%      'phasecma'
%      'cmap'
%      'locs'
%      'flipud'

% ras, 06/07: eliminated subfunctions for each view (ipGet, volGet,
%			  flatGet); all one big SWITCH statement now.
%
% JW: 2/2011: viewSet and viewSet now take the input parameter field and
%               call viewMapParameterField before the long switch/case.
%               This function returns a standardized parameter field name.
%               It removes spaces and captials, and maps multiple names
%               (aliases) onto a single name. If you add a new parameter
%               (new case) to viewSet or viewSet, please use only a single
%               standardized parameter name in the viewSet and viewSet
%               functions. You can put as many aliases as you like in
%               viewMapParameterField.
%
% For example:
%    viewMapParameterField('curdt')
% and
%   viewMapParameterField('Current Data TYPE')
% both return 'curdt'. This means that
%    viewSet(vw, 'curdt')
% and
%    viewSet(vw, 'Current Data TYPE')
% are equivalent. Hence viewSet and viewSet should have the case 'curdt'.
% They do not need the case 'Current Data TYPE' or 'currentdatatype'.
%
%

if ~exist('vw', 'var'),  error('No view defined.'); end
if notDefined('val'), val = []; end

%%%%%%%%%%%%%%%%%%%%%%%%
% Big SWITCH Statement %
%%%%%%%%%%%%%%%%%%%%%%%%
% This statement will check for all possible properties,
% for all view types.
mrGlobals;

if ischar(vw)
    %This means that we are using the new functionality to list the
    %parameter set

	vw = mrvParamFormat(vw);

    %Check to see if we are asking for just one or all parameters:
    if ~exist('param','var'), viewMapParameterField(vw, 1);
    else viewMapParameterField(vw,1,mrvParamFormat(param));
    end
 %Using the new functionality
    return %early since we don't want to go through the rest
end %if

if notDefined('param'), error('No parameter defined'); end

%Format the parameter as lowercase and without spaces
param = mrvParamFormat(param);

% Standardize the name of the parameter field with name-mapping function
param = viewMapParameterField(param);
splitParam = viewMapParameterSplit(param);

switch splitParam
    case 'anatomy'
        vw = viewSetAnatomy(vw,param,val,varargin{:});
    case 'colorbar'
        vw = viewSetColorbar(vw,param,val,varargin{:});
    case 'em'
        vw = viewSetEm(vw,param,val,varargin{:});
    case 'flat'
        vw = viewSetFlat(vw,param,val,varargin{:});
    case 'map'
        vw = viewSetMap(vw,param,val,varargin{:});
    case 'mesh'
        vw = viewSetMesh(vw,param,val,varargin{:});
    case 'retinotopy'
        vw = viewSetRetinotopy(vw,param,val,varargin{:});
    case 'roi'
        vw = viewSetROI(vw,param,val,varargin{:});
    case 'session'
        vw = viewSetSession(vw,param,val,varargin{:});
    case 'timeseries'
        vw = viewSetTimeSeries(vw,param,val,varargin{:});
    case 'travelingwave'
        vw = viewSetTravelingWave(vw,param,val,varargin{:});
    case 'ui'
        vw = viewSetUI(vw,param,val,varargin{:});
    case 'volume'
        vw = viewSetVolume(vw,param,val,varargin{:});
    otherwise
        error('Unknown viewSet parameter');
end %switch

return;


text:
function [C,OutputFiles] = cifti(obj)
% Extract CIFTI-2 extension from a NIfTI-2 file and export data
%__________________________________________________________________________
% Copyright (C) 2019 Wellcome Trust Centre for Neuroimaging

C = [];
OutputFiles = {};

%-Check input is a CIFTI file
%--------------------------------------------------------------------------
if ~isfield(obj.hdr,'ext') || isempty(obj.hdr.ext) || obj.hdr.ext.ecode ~= 32
    return;
end

%-Parse CIFTI extension
%--------------------------------------------------------------------------
C.xml = char(obj.hdr.ext.edata(:)');
C.hdr = convert(xmltree(C.xml),[],'attributes',true);
s     = size(obj.dat); s = s(5:end);
C.dat = reshape(obj.dat,s);

%-Check input is a CIFTI-2 file
%--------------------------------------------------------------------------
switch C.hdr.attributes.Version
    case {'2'}
        % noop
    case {'1','1.0'}
        warning('CIFTI-1 file not supported.');
        return;
    otherwise
        warning('Unknown CIFTI version: %s.',C.hdr.attributes.Version);
        return;
end

%-Standard CIFTI Mapping Combinations
%==========================================================================

%-Dense Scalar (.dscalar.nii) or Dense Data Series (.dtseries.nii)
%--------------------------------------------------------------------------
if ismember(obj.hdr.intent_code,[3002, 3006]) % CONNECTIVITY_DENSE_SCALARS
                                              % CONNECTIVITY_DENSE_SERIES
    MapNames = {};
    Surfaces = {};
    Volumes  = {};
    vol      = struct;

    %-Read mapping information for each dimension
    %----------------------------------------------------------------------
    for i=1:numel(C.hdr.Matrix.MatrixIndicesMap)
        map = C.hdr.Matrix.MatrixIndicesMap{i};
        switch map.attributes.AppliesToMatrixDimension
            case '0'
                switch map.attributes.IndicesMapToDataType
                    case 'CIFTI_INDEX_TYPE_SCALARS'
                        MapNames = map.NamedMap;
                        if isstruct(MapNames), MapNames = {MapNames}; end
                        MapNames = cellfun(@(x) x.MapName,MapNames,'UniformOutput',false);
                    case 'CIFTI_INDEX_TYPE_SERIES'
                        MapNames = str2double(map.attributes.NumberOfSeriesPoints);
                        MapNames = arrayfun(@(x) sprintf('%04d',x),1:MapNames,'UniformOutput',false);
                    otherwise
                        error('Mismatch between NIfTI intent code and CIFTI index type.');
                end
            case '1'
                if ~strcmp(map.attributes.IndicesMapToDataType,'CIFTI_INDEX_TYPE_BRAIN_MODELS')
                    error('Mismatch between NIfTI intent code and CIFTI index type.');
                end
                if isstruct(map.BrainModel), map.BrainModel = {map.BrainModel}; end
                for j=1:numel(map.BrainModel)
                    switch map.BrainModel{j}.attributes.ModelType
                        case 'CIFTI_MODEL_TYPE_SURFACE'
                            nV  = str2double(map.BrainModel{j}.attributes.SurfaceNumberOfVertices);
                            off = str2double(map.BrainModel{j}.attributes.IndexOffset);
                            iV  = str2num(map.BrainModel{j}.VertexIndices) + 1;
                            if numel(iV) ~= str2double(map.BrainModel{j}.attributes.IndexCount)
                                error('Problem with number of vertices.');
                            end
                            BrainStructure = map.BrainModel{j}.attributes.BrainStructure(17:end);
                            Surfaces{end+1} = struct('nV',nV,'off',off,'iV',iV,'brain',BrainStructure);
                        case 'CIFTI_MODEL_TYPE_VOXELS'
                            off = str2double(map.BrainModel{j}.attributes.IndexOffset);
                            BrainStructure = map.BrainModel{j}.attributes.BrainStructure(17:end);
                            iV  = str2num(map.BrainModel{j}.VoxelIndicesIJK) + 1;
                            iV  = reshape(iV,3,[]);
                            if size(iV,2) ~= str2double(map.BrainModel{j}.attributes.IndexCount)
                                error('Problem with number of voxels.');
                            end
                            Volumes{end+1} = struct('off',off,'iV',iV,'brain',BrainStructure);
                    end
                end
                if ~isempty(Volumes)
                    vol.dim = str2num(map.Volume.attributes.VolumeDimensions);
                    vol.mat = str2num(map.Volume.TransformationMatrixVoxelIndicesIJKtoXYZ);
                    vol.mat = reshape(vol.mat,[4 4])' * [eye(4,3) [-1 -1 -1 1]'];
                    % assumes MeterExponent is -3
                end
            otherwise
                error('Data have to be dense scalars or series.');
        end
    end

    %-Export as GIfTI and NIfTI
    %----------------------------------------------------------------------
    for i=1:numel(MapNames)
        for j=1:numel(Surfaces)
            D = NaN(Surfaces{j}.nV,1);
            D(Surfaces{j}.iV) = C.dat(i,Surfaces{j}.off+(1:numel(Surfaces{j}.iV)));
            g = gifti(D);
            OutputFiles{end+1} = sprintf('%s_%s.gii',MapNames{i},Surfaces{j}.brain);
            save(g,OutputFiles{end},'ExternalFileBinary');
        end
        if ~isempty(Volumes)
            OutputFiles{end+1} = sprintf('%s.nii',MapNames{i});
            N = nifti; % using subsasgn as being used within the @nifti class
            N = subsasgn(N,substruct('.','dat'),file_array(OutputFiles{end},vol.dim,'float32-le',352));
            N = subsasgn(N,substruct('.','mat'),vol.mat);
            N = subsasgn(N,substruct('.','mat0'),vol.mat);
            N = subsasgn(N,substruct('.','mat_intent'),'Aligned');
            N = subsasgn(N,substruct('.','mat0_intent'),'Aligned');
            create(N);
            dat = NaN(vol.dim);
            for j=1:numel(Volumes)
                ind = sub2ind(vol.dim,Volumes{j}.iV(1,:),Volumes{j}.iV(2,:),Volumes{j}.iV(3,:));
                dat(ind) = C.dat(i,Volumes{j}.off+(1:numel(ind)));
            end
            N.dat = subsasgn(N.dat,substruct('()',repmat({':'},1,numel(vol.dim))),dat);
        end
    end
end

author: spm


text:
function OBJ = objFVN(FV, N)
%Convert Matlab faces/vertices/normals to a wavefront OBJ format
%
%  This appears to work for a cortical mesh.
%  We have also made it work for a fiber produced by AFQ.
%
% Example
%  OBJ = objFVN(FV,N);
%  objWrite(OBJ,fname);
%
% LMP/BW Vistasoft Team, 2015

%% Check FV and Normals

%% Figure out the general things to do here

% Need to understand more about the OBJ file format to make this work
% better.

% Make material structures.
% Ka is the ambient illuminant
% Kd is the diffuse reflectance
% Ks is the specular reflectance
% I don't know what the values in illum (presumably illuminant) and Ns (not
% sure what that is) refer to.
material(1).type='newmtl';
material(1).data='skin';
material(2).type='Ka';
material(2).data=[0.8 0.4 0.4];
material(3).type='Kd';
material(3).data=[0.8 0.4 0.4];
material(4).type='Ks';
material(4).data=[1 1 1];
material(5).type='illum';
material(5).data = 2;
material(6).type = 'Ns';
material(6).data = 27;

%% Make alias wavefront object structure

% This is the OBJ data structure that will be written out by the routine
% objWrite
clear OBJ
OBJ.material = material;

% Vertices and normals
OBJ.vertices = FV.vertices;
OBJ.vertices_normal = N;

% No idea what these control.  They don't seem aligned with the materials
% above.  Maybe they should be and then something good would happen.
OBJ.objects(1).type='g';
OBJ.objects(1).data='skin';

OBJ.objects(2).type='usemtl';
OBJ.objects(2).data='skin';

OBJ.objects(3).type='f';
OBJ.objects(3).data.vertices=  FV.faces;
OBJ.objects(3).data.normal  =  FV.faces;

end


text:
function IW = get_BMS_capacity(W)
%This function is employed here to verify the correctness of the program
%that is being writen
HYX = sum(W(1, :).*log2(W(1, :)));
PY = (W(1, :) + W(2, :))/2;
HY = sum(-PY.*log2(PY));
IW = HY + HYX;
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
function KG = kappa(I)
% get curvature information of input image
% input: 2D image I
% output: curvature matrix KG

% Copyright (c) 2009,
% Yue Wu @ ECE Department, Tufts University
% All Rights Reserved

I = double(I);
[m,n] = size(I);
P = padarray(I,[1,1],1,'pre');
P = padarray(P,[1,1],1,'post');

% central difference
fy = P(3:end,2:n+1)-P(1:m,2:n+1);
fx = P(2:m+1,3:end)-P(2:m+1,1:n);
fyy = P(3:end,2:n+1)+P(1:m,2:n+1)-2*I;
fxx = P(2:m+1,3:end)+P(2:m+1,1:n)-2*I;
fxy = 0.25.*(P(3:end,3:end)-P(1:m,3:end)+P(3:end,1:n)-P(1:m,1:n));
G = (fx.^2+fy.^2).^(0.5);
K = (fxx.*fy.^2-2*fxy.*fx.*fy+fyy.*fx.^2)./((fx.^2+fy.^2+eps).^(1.5));
KG = K.*G;
KG(1,:) = eps;
KG(end,:) = eps;
KG(:,1) = eps;
KG(:,end) = eps;
KG = KG./max(max(abs(KG)));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [p,c,v] = matrixpolynomial(x,n,dmax,dmin)
%MATRIXPOLYNOMIAL Creates parameterized polynomial
%
% [p,c,v] = matrixpolynomial(x,n,dmax,dmin)
%
% MATRIXPOLYNOMIAL is a quick way to define a parameterized polynomial
% p=sum C_i v_i(x) with all monomials of dmin <= degree(p,x) <= dmax. The
% coefficients in the polynomial are C_i while v is the monomial basis.
%
% Example:
%
% Paramterized quartic 2x2 matrix
%  x = sdpvar(2,1);
%  p = matrixpolynomial(x,2,4);
%
% See also MONOLIST, COEFFICIENTS

if (length(dmax) > 1) && (length(dmax) ~= length(x))
    error('Dimension mismatch: The third argument should be the max degree for each variable, or a sclar');
end

if nargin > 3
    if (length(dmin) > 1) && (length(dmin) ~= length(x))
        error('Dimension mismatch: The third argument should be the max degree for each variable, or a sclar');
    end
end

if any(dmax < 0)
    error('Only non-negative polynomial degrees possible')
end

if nargin<4
    dmin = 0;
end

if any(dmin > dmax)
    error('Fourth argument (dmin) should not be larger than second argument (dmax)');
end

if any(dmin < 0)
    error('Only non-negative polynomial degrees possible')
end

if length(n)==1
    n = [n n];
end

v = monolist(x,dmax);

if dmin <= dmax & dmin>0
    s = nchoosek(length(x) + dmin-1,dmin-1);
    v = extsubsref(v,s+1:length(v));
end

p = 0;
c = [];
for i = 1:length(v)
    Ci = sdpvar(n(1),n(2));
    c = [c reshape(Ci,[],1)];
end
p = reshape(c*v,n(1),n(2));

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% BV_CLOSE  Closes a TCP connection to BrainVision Recorder
%     BV_CLOSE(H)
%     In:
%     h : handle to an existing BrainVision connection

% Author: Hal Greenwald, The MITRE Corporation, 29-NOV-2011

% 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 bv_close(h)
if ~h.initialized
    return;
end
disp('Cleaning up connection to BrainVision Recorder');
pnet(h.handle, 'close');
if evalin('base', ['exist(' 'sprintf(''%s'', h.name)' ')'])
    evalin('base',['clear ' sprintf('%s',h.name) ';']);
end

author: PatternRecognition

repo: OpenBMI


text:
function theCombos=genAllRDCombinations(n,t,startVal)
%%GENALLRDCOMBINATIONS Generate all length t combinations of n elements
%                   (n>=t>1) in a resolving door gray-code order. This is
%                   such that the difference between each successiuve
%                   combination involves changing a single number.
%
%INPUTS: n The number of items from which t items are chosen; n>0.
%        t The number of items chosen, t<=n.
% startVal This is zero or 1, indicating which value the value at which the
%          elements in the combinations can start. The default if omitted
%          or an empty matrix is passed is 0;
%
%OUTPUTS: theCombos A tXnumCombos matrix containing all possible
%                   combinations of values in the "revolving door" order.
%
%This function implements algorithm R of Chapter 7.2.1.3 of [1].
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4A: Combinatorial
%    Algorithms, Part I, Boston: Addison-Wesley, 2011.
%
%October 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(startVal))
    startVal=0;
end

if(n==t&&t==1)
    theCombos=startVal;
    return
end

numCombos=binomial(n,t);
theCombos=zeros(t,numCombos);

%Step R1, Initialization
c=[(0:(t-1)).';n];

for curCombo=1:numCombos
    %Step R2, visit the combination.
    theCombos(:,curCombo)=c(1:t);

    %Step R3
    if(mod(t,2)==1)
        if(c(1)+1<c(2))
            c(1)=c(1)+1;
            continue;
        else
            j=2;
            skip4=false;
        end
    else
        if(c(1)>0)
            c(1)=c(1)-1;
            continue;
        else
           j=2;
           skip4=true;
        end
    end

    while(1)
        if(skip4==false)
            %Step R4
            if(c(j)>=j)
                c(j)=c(j-1);
                c(j-1)=j-2;
                break;
            else
                j=j+1;
            end
        end
        skip4=false;

        %Step R5.
        if(c(j)+1<c(j+1))
            c(j-1)=c(j);
            c(j)=c(j)+1;
            break;
        end

        j=j+1;

        if(j<=t)
            continue;
        else
            if(startVal~=0)
                theCombos=theCombos+1;
            end

            return;
        end
    end
end

if(startVal~=0)
    theCombos=theCombos+1;
end

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary