text:
function p = stat_fisher_pcomb(p)
% Fisher's (1925) method for combination of independent p-values [1]
% Code adapted from Bailey and Gribskov (1998) [2]
%
% [1] Fisher RA (1925). Statistical methods for research workers (13th edition). London: Oliver and Boyd.
% [2] Bailey TL, Gribskov M (1998). Combining evidence using p-values: application to sequence homology searches. Bioinformatics, 14 (1) 48-54.
%
% Author: Peter Watson, http://imaging.mrc-cbu.cam.ac.uk/statswiki/FAQ/CombiningPvalues

product=prod(p);
n=length(p);
if n<=0
    error('pfast was passed an empty array of p-values')
elseif n==1
    p = product;
    return
elseif product == 0
    p = 0;
    return
else
    x = -log(product);
    t=product;
    p=product;
    for i = 1:n-1
        t = t * x / i;
        p = p + t;
    end
end

repo: REST


text:
% DC to Dc Buck Converter Parameters
% Dick Benson
% Copyright 2013 The MathWorks, Inc.

FET_Ron=.004;    % FET
FET_Snub= 1000;

L1=250e-9;       % Inductor
L1_R = 0.01;

Cout=1000e-6;    % Filter Cap
Cout_R = 0.0012;

Rload=1.0;       % Load Resistance

Vinput=12;       % Vin to converter

Fswitch=500e3;   % PWM switching frequency
Comp_Tau=100e-9; % PWM Comparitor time constant

Fs_ADC=2.5e6;    % Sampling rate of discrete time controller

author: Sable

repo: mcbench-benchmarks


text:
function PRshowx(x, ProbInfo,h)
%PRshowx Show the solution x in IR Tools
%
% PRshowx(x,ProbInfo)
% PRshowx(x,ProbInfo,h)
%
% This function uses the information stored in ProbInfo to display the
% solution x in the correct way for that particular test problem.
%
% An optional third parameters specifies a figure handle.
%
% See also: PRshowb

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% This file is part of the IR Tools package and is distributed under the
% 3-Clause BSD License. A separate license file should be provided as part
% of the package.

if nargin < 3
  h = [];
end
if isempty(h)
  h = gcf;
end

switch ProbInfo.xType
    case {'image2D'}
        switch ProbInfo.problemType
            case {'deblurring','tomography'}
                figure(h)
                ImageDisplayGui2D(reshape(x, ProbInfo.xSize), h)
        end
    case {'surf2D'}
        switch ProbInfo.problemType
            case{'nmr'}
                [T2g, T1g] = meshgrid(ProbInfo.T2, ProbInfo.T1);
                figure(h)
                surf(log10(T2g), log10(T1g), reshape(x, ProbInfo.xSize))
                shading interp
                axis tight
            case {'invinterp2'}
                figure(h)
                n = ProbInfo.xSize(1);
                points = linspace(0, 1, n);
                [Xg,Yg] = meshgrid(points);
                plot3(Xg(:), Yg(:), x, '.', 'color', [0 0 0.7])
        end
    case {'fem'}
        trisurf(ProbInfo.elmtab,ProbInfo.elmX,ProbInfo.elmY,x)
        axis([0 1 0 1 -0.1 1])
        shading interp
    otherwise
        error('Can only do image2D, surf2D, and fem')
end

repo: IRtools


text:
function state = asr_calibrate(X,srate,cutoff,blocksize,B,A,window_len,window_overlap,max_dropout_fraction,min_clean_fraction)
% Calibration function for the Artifact Subspace Reconstruction (ASR) method.
% State = asr_calibrate(Data,SamplingRate,Cutoff,BlockSize,FilterB,FilterA,WindowLength,WindowOverlap,MaxDropoutFraction,MinCleanFraction)
%
% The input to this data is a multi-channel time series of calibration data. In typical uses the
% calibration data is clean resting EEG data of ca. 1 minute duration (can also be longer). One can
% also use on-task data if the fraction of artifact content is below the breakdown point of the
% robust statistics used for estimation (50% theoretical, ~30% practical). If the data has a
% proportion of more than 30-50% artifacts then bad time windows should be removed beforehand. This
% data is used to estimate the thresholds that are used by the ASR processing function to identify
% and remove artifact components.
%
% The calibration data must have been recorded for the same cap design from which data for cleanup
% will be recorded, and ideally should be from the same session and same subject, but it is possible
% to reuse the calibration data from a previous session and montage to the extent that the cap is
% placed in the same location (where loss in accuracy is more or less proportional to the mismatch
% in cap placement).
%
% The calibration data should have been high-pass filtered (for example at 0.5Hz or 1Hz using a
% Butterworth IIR filter).
%
% In:
%   Data : Calibration data [#channels x #samples]; *zero-mean* (e.g., high-pass filtered) and
%          reasonably clean EEG of not much less than 30 seconds length (this method is typically
%          used with 1 minute or more).
%
%   SamplingRate : Sampling rate of the data, in Hz.
%
%
%   The following are optional parameters (the key parameter of the method is the RejectionCutoff):
%
%   RejectionCutoff: Standard deviation cutoff for rejection. Data portions whose variance is larger
%                    than this threshold relative to the calibration data are considered missing
%                    data and will be removed. The most aggressive value that can be used without
%                    losing too much EEG is 2.5. A quite conservative value would be 5. Default: 5.
%
%   Blocksize : Block size for calculating the robust data covariance and thresholds, in samples;
%               allows to reduce the memory and time requirements of the robust estimators by this
%               factor (down to Channels x Channels x Samples x 16 / Blocksize bytes). Default: 10
%
%   FilterB, FilterA : Coefficients of an IIR filter that is used to shape the spectrum of the signal
%                      when calculating artifact statistics. The output signal does not go through
%                      this filter. This is an optional way to tune the sensitivity of the algorithm
%                      to each frequency component of the signal. The default filter is less
%                      sensitive at alpha and beta frequencies and more sensitive at delta (blinks)
%                      and gamma (muscle) frequencies. Default:
%                      [b,a] = yulewalk(8,[[0 2 3 13 16 40 min(80,srate/2-1)]*2/srate 1],[3 0.75 0.33 0.33 1 1 3 3]);
%
%   WindowLength : Window length that is used to check the data for artifact content. This is
%                  ideally as long as the expected time scale of the artifacts but short enough to
%				   allow for several 1000 windows to compute statistics over. Default: 0.5.
%
%   WindowOverlap : Window overlap fraction. The fraction of two successive windows that overlaps.
%                   Higher overlap ensures that fewer artifact portions are going to be missed (but
%                   is slower). Default: 0.66
%
%   MaxDropoutFraction : Maximum fraction of windows that can be subject to signal dropouts
%                        (e.g., sensor unplugged), used for threshold estimation. Default: 0.1
%
%   MinCleanFraction : Minimum fraction of windows that need to be clean, used for threshold
%                      estimation. Default: 0.25
%
%
% Out:
%   State : initial state struct for asr_process
%
% Notes:
%   This can run on a GPU with large memory and good double-precision performance for faster processing
%   (e.g., on an NVIDIA GTX Titan or K20), but requires that the Parallel Computing toolbox is
%   installed.
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2012-08-31

% asr_calibrate_version<1.05> -- for the cache

% UC Copyright Notice
% This software is Copyright (C) 2013 The Regents of the University of California. All Rights Reserved.
%
% Permission to copy, modify, and distribute this software and its documentation for educational,
% research and non-profit purposes, without fee, and without a written agreement is hereby granted,
% provided that the above copyright notice, this paragraph and the following three paragraphs appear
% in all copies.
%
% Permission to make commercial use of this software may be obtained by contacting:
% Technology Transfer Office
% 9500 Gilman Drive, Mail Code 0910
% University of California
% La Jolla, CA 92093-0910
% (858) 534-5815
% invent@ucsd.edu
%
% This software program and documentation are copyrighted by The Regents of the University of
% California. The software program and documentation are supplied "as is", without any accompanying
% services from The Regents. The Regents does not warrant that the operation of the program will be
% uninterrupted or error-free. The end-user understands that the program was developed for research
% purposes and is advised not to rely exclusively on the program for any reason.
%
% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT,
% SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF
% THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES,
% INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
% PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF
% CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
% MODIFICATIONS.

[C,S] = size(X);

if nargin < 3 || isempty(cutoff)
    cutoff = 5; end
if nargin < 4 || isempty(blocksize)
    blocksize = 10; end
if nargin < 6 || isempty(A) || isempty(B)
    try
        % try to use yulewalk to design the filter (Signal Processing toolbox required)
        [B,A] = yulewalk(8,[[0 2 3 13 16 40 min(80,srate/2-1)]*2/srate 1],[3 0.75 0.33 0.33 1 1 3 3]);
    catch e %#ok<NASGU>
        % yulewalk not available (maybe no toolbox installed) -- use precomputed filter
        % coefficients depending on sampling rate
        switch srate
            case 100
                [B,A] = deal([0.9314233528641650 -1.0023683814963549 -0.4125359862018213  0.7631567476327510  0.4160430392910331 -0.6549131038692215 -0.0372583518046807  0.1916268458752655  0.0462411971592346],[1.0000000000000000 -0.4544220180303844 -1.0007038682936749  0.5374925521337940  0.4905013360991340 -0.4861062879351137 -0.1995986490699414  0.1830048420730026  0.0457678549234644]);
            case 128
                [B,A] = deal([1.1027301639165037 -2.0025621813611867  0.8942119516481342  0.1549979524226999  0.0192366904488084  0.1782897770278735 -0.5280306696498717  0.2913540603407520 -0.0262209802526358],[1.0000000000000000 -1.1042042046423233 -0.3319558528606542  0.5802946221107337 -0.0010360013915635  0.0382167091925086 -0.2609928034425362  0.0298719057761086  0.0935044692959187]);
            case 200
                [B,A] = deal([1.4489483325802353 -2.6692514764802775  2.0813970620731115 -0.9736678877049534  0.1054605060352928 -0.1889101692314626  0.6111331636592364 -0.3616483013075088  0.1834313060776763],[1.0000000000000000 -0.9913236099393967  0.3159563145469344 -0.0708347481677557 -0.0558793822071149 -0.2539619026478943  0.2473056615251193 -0.0420478437473110  0.0077455718334464]);
            case 256
                [B,A] = deal([1.7587013141770287 -4.3267624394458641  5.7999880031015953 -6.2396625463547508  5.3768079046882207 -3.7938218893374835  2.1649108095226470 -0.8591392569863763  0.2569361125627988],[1.0000000000000000 -1.7008039639301735  1.9232830391058724 -2.0826929726929797  1.5982638742557307 -1.0735854183930011  0.5679719225652651 -0.1886181499768189  0.0572954115997261]);
            case 300
                [B,A] = deal([1.9153920676433143  -5.7748421104926795   9.1864764859103936 -10.7350356619363630   9.6423672437729007  -6.6181939699544277   3.4219421494177711  -1.2622976569994351   0.2968423019363821],[1.0000000000000000 -2.3143703322055491  3.2222567327379434 -3.6030527704320621  2.9645154844073698 -1.8842615840684735  0.9222455868758080 -0.3103251703648485  0.0634586449896364]);
            case 500
                [B,A] = deal([2.3133520086975823 -11.9471223009159130  29.1067166493384340 -43.7550171007238190  44.3385767452216370 -30.9965523846388000  14.6209883020737190  -4.2743412400311449   0.5982553583777899],[1.0000000000000000  -4.6893329084452580  10.5989986701080210 -14.9691518101365230  14.3320358399731820  -9.4924317069169977   4.2425899618982656  -1.1715600975178280   0.1538048427717476]);
            case 512
                [B,A] = deal([2.3275475636130865 -12.2166478485960430  30.1632789058248850 -45.8009842020820410  46.7261263011068880 -32.7796858196767220  15.4623349612560630  -4.5019779685307473   0.6242733481676324],[1.0000000000000000  -4.7827378944258703  10.9780696236622980 -15.6795187888195360  15.1281978667576310 -10.0632079834518220   4.5014690636505614  -1.2394100873286753   0.1614727510688058]);
            otherwise
                error('asr_calibrate:NoYulewalk','The yulewalk() function was not found and there is no pre-computed spectral filter for your sampling rate. If you would like to use the default spectral filter please try to resample to one of the supported rates (100,128,200,256,300,500,512) or get the appropriate toobox license (you can also disable the spectral weighting feature or supply your own precalculated IIR filter coefficients).');
        end
    end
end
if nargin < 8 || isempty(window_len)
    window_len = 0.5; end
if nargin < 9 || isempty(window_overlap)
    window_overlap = 0.66; end
if nargin < 10 || isempty(max_dropout_fraction)
    max_dropout_fraction = 0.1; end
if nargin < 11 || isempty(min_clean_fraction)
    min_clean_fraction = 0.25; end

X(~isfinite(X(:))) = 0;

% apply the signal shaping filter and initialize the IIR filter state
[X,iirstate] = filter(B,A,double(X),[],2); X = X';
if any(~isfinite(X(:)))
    error('The IIR filter diverged on your data. Please try using either a more conservative filter or removing some bad sections/channels from the calibration data.'); end

% calaulate
% calculate the sample covariance matrices U (averaged in blocks of blocksize successive samples)
U = zeros(length(1:blocksize:S),C*C);
for k=1:blocksize
    range = min(S,k:blocksize:(S+k-1));
    U = U + reshape(bsxfun(@times,reshape(X(range,:),[],1,C),reshape(X(range,:),[],C,1)),size(U));
end

% get the mixing matrix M
M = sqrtm(real(reshape(geometric_median(U/blocksize),C,C)));

% window
% window length for calculating thresholds
N = round(window_len*srate);

% get the threshold matrix T
fprintf('Determining per-component thresholds...');
[V,D] = eig(M);
X = abs(X*V);
for c = C:-1:1
    % compute RMS amplitude for each window...
    rms = X(:,c).^2;
    rms = sqrt(sum(rms(bsxfun(@plus,round(1:N*(1-window_overlap):S-N),(0:N-1)')))/N);
    % fit a distribution to the clean part
    [mu(c),sig(c)] = fit_eeg_distribution(rms,min_clean_fraction,max_dropout_fraction);
end
T = diag(mu + cutoff*sig)*V';
disp('done.');

% initialize the remaining filter state
if length(B) > 9
    % for more than 8'th order we try to initialize a second-order-sections filter implementation
    try
        [SOS,G] = tf2sos(B,A);
    catch e
        fprintf('Cannot use second-order section filter due to error: %s\n',e.message);
        [SOS,G] = deal([]);
    end
else
    [SOS,G] = deal([]);
end
state = struct('M',M,'T',T,'B',B,'A',A,'SOS',SOS,'G',G,'Zi',{repmat({[]},1,size(SOS,1))},'cov',[],'carry',[],'iir',iirstate,'last_R',[],'last_trivial',true);

function y = block_geometric_median(X,blocksize,varargin)
% Calculate a blockwise geometric median (faster and less memory-intensive
% than the regular geom_median function).
%
% This statistic is not robust to artifacts that persist over a duration that
% is significantly shorter than the blocksize.
%
% In:
%   X : the data (#observations x #variables)
%   blocksize : the number of successive samples over which a regular mean
%               should be taken
%   tol : tolerance (default: 1.e-5)
%   y : initial value (default: median(X))
%   max_iter : max number of iterations (default: 500)
%
% Out:
%   g : geometric median over X
%
% Notes:
%   This function is noticably faster if the length of the data is divisible by the block size.
%   Uses the GPU if available.
%

if nargin < 2 || isempty(blocksize)
    blocksize = 1; end

if blocksize > 1
    [o,v] = size(X);       % #observations & #variables
    r = mod(o,blocksize);  % #rest in last block
    b = (o-r)/blocksize;   % #blocks
    if r > 0
        X = [reshape(sum(reshape(X(1:(o-r),:),blocksize,b*v)),b,v); sum(X((o-r+1):end,:))*(blocksize/r)];
    else
        X = reshape(sum(reshape(X,blocksize,b*v)),b,v);
    end
end

% try
try
    y = gather(geometric_median(gpuArray(X),varargin{:}))/blocksize;
catch
    y = geometric_median(X,varargin{:})/blocksize;
end

% geometric_median
function y = geometric_median(X,tol,y,max_iter)
% Calculate the geometric median for a set of observations (mean under a Laplacian noise distribution)
% This is using Weiszfeld's algorithm.
%
% In:
%   X : the data, as in mean
%   tol : tolerance (default: 1.e-5)
%   y : initial value (default: median(X))
%   max_iter : max number of iterations (default: 500)
%
% Out:
%   g : geometric median over X

if ~exist('tol','var') || isempty(tol)
    tol = 1.e-5; end
if ~exist('y','var') || isempty(y)
    y = median(X); end
if ~exist('max_iter','var') || isempty(max_iter)
    max_iter = 500; end

for i=1:max_iter
    invnorms = 1./sqrt(sum(bsxfun(@minus,X,y).^2,2));
    [y,oldy] = deal(sum(bsxfun(@times,X,invnorms)) / sum(invnorms),y);
    if norm(y-oldy)/norm(y) < tol
        break; end
end

author: PatternRecognition

repo: OpenBMI


text:
classdef CEC2017_F13 < PROBLEM
% <single> <real> <constrained>
% CEC'2017 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% G. Wu, R. Mallipeddi, and P. N. Suganthan, Problem definitions and
% evaluation criteria for the CEC 2017 competition on constrained real-
% parameter optimization, National University of Defense Technology, China,
% 2016.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;  % Optimal decision vector
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2017.mat'),'Data');
            obj.O = Data{12}.o;
            obj.M = 1;
            if isempty(obj.D); obj.D = 10; end
            obj.D = min(obj.D,length(obj.O));
            obj.lower    = zeros(1,obj.D) - 100;
            obj.upper    = zeros(1,obj.D) + 100;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Y = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopObj = sum(100*(Y(:,1:end-1).^2-Y(:,2:end)).^2+(Y(:,1:end-1)-1).^2,2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Y = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopCon(:,1) = sum(Y.^2-10*cos(2*pi*Y)+10,2) - 100;
            PopCon(:,2) = sum(Y,2) - 2*size(Y,2);
            PopCon(:,3) = 5 - sum(Y,2);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef test1DLHS < handle

    properties (Access = public)
        tol = 1e-12;
    end

    properties (Access = private)
        mesh
        K, M
        LHS, RHS
        funL2, funH1
    end

    methods (Access = protected, Static)

        function nu = createPoissonValue()
            nu = 0;
        end
    end

    methods (Access = public)

        function obj = test1DLHS()
            obj.init()
            obj.createMesh()
            obj.createLHS();
            obj.createRHS();
            obj.solveSystem();
        end

        function p = hasPassed(obj)
            x = load('test_1DLHS', 'x');
            p = isequal(x.x, obj.funH1.fValues);
        end

    end

    methods (Access = private)

        function init(obj)
        end

        function createMesh(obj)
            n = 50;
            xCoord =  0:1/(n-1):1;
            s.coord(:,1) = xCoord;
            s.coord(:,2) = zeros(n,1);
            s.connec(:,1) = 1:n-1;
            s.connec(:,2) = 2:n;
            s.kFace = -1;
            m = Mesh(s);
%             m.plot();
            obj.mesh = m;
        end

        function createLHS(obj)
            obj.createStifnessMatrix();
            obj.createMassMatrix();
            e = 0.1*obj.mesh.computeMeanCellSize();
            obj.LHS = e*obj.K + obj.M;
        end

        function createStifnessMatrix(obj)
            l.type = 'StiffnessMatrix';
            l.fun  = P1Function.create(obj.mesh,1);
            l.mesh = obj.mesh;
            lhs = LHSintegrator.create(l);
            obj.K = lhs.compute();
        end

        function createMassMatrix(obj)
            s.type = 'MassMatrix';
            s.fun  = P1Function.create(obj.mesh,1);
            s.mesh = obj.mesh;
            lhs    = LHSintegrator.create(s);
            obj.M  = lhs.compute();
        end

        function createRHS(obj)
            xCoord = obj.mesh.coord(:,1)';
            xmean = (max(xCoord)-min(xCoord))/2;
            f = heaviside(xCoord-xmean);

            s.mesh = obj.mesh;
            s.fValues = f;
            fL2 = P1Function(s);
%             fL2.plot;
            obj.RHS = obj.M*(fL2.fValues)';
            obj.funL2 = fL2;
        end

        function solveSystem(obj)
            s.type    =  'DIRECT';
            solver    = Solver.create(s);
            fH1Values = solver.solve(obj.LHS,obj.RHS);

            s.mesh = obj.mesh;
            s.fValues = fH1Values;
            fH1 = P1Function(s);
%             fH1.plot()
            obj.funH1 = fH1;
%             obj.mesh.plot();
%             hold on
%             xCoord = obj.mesh.coord(:,1)';
%             plot(xCoord,obj.funL2.fValues,'-+');
%             plot(xCoord,obj.funH1.fValues,'-+');
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function rayPlot(ray,ord)
%+========================================================================+
%|                                                                        |
%|           OPENRAY - LIBRARY FOR TRI-DIMENSIONAL RAY TRACING            |
%|           openRay is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : rayPlot.m                                     |
%|    #    |   VERSION    : 0.41                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 01.04.2018                                    |
%| ( === ) |   SYNOPSIS   : Plot trajectory and last direction            |
%|  `---'  |                                                              |
%+========================================================================+

% Graphical configuration
figure(gcf)
off = ~ishold;
if off
    hold on
end

% Plot trajectory
if (ord > -1)
    for i = 2:ord+1
        Pim1 = ray.pos{i-1};
        Pi   = ray.pos{i};
        plot3([Pim1(:,1) Pi(:,1)]',[Pim1(:,2) Pi(:,2)]',[Pim1(:,3) Pi(:,3)]','o-b')
    end
end

% Plot direction
if (ord == -1) || (length(ray.pos) == 1)
    P = ray.pos{end};
    U = ray.dir;
    quiver3(P(:,1),P(:,2),P(:,3),U(:,1),U(:,2),U(:,3),'r')
end

% Hold
if off
    hold off
end
end


text:
function [V] = radius2volume(r)
% radius2volume returns the volume of a sphere of radius r.
% Chad Greene 2012
V = (4*pi/3).*r.^3;


text:
function [qua, DCMbn, euler] = att_update(wb, DCMbn, qua, omega_ie_n, omega_en_n, dt, att_mode)
% att_update: updates attitude using quaternion or DCM.
%
% INPUT
%   wb,         3x1 incremental turn-rates in body-frame (rad/s).
%   DCMbn,      3x3 body-to-nav DCM.
%   qua,        4x1 quaternion.
%   omega_ie_n, 3x3 skew-symmetric Earth rate matrix (rad/s).
%   omega_en_n, 3x3 skew-symmetric transport rate (rad/s).
%   dt,         1x1 IMU sampling interval (s).
%	att_mode,   attitude mode (string).
%      'quaternion': attitude updated as quaternion. Default value.
%             'dcm': attitude updated as Direct Cosine Matrix.
%
% OUTPUT
%   qua,      4x1 updated quaternion.
%   DCMbn,    3x3 updated body-to-nav DCM.
%   euler,    3x1 updated Euler angles (rad).
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% Reference:
%
%	Titterton, D.H. and Weston, J.L. (2004). Strapdown
% Inertial Navigation Technology (2nd Ed.). Institution
% of Engineering and Technology, USA.
%
%	Crassidis, J.L. and Junkins, J.L. (2011). Optimal Esti-
% mation of Dynamic Systems, 2nd Ed. Chapman and Hall/CRC, USA.
% Eq. 7.39, p. 458.
%
% Version: 004
% Date:    2022/03/06
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

if nargin < 7, att_mode  = 'quaternion'; end

%% Gyros output correction for Earth and transport rates

om_ie_n = skewm_inv(omega_ie_n);
om_en_n = skewm_inv(omega_en_n);
wb = (wb - DCMbn' * (om_ie_n + om_en_n));  % Titterton, Eq. 3.29, p. 32

if strcmp(att_mode, 'quaternion')
%% Quaternion update

    qua   = qua_update(qua, wb, dt);    % Quaternion update
    qua   = qua / norm(qua);            % Brute-force normalization
    DCMbn = qua2dcm(qua);               % DCM update
    euler = qua2euler(qua);             % Euler angles update

elseif strcmp(att_mode, 'dcm')
%% DCM update

    delta_theta = wb * dt;                  % Incremental Euler angles
    DCMbn = dcm_update(DCMbn, delta_theta); % DCM update
    euler = dcm2euler(DCMbn);               % Euler angles update
    qua   = euler2qua(euler);               % Quaternion update
    qua   = qua / norm(qua);                % Brute-force normalization

else
    error('att_update: no attitude update mode defined. Check the attitude update mode selected.')
end

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
%SurfPointFeature.support Support region of feature
%
% OUT = F.support(IM, W) is an image of the support region of the
% feature F, extracted from the image IM in which the feature appears.
% The support region is scaled to WxW and rotated so that the feature's
% orientation axis is upward.
%
% OUT = F.support(IMAGES, W) as above but if the features were extracted
% from an image sequence IMAGES then the feature is extracted from the
% appropriate image in the same sequence.
%
% [OUT,T] = F.support(IMAGES, W) as above but returns the pose of the feature
% as a 3x3 homogeneous transform in SE(2) that comprises the feature position
% and orientation.
%
% F.support(IM, W) as above but the support region is displayed.
%
% See also SurfPointFeature.

function [out,TT] = support(sf, images, N)

    if nargin < 3
        N = 50;
    end

    if isempty(sf.image_id_)
        im = images(:,:);
    else
        im = images(:,:,sf.image_id_);
    end

    d = 20*sf.scale_;

    [Uo,Vo] = imeshgrid(N, N);

    T = SE2(sf.u_, sf.v_, sf.theta_) * SE2(diag([d/N,d/N,1])) * SE2(-N/2, -N/2);

    UV = T * [Uo(:) Vo(:)]';
    U = reshape(UV(1,:), size(Uo));
    V = reshape(UV(2,:), size(Vo));

    [Ui,Vi] = imeshgrid(im);

    im2 = interp2(Ui, Vi, idouble(im), U, V);

    if nargout == 0
        idisp(im2)
    elseif nargout == 1
        out = im2;
    elseif nargout == 2
        out = im2;
        TT = T;
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
 function printf(str, varargin)
%function printf(str, varargin)
% shorthand for disp(sprintf())
% note: fprintf() is overloaded to work without a file id, so i could have
% used it instead.  oh well.
% Jeff Fessler
if nargin < 1, ir_usage, end
disp(sprintf(str, varargin{:}))

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [orientim, G_xx, G_yy, G_xy, cos_2_theta, sin_2_theta, denom] = orient(or_im, g_sigma, b_sigma, o_smooth)

    % Calculate image gradients.
    hsize = fix(6*g_sigma);
    if ~mod(hsize,2); hsize = hsize+1; end
    f = fspecial('gaussian', hsize, g_sigma);
    [f_x,f_y] = gradient(f);

    G_x = filter2(f_x, or_im);
    G_y = filter2(f_y, or_im);

    G_xx = G_x.^2;
    G_xy = G_x.*G_y;
    G_yy = G_y.^2;

    % to smooth the covariance data

    hsize = fix(6*b_sigma);
    if ~mod(hsize,2); hsize = hsize+1; end
    f = fspecial('gaussian', hsize, b_sigma);

    G_xx = filter2(f, G_xx);
    G_xy = 2*filter2(f, G_xy);
    G_yy = filter2(f, G_yy);

    denom = sqrt(G_xy.^2 + (G_xx - G_yy).^2) + eps;
    sin_2_theta = G_xy./denom;            % Sin of angle
    cos_2_theta = (G_xx-G_yy)./denom;     % Cos of angle

    %Smooth
        hsize = fix(6*o_smooth);
        if ~mod(hsize,2); hsize = hsize+1; end
        f = fspecial('gaussian', hsize, o_smooth);
        cos_2_theta = filter2(f, cos_2_theta);
        sin_2_theta = filter2(f, sin_2_theta);

    a = pi/2;
    orientim = a + atan2(sin_2_theta,cos_2_theta)/2;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [x] = spm_gen_phase(P,M,U)
% Generate state activities for trial-specific phase-coupled activity
% FORMAT [x] = spm_gen_phase(P,M,U)
%
% P - parameters
% M - model structure
% U - trial-specific effects
%
% x - states
%
%__________________________________________________________________________
% Copyright (C) 2009 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_gen_phase.m 2908 2009-03-20 14:54:03Z will $

% between-trial inputs
%==========================================================================
try, X = U.X; catch, X = sparse(1,0); end

if isfield(M,'trial')
    Mx=M.x;
end

% cycle over trials
%--------------------------------------------------------------------------
for  c = 1:size(X,1)

    % baseline parameters
    %----------------------------------------------------------------------
    Q  = P;

    % trial-specific inputs
    %----------------------------------------------------------------------
    for i = 1:size(X,2)
          Q.As = Q.As + X(c,i)*Q.Bs{i};
    end
    if isfield(Q,'Ac')
        for i = 1:size(X,2)
            Q.Ac = Q.Ac + X(c,i)*Q.Bc{i};
        end
    end

    if isfield(M,'trial')
        % Set up trial-specific initial state
        M.x=M.trial{c}.x;
    end

    x{c,1} = spm_int_ode(Q,M,U);

end

if isfield(M,'trial')
    M.x=Mx;
end

author: spm


text:
% dipfit_1_to_2() - convert dipfit 1 structure to dipfit 2 structure.
%
% Usage:
%  >> EEG.dipfit = dipfit_1_to_2(EEG.dipfit);
%
% Note:
% For non-standard BESA models (where the radii or the conductances
% have been modified, users must create a new model in Dipfit2 from
% the default BESA model.
%
% Author: Arnaud Delorme, SCCN, La Jolla 2005

% Copyright (C) Arnaud Delorme, SCCN, La Jolla 2005
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function newdipfit = dipfit_1_to_2( dipfit );

    if isfield( dipfit, 'model')
        newdipfit.model = dipfit.model;
    end;
    if isfield( dipfit, 'chansel')
        newdipfit.chansel = dipfit.chansel;
    end;

    ind = 1; % use first template (BESA)
    newdipfit.coordformat = template_models(ind).coordformat;
    newdipfit.mrifile     = template_models(ind).mrifile;
    newdipfit.chanfile    = template_models(ind).chanfile;

    if ~isfield(dipfit, 'vol')
        newdipfit.hdmfile = template_models(ind).hdmfile;
    else
        newdipfit.vol     = dipfit.vol;

        %if length(dipfit.vol) == 4
            %if ~all(dipfit.vol == [85-6-7-1 85-6-7 85-6 85]) | ...
            %   ~all(dipfit.c   == [0.33 1.00 0.0042 0.33]) | ...
            %   ~all(dipfit.o = [0 0 0])
            %    disp('Warning: Conversion from dipfit 1 to dipfit 2 can only deal');
            %    disp('         with standard (not modified) BESA model');
            %    disp('         See "help dipfit_1_to_2" to convert this model');
            %    newdipfit = [];
            %end;
        %end;
    end;


text:
function test_bug896

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_multiplotTFR ft_prepare_layout ft_datatype ft_datatype_freq

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug896.mat'));

ft_checkdata(stat_coh, 'datatype', 'freq');
ft_checkdata(stat_coh_full, 'datatype', 'freq');
ft_checkdata(stat_coh, 'datatype', 'freq', 'cmbstyle', 'full');
ft_checkdata(stat_coh_full, 'datatype', 'freq', 'cmbstyle', 'full');

if ~ft_datatype(stat_coh, 'freq')
  error('ft_datatype failed on the "stat_coh" input data');
end

if ~ft_datatype(stat_coh_full, 'freq')
  error('ft_datatype failed on the "stat_coh_full" input data');
end

figure
ft_multiplotTFR(cfg, stat_coh);
close

figure
ft_multiplotTFR(cfg, stat_coh_full);
close

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [W1,N] = UniformlyRandomlyPoint(N,M)
%UniformlyRandomlyPoint - Generate a set of uniform randomly distributed points on
%the unit hyperplane
%
%   [W,N] = UniformlyRandomlyPoint(N,M) returns N uniform randomly distributed
%   points with M objectives.
%
%   Example:
%       [W,N] = UniformlyRandomlyPoint(275,10)

%--------------------------------------------------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB Platform
% for Evolutionary Multi-Objective Optimization [Educational Forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Lucas Farias

    W1=eye(M,M);
	W1=[W1;ones(1,M)/M];

	W2=rand(5000,M);
	W2 = W2./repmat(sum(W2,2),1,size(W2,2));

	while size(W1,1) < N
		index = find_index_with_largest_distance (W1,W2);
		W1(size(W1,1)+1,:)=W2(index,:);
		W2(index,:)=[];
	end

    W1 = max(W1,1e-6);
    N = size(W1,1);

end

function index = find_index_with_largest_distance (W1,W2)
    Distance = pdist2(W2,W1);
    Temp     = sort(Distance,2);
    [~,Rank] = sortrows(Temp);
    index=Rank(length(Rank));
end


text:
function [PDR, deltaHD, deltaSEN, deltaPRO, deltaCOL, CBR] = CV2XMode4(beta,lambda,Pt,S,B);

% CV2XMode4 is the main script of the implementation of the analytical
% models of the communication performance of C-V2X or LTE-V Mode 4
% described in the following paper:
%
%    Manuel Gonzalez-Martín, Miguel Sepulcre, Rafael Molina-Masegosa, Javier Gozalvez,
%    "Analytical Models of the Performance of C-V2X Mode 4 Vehicular Communications",
%    IEEE Transactions on Vehicular Technology, Vol. 68, Issue 2, Feb. 2019. DOI: 10.1109/TVT.2018.2888704
%    Final version available at: https://ieeexplore.ieee.org/document/8581518
%    Post-print version available at: https://arxiv.org/abs/1807.06508
%
% The paper presents analytical models for the average PDR (Packet Delivery Ratio) as a
% function of the distance between transmitter and receiver, and for the four different
% types of transmission errors that can be encountered in C-V2X or LTE-V Mode 4. The models
% have been validated for a wide range of transmission parameters and traffic densities by
% comparing the results obtained with the analytical models to those obtained with a C-V2X
% or LTE-V Mode 4 simulator implemented by the authors over the Veins simulation platform.
%
% CV2XMode4.m is the main script you have to run to get the PDR curve as a function of the
% distance for a given set of parameters, and the probability of each of the four
% transmission errors.
%
% The resulting figures are compared with simulations when the same configuration
% is available in the ./simulations folder.
%
% The resulting figures are stored in the ./fig folder.
%
% Input parameters:
%    beta: traffic density in veh/m. Values tested: 0.1, 0.2 and 0.3.
%    lambda: packet transmission frequency in Hz. Values tested: 10 and 20.
%    Pt: transmission power in dBm. Values tested: 20 and 23.
%    S: number of sub-channels. Values tested: 2 and 4.
%    B: packet size in bytes. Values tested: 190.
%
% Output metrics:
%    PDR: Packet Delivery Ratio for different Tx-Rx distances
%    deltaHD: probability of packet loss due to half-duplex transmissions for different Tx-Rx distances
%    deltaSEN: probability of packet loss due to a received signal power below the sensing power threshold for different Tx-Rx distances
%    deltaPRO: probability of packet loss due to propagation effects for different Tx-Rx distances
%    deltaCOL: probability of packet loss due to packet collisions for different Tx-Rx distances
%    CBR: Channel Busy Ratio between 0 and 1
%
% Overall code structure:
%     CV2XMode4.m
%         |---->   CV2XMode4_common.m   ----> get_PL_SH.m, get_SINRdistribution.m, get_BLER.m
%         |---->   CV2XMode4_Step2.m    ----> get_PL_SH.m, get_SINRdistribution.m, get_BLER.m
%         |---->   CV2XMode4_Step3.m    ----> get_PL_SH.m, get_SINRdistribution.m, get_BLER.m
%
% The equations that are identified with a number between brackets in this script are the ones
% that also appear in the paper so that they can be easily identified.

    disp('=========================================================')
    disp('Input parameters:')
    fprintf('  beta   = %f veh/m \n', beta)
    fprintf('  lambda = %d Hz \n', lambda)
    fprintf('  Pt     = %d dBm \n', Pt)
    fprintf('  S      = %d subchannels \n', S)
    fprintf('  B      = %d bytes \n', B)

    distance = [0:25:500];  % Tx-Rx distances to evaluate (m)

    Psen = -90.5;               % Sensing threshold (dBm)

    step_dB = 0.1;              % Discrete steps to compute the PDF of the SNR and SINR (dB)

    % Calculate the number of RBs that are needed to transmit each message
    % and the coding used based on the number of sub-channels and packet size:
    switch B
        case 190
            switch S
                case 4
                    coding = 1;   % Used to identify the BLER vs SINR curve to be used (190 Bytes, QPSK r=0.7, Vr = 280 km/h)
                    RBs = 10;     % Number of RBs needed to transmit the DATA field of each message
                case 2
                    coding = 2;   % Used to identify the BLER vs SINR curve to be used (190 Bytes, QPSK r=0.5, Vr = 280 km/h)
                    RBs = 12;     % Number of RBs needed to transmit the DATA field of each message
            end
    end

    noise = -95 - 10*log10(50/RBs);     % Noise corresponding to the DATA field of each message. Assumes a noise figure of 9dB and 10MHz channel (background noise of -95dBm). The total number of RBs in 10MHz is 50.

    % Calculate errors associated to HD, SEN and PRO:
    [ deltaHD_pre , deltaSEN_pre , deltaPRO_pre ] = CV2XMode4_common( lambda , Pt , distance, Psen , step_dB , noise , coding );

    % Calculate probability of collision considering only Step 2 and CBR:
    [ deltaCOL2_pre , CBR ] = CV2XMode4_Step2( beta , lambda , Pt , S , distance , Psen , step_dB , noise , coding , deltaPRO_pre );

    % Calculate weighting factor alpha using equation (22):
    if CBR < 0.2
        alpha = 0;
    elseif CBR <= 0.7
        alpha = 2*CBR - 0.4;
    else
        alpha = 1;
    end

    % Calculate probability of collision considering only Step 3:
    if alpha < 1
         [ deltaCOL3_pre ] = CV2XMode4_Step3( beta , lambda , Pt , S , distance , Psen , step_dB , noise , coding , deltaPRO_pre );
    else
        deltaCOL3_pre = 0;
    end

    % Calculate final probabilities for each type of error:
    deltaHD   = deltaHD_pre;                               % Equation (6.1)
    deltaSEN  = deltaSEN_pre .* (1 - deltaHD);             % Equation (6.2)
    deltaPRO  = deltaPRO_pre .* (1 - deltaHD_pre) .* (1 - deltaSEN_pre);  % Equation (6.3)
    deltaCOL2 = deltaCOL2_pre .* (1 - deltaHD_pre) .* (1 - deltaSEN_pre) .* (1 - deltaPRO_pre); % Equation (6.4)
    deltaCOL3 = deltaCOL3_pre .* (1 - deltaHD_pre) .* (1 - deltaSEN_pre) .* (1 - deltaPRO_pre); % Equation (6.5)
    deltaCOL = alpha*deltaCOL2 + (1-alpha)*deltaCOL3; % Equation (21)

    % Calculate PDR:
    PDR = 1 - deltaHD - deltaSEN - deltaPRO - deltaCOL; % Equation (6)

    % Presentation of the obtained results:

    % Load simulation results (if available):
    simulation_path = [pwd '\simulations\' num2str(S) 'subchannels\'  num2str(Pt) 'dBm'];
    sim_file = [ simulation_path '\FALLOS_LOS_D2D_' num2str(beta*1000) 'vehpkm_16Alg_' num2str(1/lambda) 's_' num2str(8*B) 'bit.fig' ];
    fig_name = ['CV2XMode4_beta' num2str(beta) '-lambda' num2str(lambda) '-Pt' num2str(Pt) '-S' num2str(S) '-B' num2str(B)];

    if exist(sim_file,'file')==2
        open(sim_file)
        lh = findall(gca, 'type', 'line');
        X = get(lh,'xdata');
        Y = get(lh,'ydata');
        deltaHD_sim = Y{5};
        deltaSEN_sim = Y{1};
        deltaPRO_sim = Y{4};
        deltaCOL_sim = Y{3};
        close

        figure; hold on; grid on; box on
        plot(distance , deltaHD_sim/100,'b-','LineWidth',2)
        plot(distance , deltaSEN_sim/100,'m-','LineWidth',2)
        plot(distance , deltaPRO_sim/100,'r-','LineWidth',2)
        plot(distance , deltaCOL_sim/100,'k-','LineWidth',2)

    else
        figure; hold on; grid on; box on
    end

    % Plot errors:
    ylabel('Error probability')
    xlabel('Distance [m]')
    plot(distance , deltaHD,'b--','LineWidth',2)
    plot(distance , deltaSEN,'m--','LineWidth',2)
    plot(distance , deltaPRO,'r--','LineWidth',2)
    plot(distance , deltaCOL,'k--','LineWidth',2)
    ylim([0 1])

    if exist(sim_file,'file')==2
        legend('\delta_{HD} Simulation','\delta_{SEN} Simulation','\delta_{PRO} Simulation','\delta_{COL} Simulation','\delta_{HD} Analytical','\delta_{SEN} Analytical','\delta_{PRO} Analytical','\delta_{COL} Analytical','Location','northwest')
    else
        legend('\delta_{HD} Analytical','\delta_{SEN} Analytical','\delta_{PRO} Analytical','\delta_{COL} Analytical','Location','northwest')
    end
    hgsave(['fig/' fig_name '_errors.fig'])

    % Plot PDR:
    figure; hold on; grid on; box on
    if exist(sim_file,'file')==2
        PDR_sim = (100 - deltaHD_sim - deltaSEN_sim - deltaPRO_sim - deltaCOL_sim)/100;
        plot(distance,PDR_sim,'b-','LineWidth',2)
        plot(distance,PDR,'b--','LineWidth',2)
        legend('PDR Simulation','PDR Analytical')
    else
        plot(distance,PDR,'--')
        legend('PDR Analytical')
    end
    ylim([0 1])
    ylabel('PDR')
    xlabel('Distance [m]')
    hgsave(['fig/' fig_name '_PDR.fig'])

    fprintf('Channel load: CBR = %.2f, alpha = %f \n',CBR, alpha)

    if exist(sim_file,'file')==2
        MAD_HD  = mean( abs(deltaHD*100 - deltaHD_sim) );    % Equation (35)
        MAD_SEN = mean( abs(deltaSEN*100 - deltaSEN_sim) );  % Equation (35)
        MAD_PRO = mean( abs(deltaPRO*100 - deltaPRO_sim) );  % Equation (35)
        MAD_COL = mean( abs(deltaCOL*100 - deltaCOL_sim) );  % Equation (35)
        MAD     = mean( abs(PDR - PDR_sim) * 100 );          % Equation (35)

        disp('Mean Absolute Deviation results : ')
        fprintf('PDR \tHD  \tSEN \tPRO \tCOL \n')
        fprintf('%.2f\t%.2f\t%.2f\t%.2f\t%.2f \n', MAD, MAD_HD, MAD_SEN, MAD_PRO, MAD_COL)
    end

    disp('=========================================================')

 return

author: msepulcre


text:
%--- help for embed ---
%
%  INTERNAL FUNCTION: Maps the data into proper lag/stacked form
%
%

author: jmaih


text:
function [ n_data, n, a, b, x, fx ] = j_polynomial_values ( n_data )

%*****************************************************************************80
%
%% J_POLYNOMIAL_VALUES returns some values of the Jacobi polynomial J(n,a,b,x).
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      JacobiP[ n, a, b, x ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the degree of the polynomial.
%
%    Output, real A, B, parameters of the function.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 26;

  a_vec = [ ...
     0.0, 0.0, 0.0, 0.0, ...
     0.0, 0.0, 1.0, 2.0, ...
     3.0, 4.0, 5.0, 0.0, ...
     0.0, 0.0, 0.0, 0.0, ...
     0.0, 0.0, 0.0, 0.0, ...
     0.0, 0.0, 0.0, 0.0, ...
     0.0, 0.0 ];

  b_vec = [ ...
    1.0, 1.0, 1.0, 1.0, ...
    1.0, 1.0, 1.0, 1.0, ...
    1.0, 1.0, 1.0, 2.0, ...
    3.0, 4.0, 5.0, 1.0, ...
    1.0, 1.0, 1.0, 1.0, ...
    1.0, 1.0, 1.0, 1.0, ...
    1.0, 1.0 ];

  fx_vec = [ ...
      1.000000000000000, ...
      0.2500000000000000, ...
     -0.3750000000000000, ...
     -0.4843750000000000, ...
     -0.1328125000000000, ...
      0.2753906250000000, ...
     -0.1640625000000000, ...
     -1.174804687500000, ...
     -2.361328125000000, ...
     -2.616210937500000, ...
      0.1171875000000000, ...
      0.4218750000000000, ...
      0.5048828125000000, ...
      0.5097656250000000, ...
      0.4306640625000000, ...
     -6.000000000000000, ...
      0.03862000000000000, ...
      0.8118400000000000, ...
      0.03666000000000000, ...
     -0.4851200000000000, ...
     -0.3125000000000000, ...
      0.1891200000000000, ...
      0.4023400000000000, ...
      0.01216000000000000, ...
     -0.4396200000000000, ...
      1.000000000000000 ];

  n_vec = [ ...
     0, 1, 2, 3, ...
     4, 5, 5, 5, ...
     5, 5, 5, 5, ...
     5, 5, 5, 5, ...
     5, 5, 5, 5, ...
     5, 5, 5, 5, ...
     5, 5 ];

  x_vec = [ ...
      0.5, ...
      0.5, ...
      0.5, ...
      0.5, ...
      0.5, ...
      0.5, ...
      0.5, ...
      0.5, ...
      0.5, ...
      0.5, ...
      0.5, ...
      0.5, ...
      0.5, ...
      0.5, ...
      0.5, ...
     -1.0, ...
     -0.8, ...
     -0.6, ...
     -0.4, ...
     -0.2, ...
      0.0, ...
      0.2, ...
      0.4, ...
      0.6, ...
      0.8, ...
      1.0 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    a = 0.0;
    b = 0.0;
    x = 0.0;
    fx = 0.0;
  else
    n = n_vec(n_data);
    a = a_vec(n_data);
    b = b_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ixZ_vec, endleaves_vec, minBranche] = recursivesubtreePlus(Z, ixZ, ixZ_vec, endleaves_vec, minBranche)
% [ixZ_vec, endleaves_vec] = recursivesubtree(Z, ixZ, ixZ_vec, endleaves_vec)
% recursively finds all subnodes and leaves from the node with index ixZ in
% matrix Z
% initialized as:
% [ixZ_vec, endleaves_vec] = recursivesubtree(Z, ixZ,[],[])
% Z output from 'linkage' function
% ixZ -  index of the node in matrix Z (eg ixZ=(find(Z(:,3)==distance_value)))
% Z = linkage(ccds,'complete');
% [H,T,perm] = dendrogram(Z,0, 'colorthreshold',0.3);
% then subtree of node with ixZ = 50 in dendrogram
% [H,T,perm] = dendrogram(Z,0);
% can be enhanced by eg
% set(H(recursivesubtree(Z,50,[],[])),'LineWidth',2)
ixZ_vec(end+1) = ixZ;
value = Z(ixZ, 1:2) - (length(Z)+1);
if and(value(1)>0, value(1)>0)
    if isempty(minBranche)
        minBranche = ixZ; %to initialize
    end
    if Z(ixZ,3)<Z(minBranche,3)
        minBranche = ixZ;
    end
end
if value(1) > 0
    ixZ = value(1);
    ixZ_vec(end+1) = ixZ;
    [ixZ_vec, endleaves_vec, minBranche ] = recursivesubtreePlus(Z, ixZ, ixZ_vec, endleaves_vec, minBranche );
else
    endleaves_vec(end+1) = Z(ixZ, 1);
end
if value(2) > 0
    ixZ = value(2);
    ixZ_vec(end+1) = ixZ;
    [ixZ_vec, endleaves_vec, minBranche ] = recursivesubtreePlus(Z, ixZ, ixZ_vec, endleaves_vec, minBranche );
else
    endleaves_vec(end+1) = Z(ixZ, 2);
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function out = disp2str(x)
%DISP2STR     Output of intval x into out
%
%The output produced by disp_(x(:)), infsup(x(:)), or midrad(x(:)),
%  according to the format in use, is returned into  out  as follows:
%
% out.exp  empty if no common exponent printed, otherwise
%            string of common exponent
% out.str  column array of strings representing to x(:)
%

% written  08/29/00     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 08/26/12     S.M. Rump  global variables removed
%

  INTLAB_INTVAL_DISPLAY = getappdata(0,'INTLAB_INTVAL_DISPLAY');
%VVVV x = x(:);
  x = reshape(x,prod(size(x)),1);
%AAAA Matlab V5.2 bug fix
  switch INTLAB_INTVAL_DISPLAY
    case 'DisplayInfsup', out = infsup(x,[],[]);
    case 'DisplayMidrad', out = midrad(x,[],[]);
    case 'Display_', out = disp_(x,[],[]);
  end

repo: OpenMAS


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

% Edit the above text to modify the response to help avg_data_quality

% Last Modified by GUIDE v2.5 17-Jan-2023 13:20:16

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @avg_data_quality_OpeningFcn, ...
    'gui_OutputFcn',  @avg_data_quality_OutputFcn, ...
    'gui_LayoutFcn',  [] , ...
    'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

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

handles.output = hObject;

disp('Awaiting Data Quality settings...');

% Parse input
handles.DQ_spec = varargin{1};
handles.timelimits = varargin{2};

if isempty(varargin{1})

    preAvg_info = make_DQ_spec(handles.timelimits);
    sme_tw_labels = preAvg_info(3).time_window_labels;
    sme_tw = num2cell(preAvg_info(3).times);
    preAvg_labels = strrep(sme_tw_labels,'aSME','(Metric)');
    sme_tw3= [preAvg_labels' sme_tw];

else
    preAvg_info = varargin{1};

    nDQ = length(preAvg_info);
    times_ind = zeros(1, nDQ);

    %search "times" struct
    %time_present = 0;
    for ti = 1:nDQ

        if ~isempty(preAvg_info(ti).times)
            times_ind(ti) = 1;
        end

    end

    if ~any(times_ind) %no valid times
        sme_tw3 = [];
%         preAvg_labels_old = preAvg_info(3).time_window_labels';
%         preAvg_info_times = num2cell(preAvg_info(3).times);
    else
         try %in case there is only one time-window metric
             preAvg_info_times = num2cell(preAvg_info.times);
             %             preAvg_info_times = num2cell(preAvg_info.times);
             %

%              sme_tw2 = cellstr(strcat({'(Metric) at '}, ...
%                 string(preAvg_info_times(:,1)),{' to '},string(preAvg_info_times(:,2))));
            sme_tw2 = preAvg_info.time_window_labels;

            sme_tw2 = strrep(sme_tw2,'aSME','');
            sme_tw2 = strrep(sme_tw2,'(Corrected)','');
            sme_tw2 = strrep(sme_tw2,'aSD','');
            sme_tw2 = strcat('(Metric)',sme_tw2);

            try
                sme_tw3 = [sme_tw2 preAvg_info_times];
            catch
                sme_tw3 = [sme_tw2' preAvg_info_times];
            end

         catch
             %find usable inddex
             use_row = find(times_ind ==1);
             preAvg_info_times = num2cell(preAvg_info(use_row(1)).times);

             %retain time_window_labels if previously set

%             sme_tw2 = cellstr(strcat({'(Metric) at '}, ...
%                 string(preAvg_info_times(:,1)),{' to '},string(preAvg_info_times(:,2))));

            sme_tw2 = preAvg_info(use_row(1)).time_window_labels;
            sme_tw2 = strrep(sme_tw2,'aSME','');
            sme_tw2 = strrep(sme_tw2,'(Corrected)','');
            sme_tw2 = strrep(sme_tw2,'aSD','');
          %  sme_tw2 = strrep(sme_tw2,'aSD (Corrected)','(');

            sme_tw2 = strcat('(Metric)',sme_tw2);

            try
                sme_tw3 = [sme_tw2 preAvg_info_times];
            catch
                sme_tw3 = [sme_tw2' preAvg_info_times];
            end
         end

    end

%     try
%
%     catch ME
%         switch ME.identifier
%             case 'MATLAB:nonExistentField'
%                 preAvg_labels_old = 'EMPTY';
%                 preAvg_info_times = 'NAN';
%             otherwise
%                 preAvg_labels_old = preAvg_info.time_window_labels';
%                 preAvg_info_times = num2cell(preAvg_info.times);
%         end
%
%     end
%     preAvg_labels = strrep(preAvg_labels_old,'aSME','(Metric)');
%
%
%     try
%         sme_tw2= [preAvg_labels' preAvg_info_times];
%     catch
%         sme_tw2= [preAvg_labels preAvg_info_times];
%     end

end

% Choose default command line output for avg_data_quality

handles.DQout = zeros(2,4);
handles.tupdate = 0;
handles.num_tests = 1;

handles.paraSME = 1;
handles.tdata = handles.SME_table.Data;
handles.SME_table.Data = sme_tw3;

handles.Tout = sme_tw3;
disp(handles.Tout);

handles.sel_row = size(handles.Tout,1); % set to max on load

%read from default/working memory values
type_struct = {preAvg_info.type};

for curtype = type_struct
    switch char(curtype)

        case 'Baseline Measure - SD'
            set(handles.checkbox_baseline_noise,'Value',1);
            set(handles.radiobutton_basel_sd,'Value',1);

        case 'Baseline Measure - SD (Corrected)'
            set(handles.checkbox_baseline_noise,'Value',1);
            set(handles.radiobutton_basel_sdcorr,'Value',1);

        case 'Baseline Measure - RMS'
            set(handles.checkbox_baseline_noise,'Value',1);
            set(handles.radiobutton_basel_rms,'Value',1);

        case 'Point-wise SEM'
            set(handles.checkbox_SEM,'Value',1);
            set(handles.radiobutton12, 'Value',1);

        case 'Point-wise SEM (Corrected)'
            set(handles.checkbox_SEM,'Value',1);
            set(handles.radiobutton_sem_corr,'Value',1);

        case 'aSME'
            set(handles.checkbox1_paraSME,'Value',1);
            set(handles.radiobutton_SD,'Value',1);

        case 'aSME (Corrected)'
            set(handles.checkbox1_paraSME,'Value',1);
            set(handles.radiobutton_SD_corr,'Value',1);

        case 'aSD'
            set(handles.checkbox6_paraSD,'Value',1);
            set(handles.radiobutton_SD,'Value',1);

        case 'aSD (Corrected)'
            set(handles.checkbox6_paraSD,'Value',1);
            set(handles.radiobutton_SD_corr,'Value',1);
    end

end

% Update handles structure
guidata(hObject, handles);

%helpbutton
initialize_gui(hObject, handles, false);

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

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

% Get default command line output from handles structure
pause(0.1)
% Update handles structure
guidata(hObject, handles);

% Get default command line output from handles structure
varargout{1} = handles.output;
% The figure can be deleted now
pause(0.1)
delete(handles.avg_dq);
pause(0.1)

% --- Executes on button press in cancel.
function cancel_Callback(hObject, eventdata, handles)
% hObject    handle to cancel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isequal(get(handles.avg_dq, 'waitstatus'), 'waiting')
    %The GUI is still in UIWAIT, us UIRESUME
    handles.output = '';
    %Update handles structure
    guidata(hObject, handles);
    uiresume(handles.avg_dq);
else
    % The GUI is no longer waiting, just close it
    delete(handles.avg_dq);
end

% --- Executes on button press in save.
function save_Callback(hObject, eventdata, handles)
% hObject    handle to save (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% First, check time windows are valid
try
    times_here = cell2mat(handles.Tout(:,2:3));
    if any(times_here(:) < min(handles.timelimits)) || any(times_here(:) > max(handles.timelimits))
        beep
        time_off_str = ['It appears that some of the listed times are out of bounds of the epoch time of this EEG set'];
        disp(time_off_str)
        time_off2 = ['Time limits here are: ' num2str(handles.timelimits)];
        disp(time_off2)
        pause(0.1)
        return
    end
catch
    times_here = [];

end

disp('Saving Data Quality settings...');

if handles.paraSME == 1
    num_tests = 1;
    type = 'SME';
end

baseline_on = get(handles.checkbox_baseline_noise,'Value');
sem_on = get(handles.checkbox_SEM,'Value');
asme_on = get(handles.checkbox1_paraSME,'Value');
aSD_on = get(handles.checkbox6_paraSD,'Value');

%num_tests = baseline_on + sem_on + asme_on;

% Make DQout
DQout = [];
dq_slot = 0;
if baseline_on
    dq_slot = dq_slot + 1;

    basel_sd = get(handles.radiobutton_basel_sd,'Value');
    basel_rms = get(handles.radiobutton_basel_rms,'Value');
    %basel_sdcorr = get(handles.radiobutton_basel_sdcorr,(
    if basel_sd
        DQout(dq_slot).type = 'Baseline Measure - SD';
    elseif basel_rms
        DQout(dq_slot).type = 'Baseline Measure - RMS';
    else
        DQout(dq_slot).type = 'Baseline Measure - SD (Corrected)';
    end

    default_t = get(handles.checkbox_basel_default_times,'Value');
    if default_t
        DQout(dq_slot).times = [];
    else
        custom_t = [1 str2double(get(handles.baseline_start,'String')) str2double(get(handles.baseline_end,'String'))];
        DQout(dq_slot).times = custom_t;
    end
end

if sem_on
    dq_slot = dq_slot + 1;
    sem_corr = get(handles.radiobutton_sem_corr,'Value');

    if sem_corr
        DQout(dq_slot).type = 'Point-wise SEM (Corrected)';
    else
        DQout(dq_slot).type = 'Point-wise SEM';
    end
end

if asme_on
    dq_slot = dq_slot + 1;

    dq_sd = get(handles.radiobutton_SD,'Value');
    if dq_sd
        DQout(dq_slot).type = 'aSME';
        DQout(dq_slot).times = times_here;
        DQout(dq_slot).time_window_labels = strrep(handles.Tout(:,1),'(Metric)','aSME');
    else
        DQout(dq_slot).type = 'aSME (Corrected)';
        DQout(dq_slot).times = times_here;
        DQout(dq_slot).time_window_labels = strrep(handles.Tout(:,1),'(Metric)','aSME (Corrected)');
    end

   % DQout(dq_slot).time_window_labels = strrep(handles.Tout(:,1),'(Metric)','aSME');

end

if aSD_on
    dq_slot = dq_slot + 1;

    dq_sd = get(handles.radiobutton_SD,'Value');
    if dq_sd
        DQout(dq_slot).type = 'aSD';
        DQout(dq_slot).times = times_here;
        DQout(dq_slot).time_window_labels = strrep(handles.Tout(:,1),'(Metric)','aSD');
    else
        DQout(dq_slot).type = 'aSD (Corrected)';
        DQout(dq_slot).times = times_here;
        DQout(dq_slot).time_window_labels = strrep(handles.Tout(:,1),'(Metric)','aSD (Corrected)');
    end

   % DQout(dq_slot).time_window_labels = strrep(handles.Tout(:,1),'(Metric)','aSD');

end

%DQout.num_tests = num_tests;

handles.output = DQout;

% Update handles structure
guidata(hObject, handles);
uiresume(handles.avg_dq);
%avg_dq_CloseRequestFcn(hObject, eventdata, handles);

% --- Executes when selected object changed in unitgroup.
function unitgroup_SelectionChangedFcn(hObject, eventdata, handles)
% hObject    handle to the selected object in unitgroup
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function initialize_gui(fig_handle, handles, isreset)
% If the metricdata field is present and the save flag is false, it means
% we are we are just re-initializing a GUI by calling it from the cmd line
% while it is up. So, bail out as we dont want to save the data.
if isfield(handles, 'metricdata') && ~isreset
    return;
end

% Update handles structure
guidata(handles.avg_dq, handles);

% --- Executes on button press in checkbox1_paraSME.
function checkbox1_paraSME_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox1_paraSME (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox1_paraSME
SME_here = get(hObject,'Value');
%disp(SME_here);

temp_times = handles.Tout;

if SME_here
    set(handles.SME_table,'Enable','on');
    set(handles.radiobutton_SD,'Enable','on');
    set(handles.radiobutton_SD_corr,'Enable','on');

    if isempty(temp_times)
        temp_DQ_spec = make_DQ_spec(handles.timelimits);
        sme_tw_labels_old = temp_DQ_spec(3).time_window_labels;
        sme_tw_labels = strrep(sme_tw_labels_old,'aSME','(Metric)');
        sme_tw = num2cell(temp_DQ_spec(3).times);
        sme_tw2 = [sme_tw_labels' sme_tw];

        handles.Tout = sme_tw2;
        set(handles.SME_table, 'Data', sme_tw2);
        pause(0.3)
    end

else

    SD_check = get(handles.checkbox6_paraSD,'Value');
    if SD_check == 0
        set(handles.SME_table,'Enable','off');
        set(handles.radiobutton_SD,'Enable','off');
        set(handles.radiobutton_SD_corr,'Enable','off');

    end
end

guidata(hObject, handles);

function pushbutton_help_Callback(hObject, eventdata, handles)
% doc pop_averager
web('https://github.com/lucklab/erplab/wiki/Computing-Averaged-ERPs#data-quality-measures', '-browser');

% --- Executes when user attempts to close avg_dq.
function avg_dq_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to avg_dq (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isequal(get(handles.avg_dq, 'waitstatus'), 'waiting')
    %The GUI is still in UIWAIT, us UIRESUME
    handles.output = '';
    %Update handles structure
    guidata(hObject, handles);
    uiresume(handles.avg_dq);
else
    % The GUI is no longer waiting, just close it
    delete(handles.avg_dq);
end
% Hint: delete(hObject) closes the figure
%delete(hObject);

% --- Executes during object creation, after setting all properties.
function SME_table_CreateFcn(hObject, eventdata, handles)
% hObject    handle to SME_table (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

Tout = hObject.Data;
handles.Tout = Tout;
guidata(hObject, handles);

% --- Executes when entered data in editable cell(s) in SME_table.
function SME_table_CellEditCallback(hObject, eventdata, handles)
% hObject    handle to SME_table (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.CONTROL.TABLE)
%	Indices: row and column indices of the cell(s) edited
%	PreviousData: previous data for the cell(s) edited
%	EditData: string(s) entered by the user
%	NewData: EditData or its converted form set on the Data property. Empty if Data was not changed
%	Error: error string when failed to convert EditData to appropriate value for Data
% handles    structure with handles and user data (see GUIDATA)
% handles.tdata = eventdata.NewData;
% handles.SME_table.Data = handles.tdata;

Tout = hObject.Data;
handles.Tout = Tout;
guidata(hObject, handles);

% --- Executes on button press in SME_row_minus.
function SME_row_minus_Callback(hObject, eventdata, handles)
% hObject    handle to SME_row_minus (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

curr_rows = size(handles.Tout,1);

row_del = handles.sel_row;

if curr_rows <= 1
    beep
    disp('Already at 1 rows')
else

    % notify
    row_del_text = ['Now removing row ' num2str(row_del)];
    disp(row_del_text)

    new_rows = curr_rows - 1;

    new_Tout = handles.Tout;
    new_Tout(row_del,:) = []; % pop the selected row out
    handles.Tout = new_Tout;

    %     row_drop_txt = ['The number of rows was ' num2str(curr_rows) '. Now dropping to ' num2str(new_rows)];
    % disp(row_drop_txt)
    set(handles.SME_table,'Data',new_Tout)
    pause(0.3)
    handles.sel_row = new_rows;

    %

end

guidata(hObject, handles);

% --- Executes on button press in SME_row_plus.
function SME_row_plus_Callback(hObject, eventdata, handles)
% hObject    handle to SME_row_plus (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

curr_rows = size(handles.Tout,1);
new_rows = curr_rows + 1;

old_Tout = handles.Tout;
win_size = old_Tout{end,3} - old_Tout{end,2};

new_row_str = ['(Metric) Time Window custom ' num2str(new_rows)];
new_row_cell = {new_row_str,old_Tout{end,3},old_Tout{end,3}+win_size};
new_Tout = [old_Tout;new_row_cell];

handles.Tout = new_Tout;

set(handles.SME_table,'Data',new_Tout);

guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function text_time_range_CreateFcn(hObject, eventdata, handles)
% hObject    handle to text_time_range (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% --- Executes on button press in checkbox_baseline_noise.
function checkbox_baseline_noise_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox_baseline_noise (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox_baseline_noise
inc_baseline = get(hObject,'Value');
if inc_baseline == 1
    set(handles.checkbox_basel_default_times,'Value',1);
    set(handles.radiobutton_basel_custom_times,'Value',1);
    set(handles.radiobutton_basel_sd,'Value',1);
    set(handles.radiobutton_basel_sdcorr,'Value',1);
    set(handles.radiobutton_basel_rms,'Value',1);
else
    set(handles.checkbox_basel_default_times,'Value',0);
    set(handles.radiobutton_basel_custom_times,'Value',0);
    set(handles.radiobutton_basel_sd,'Value',0);
    set(handles.radiobutton_basel_sdcorr,'Value',0);
    set(handles.radiobutton_basel_rms,'Value',0);
end
guidata(hObject, handles);

% --- Executes on button press in checkbox_SEM.
function checkbox_SEM_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox_SEM (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox_SEM
inc_SEM = get(hObject,'Value');
if inc_SEM == 1
    set(handles.radiobutton12,'Value',1);
    set(handles.radiobutton_sem_corr,'Value',1);
else
    set(handles.radiobutton12,'Value',0);
    set(handles.radiobutton_sem_corr,'Value',0);
end
guidata(hObject, handles);

function baseline_start_Callback(hObject, eventdata, handles)
% hObject    handle to baseline_start (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of baseline_start as text
%        str2double(get(hObject,'String')) returns contents of baseline_start as a double

% --- Executes during object creation, after setting all properties.
function baseline_start_CreateFcn(hObject, eventdata, handles)
% hObject    handle to baseline_start (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function baseline_end_Callback(hObject, eventdata, handles)
% hObject    handle to baseline_end (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of baseline_end as text
%        str2double(get(hObject,'String')) returns contents of baseline_end as a double

% --- Executes during object creation, after setting all properties.
function baseline_end_CreateFcn(hObject, eventdata, handles)
% hObject    handle to baseline_end (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in checkbox_baseline_subtract.
function checkbox_baseline_subtract_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox_baseline_subtract (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox_baseline_subtract

% --- Executes on button press in radiobutton_basel_sd.
function radiobutton_basel_sd_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton_basel_sd (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton_basel_sd
basel_sd = get(hObject,'Value');

if basel_sd
    set(handles.radiobutton_basel_rms,'Value',0);
end

% --- Executes on button press in radiobutton_basel_rms.
function radiobutton_basel_rms_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton_basel_rms (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton_basel_rms
basel_rms = get(hObject,'Value');

if basel_rms
    set(handles.radiobutton_basel_sd,'Value',0);
end

% --- Executes on button press in checkbox_basel_default_times.
function checkbox_basel_default_times_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox_basel_default_times (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox_basel_default_times
default_times = get(hObject,'Value');

if default_times
    set(handles.baseline_start,'Enable','off');
    set(handles.baseline_end,'Enable','off');
    %set(handles.radiobutton_basel_custom_times,'Enable','off');
    set(handles.text_basel2,'Enable','off');
    set(handles.text_basel3,'Enable','off');

else
    set(handles.baseline_start,'Enable','on');
    set(handles.baseline_end,'Enable','on');
    %set(handles.radiobutton_basel_custom_times,'Enable','on');
    set(handles.text_basel2,'Enable','on');
    set(handles.text_basel3,'Enable','on');
end

% --- Executes on button press in radiobutton_basel_custom_times.
function radiobutton_basel_custom_times_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton_basel_custom_times (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton_basel_custom_times
custom_times = get(hObject,'Value');

if custom_times
    set(handles.baseline_start,'Enable','on');
    set(handles.baseline_end,'Enable','on');
    %set(handles.radiobutton_basel_custom_times,'Enable','on');
    set(handles.text_basel2,'Enable','on');
    set(handles.text_basel3,'Enable','on');
else
    set(handles.baseline_start,'Enable','off');
    set(handles.baseline_end,'Enable','off');
    %set(handles.radiobutton_basel_custom_times,'Enable','off');
    set(handles.text_basel2,'Enable','off');
    set(handles.text_basel3,'Enable','off');
end

% --- Executes when selected cell(s) is changed in SME_table.
function SME_table_CellSelectionCallback(hObject, eventdata, handles)
% hObject    handle to SME_table (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.CONTROL.TABLE)
%	Indices: row and column indices of the cell(s) currently selecteds
% handles    structure with handles and user data (see GUIDATA)

%disp(eventdata)
% if the selected cell info exists, save it to handles
if numel(eventdata.Indices)
    row_here = eventdata.Indices(1);
    old_row = handles.sel_row;
    handles.sel_row = row_here;
    if isequal(row_here,old_row) == 0
        row_text = ['Row number ' num2str(row_here) ' now selected'];
        disp(row_text)
    end
    guidata(hObject, handles);
end

% --- Executes on button press in restbutton.
function restbutton_Callback(hObject, eventdata, handles)
% hObject    handle to restbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% get default timewindows
temp_DQ_spec = make_DQ_spec(handles.timelimits);
sme_tw_labels_old = temp_DQ_spec(3).time_window_labels;
sme_tw_labels = strrep(sme_tw_labels_old,'aSME','(Metric)');
sme_tw = num2cell(temp_DQ_spec(3).times);
sme_tw2 = [sme_tw_labels' sme_tw];

handles.Tout = sme_tw2;
set(handles.SME_table, 'Data', sme_tw2);
pause(0.3)

guidata(hObject, handles);

% --- Executes on button press in checkbox7.
function checkbox7_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox7 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox7

% --- Executes on button press in pushbutton9.
function pushbutton9_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton9 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in pushbutton10.
function pushbutton10_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton10 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in pushbutton11.
function pushbutton11_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton11 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in checkbox6_paraSD.
function checkbox6_paraSD_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox6_paraSD (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
SD_check = get(hObject,'Value');
%disp(SME_here);
temp_times = handles.Tout;

if SD_check
    set(handles.SME_table,'Enable','on');
    set(handles.radiobutton_SD,'Enable','on');
    set(handles.radiobutton_SD_corr,'Enable','on');

    if isempty(temp_times)
        temp_DQ_spec = make_DQ_spec(handles.timelimits);
        sme_tw_labels_old = temp_DQ_spec(3).time_window_labels;
        sme_tw_labels = strrep(sme_tw_labels_old,'aSME','(Metric)');
        sme_tw = num2cell(temp_DQ_spec(3).times);
        sme_tw2 = [sme_tw_labels' sme_tw];

        handles.Tout = sme_tw2;
        set(handles.SME_table, 'Data', sme_tw2);
        pause(0.3)
    end

else
    SME_check = get(handles.checkbox1_paraSME,'Value');
    if SME_check == 0
    set(handles.SME_table,'Enable','off');
    set(handles.radiobutton_SD,'Enable','off');
    set(handles.radiobutton_SD_corr,'Enable','off');
    end
end

guidata(hObject, handles);

% Hint: get(hObject,'Value') returns toggle state of checkbox6_paraSD

% --- Executes on button press in radiobutton_SD_corr.
function radiobutton_SD_corr_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton_SD_corr (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton_SD_corr

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
%% Data Structure: Boundary Conditions
%
% We use |bdFlag(1:NT,1:3)| to record the type of three edges of each
% triangle. Similarly in 3-D, we use |bdFlag(1:NT,1:4)| to record the type
% of four faces of each tetrahedron. The value is the type of boundary
% condition.
%
% * 0: non-boundary, i.e., an interior edge or face.
% * 1: first type, i.e., a Dirichlet boundary edge or face.
% * 2: second type, i.e., a Neumann boundary edge or face.
% * 3: third type, i.e., a Robin boundary edge or face.

%% Local labeling of edges and faces
% We label three edges of a triangle such that |bdFlag(t,i)| is the edge
% opposite to the i-th vertex. Similarly |bdFlag(t,i)| is the face opposite
% to the i-th vertex.

node = [1,0; 1,1; 0,0];
elem = [1 2 3];
edge = [2 3; 1 3; 1 2];
showmesh(node,elem);
findnode(node);
findedge(node,edge);

%% Set up boundary conditions
%
% The function |setboundary| is to set up the bdFlag matrix for a 2-D
% triangulation and |setboundary3| for a 3-D triangulation.
%
help setboundary

%%
% Note that if the i-th edge of t is on the boundary but |bdFlag(t,i)=0|,
% it is equivalent to use homogenous Neumann boundary condition (zero
% flux).

%% Example: Crack Domain
node = [1,0; 0,1; -1,0; 0,-1; 0,0; 1,0];        % nodes
elem = [5,1,2; 5,2,3; 5,3,4; 5,4,6];            % elements
elem = label(node,elem);                        % label the mesh
figure;
showmesh(node,elem);                            % plot mesh
findelem(node,elem);                            % plot element indices
findnode(node,2:6);                             % plot node indices
text(node(6,1),node(6,2)+0.075,int2str(1),'FontSize',16,'FontWeight','bold');
hold on;
plot([node(1,1),node(5,1)], [node(1,2),node(5,2)],'r-', 'LineWidth',3);
bdFlag = setboundary(node,elem,'Dirichlet');                % Dirichlet boundary condition
display(elem)
display(bdFlag)
%%
% node 1 and node 6 are the same point (1,0)

%%
bdFlag = setboundary(node,elem,'Dirichlet','abs(x) + abs(y) == 1','Neumann','y==0');
display(bdFlag)

%% Example: Prism Domain
node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
elem = [1,2,3,7; 1,6,2,7; 1,5,6,7];
elem = label3(node,elem);
figure;
showmesh3(node,elem);
view([-53,8]);
findnode3(node,[1 2 3 5 6 7]);
findelem3(node,elem);
bdFlag = setboundary3(node,elem,'Dirichlet','(z==1) | (z==-1)');
display(elem)
display(bdFlag)

%%
% The top and bottom of the prism is set as Dirichlet boundary condition
% and other faces are zero flux boundary condition.

%% Remark
% It would save storage if we record boundary edges or faces only. The
% current data structure is convenient for the local refinement and
% coarsening since the boundary can be easily update along with the change
% of elements. The matrix |bdFlag| is sparse but we use a dense matrix to
% store it. We do not save |bdFlag| as a sparse matrix since updating
% sparse matrix is time consuming. We set up the type of |bdFlag| or
% |bdFlag| to |uint8| to minimize the waste of spaces.

repo: Aero-matlab


text:
%Detection
%=========
%
%Calculation of the likelihood ratio of ARMA models.
%
%  likelihoodR         - Likelihood ratio
%  likelihoodR_mod     - modified likelihood ratio: includes
%                        a correction for the number of
%                        estimated parameters. This can be used
%                        for detection with models that have been
%                        selected using statistical order selection
%                        criteria.


text:
function [grad] = scaleGrad(grad, scale)
% scaleGrad just scales the gradint with the scalar
    if strcmp(grad.type,'trap')
        grad.amplitude=grad.amplitude*scale;
        grad.area=grad.area*scale;
        grad.flatArea=grad.flatArea*scale;
    else
        grad.waveform=grad.waveform*scale;
        grad.first=grad.first*scale;
        grad.last=grad.last*scale;
    end
    if isfield(grad,'id')
        grad = rmfield(grad,'id');
    end
end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function value=moxunit_fieldtrip_util_parse_dependency(line)
% parse WALLTIME property in FieldTrip test
%
% value=moxunit_fieldtrip_util_parse_walltime(line)
%
% Input:
%   line                string which represents one or more alphanumeric
%                       strings
%
% Output:
%   value               cell string with the strings in line separated by
%                       white space
%
% #   For MOxUnit_fieldtrip's copyright information and license terms,   #
% #   see the COPYING file distributed with MOxUnit_fieldtrip.           #

    if ~ischar(line)
        error('input must be numeric');
    end

    value=regexp(line,'(\w+)','match');


text:
function cdf = f_cdf ( x, m, n )

%*****************************************************************************80
%
%% F_CDF evaluates the F central CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Formula 26.5.28
%    Abramowitz and Stegun,
%    Handbook of Mathematical Functions.
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%
%    Input, integer M, N, the parameters of the PDF.
%    1 <= M,
%    1 <= N.
%
%    Output, real CDF, the value of the CDF.
%
  if ( x <= 0.0 )

    cdf = 0.0;

  else

    arg1 = 0.5 * n;
    arg2 = 0.5 * m;
    arg3 = n / ( n + m * x );

    cdf = beta_inc ( arg1, arg2, arg3 );

  end

  return
end