text:
function [lf] = eeg_infinite_monopole(monpos, elc, vol)

% EEG_INFINITE_MONOPOLE calculate the infinite medium potential for a monopole
%
% Use as
%   [lf] = eeg_infinite_monopole(monpos, elc, vol)
%
% Implemented from Malmivuo J, Plonsey R, Bioelectromagnetism (1993)
% http://www.bem.fi/book/08/08.htm
%
% See also EEG_INFINITE_DIPOLE, EEG_HALFSPACE_DIPOLE, EEG_HALFSPACE_MONOPOLE

% Copyright (C) 2011, Cristiano Micheli
% Copyright (C) 2019, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if ~isstruct(vol)
  % it only represents the conductivity, make a structure out of it
  vol = struct('cond', vol);
end

siz = size(monpos);
if any(siz==1)
  % positions are specified as a single vector
  Npoles = prod(siz)/3;
  monpos = monpos(:)'; % ensure that it is a row vector
elseif siz(2)==3
  % positions are specified as a Nx3 matrix -> reformat to a single vector
  Npoles = siz(1);
  monpos = monpos';
  monpos = monpos(:)'; % ensure that it is a row vector
else
  ft_error('incorrect specification of monopole locations');
end

cond     = vol.cond;
Nelc     = size(elc,1);
lf       = zeros(Nelc,Npoles);

mu0   = 4*pi*1e-7;         % Permeability of free space
c     = 2.99792458 * 1e8;  % Speed of light
e0    = 1 / (mu0*c^2);     % Permittivity of Free Space

for i=1:Npoles
  % this is the position of monopole "i"
  monopole = monpos((1:3) + 3*(i-1));

  % distances from electrodes to monopole
  r = elc - ones(Nelc,1) * monopole;
  r = sqrt(sum(r.^2,2));

  lf(:,i) = 1 ./ (4*pi*cond*r);
end


text:
function [metrics, modelEssential] = evalGeneEssentialityPred(model, expData, metabTasks, modelEssential)
% Compare GEM-based gene essentiality predictions with experimental data
%
% Inputs:
%
%   model       Genome-scale metabolic model structure.
%
%               NOTE: The model should include boundary metabolites! These
%               can be added using the addBoundaryMets function.
%
%   expData     A two-column cell array, where the first column contains
%               gene names or IDs (of the same type as those used in the
%               model), and the second column indicates whether the gene
%               was found as essential or non-essential.
%
%               The second column can be text (e.g., 'essential', 'non') or
%               numeric (0 = non-essential, 1 = essential). For example:
%
%               {'gene1'   'essential'  }       {'gene1'   [1]}
%               {'gene2'   'non'        }       {'gene2'   [0]}
%               {'gene3'   'conditional'}       {'gene3'   [1]}
%               {...       ...          }       {...       ...}
%
%               NOTE: If essentiality is provided as strings, only genes
%               starting with "non" will be treated as non-essential. This
%               means that genes annotated as "essential" or "conditional"
%               will both be treated as essential.
%
%               NOTE: ALL genes that were tested for essentiality should be
%               included in expData. If expData contains only essential
%               genes, then it will be assumed that all genes in the genome
%               (and thus all genes in the model) were tested.
%
%   metabTasks  Either the filename or structure containing the metabolic
%               tasks that will be tested when using the model to predict
%               gene essentiality. The task structure can be loaded using
%               the RAVEN "parseTaskList" function.
%
%               Genes will be predicted as essential if their deletion
%               impairs ANY of the metabolic tasks in metabTasks. If the
%               model can perform ALL metabolic tasks upon the deletion of
%               a gene, that gene will be predicted as non-essential.
%
%   modelEssential  (Optional) A list of model-predicted essential genes.
%                   If provided, the model gene essentiality prediction
%                   will be skipped and the provided list used instead.
%
% Outputs:
%
%   metrics     A result structure with the following fields:
%                   sensitivity
%                   specificity
%                   accuracy
%                   F1 statistic
%                   Matthew's Correlation Coefficient (MCC)
%                   p-value associated with a hypergeometric test of the
%                           true and false positives and negatives
%                           (2x2 contingency table).
%
%   modelEssential   List of model-predicted essential genes.
%

if nargin < 4
    modelEssential = [];
end

%% Pre-process expData

% extract information from expData and convert essentiality to 0, 1
ex.genes = expData(:,1);
if length(ex.genes) > length(unique(ex.genes))
    error('Gene names or IDs in expData are not unique. Duplicated entries must be removed.');
end
if isnumeric(expData{1,2})
    ex.essentiality = double(cell2mat(expData(:,2)));
else
    ex.essentiality = double(~startsWith(lower(expData(:,2)), 'non'));
end

% check that gene IDs/names in the model are the same type as in expData
if sum(ismember(ex.genes, model.genes)) < 3
    error('The gene name or ID type in expData seem to differ from those used in the model.');
end

% define essential and non-essential gene lists
ex.essential = ex.genes(ex.essentiality == 1);
if all(ex.essentiality == 1)
    fprintf('NOTE: All genes in expData are marked as essential; it will therefore be assumed that all genes in the genome were tested.\n\n');
    ex.nonessential = setdiff(model.genes, ex.essential);
    ex.genes = [ex.genes; ex.nonessential];
else
    ex.nonessential = ex.genes(ex.essentiality == 0);
end

%% Run gene essentiality predictions with the model

if isempty(modelEssential)

    if ischar(metabTasks)
        taskStruct = parseTaskList(metabTasks);
    else
        taskStruct = metabTasks;
    end

    % first confirm that the model can perform all the tasks
    taskReport = checkTasks(model,[],false,false,false,taskStruct);
    if ~all(taskReport.ok)
        fprintf('\nThe provided model could not perform the following tasks:\n');
        fprintf('\t> %s\n', taskReport.description{~taskReport.ok});
        fprintf('\n');
        error('Model failed task(s) before gene deletion!');
    end

    % get essential genes for each task in taskStruct
    [~,essentialGeneMat] = checkTasksGenes(model,[],false,false,true,taskStruct);

    % essential genes are counted as those that are essential for ANY task
    essentialGeneVect = any(essentialGeneMat > 0, 2);
    pred.essential = model.genes(essentialGeneVect);
    pred.nonessential = setdiff(model.genes, pred.essential);

else
    pred.essential = modelEssential;
    pred.nonessential = setdiff(model.genes, modelEssential);
end

%% Evaluate prediction performance

TP = sum(ismember(pred.essential, ex.essential));        % true positives
TN = sum(ismember(pred.nonessential, ex.nonessential));  % true negatives
FP = sum(ismember(pred.essential, ex.nonessential));     % false positives
FN = sum(ismember(pred.nonessential, ex.essential));     % false negatives

% calculate some metrics
sensitivity = TP./(TP + FN);
specificity = TN./(TN + FP);
accuracy = (TP + TN)./(TP + TN + FP + FN);
F1 = 2*TP./(2*TP + FP + FN);
MCC = ((TP.*TN) - (FP.*FN))./sqrt((TP+FP).*(TP+FN).*(TN+FP).*(TN+FN));
[~, p_hyper] = fishertest([TP, FP; FN, TN], 'tail', 'right');

% combine metrics into output structure
metrics.sensitivity = sensitivity;
metrics.specificity = specificity;
metrics.accuracy = accuracy;
metrics.F1 = F1;
metrics.MCC = MCC;
metrics.p_hyper = p_hyper;

% assign output
modelEssential = pred.essential;

end

author: SysBioChalmers

repo: Human-GEM


text:
classdef C10MOP6 < EvoXBenchProblem
% <multi> <real> <large/none> <expensive/none>

%------------------------------- Reference --------------------------------
% Z. Lu, R. Cheng, Y. Jin, K. C. Tan, and K. Deb, Neural architecture
% search as multiobjective optimization benchmarks: Problem formulation and
% performance assessment, IEEE Transactions on Evolutionary Computation,
% 2023.
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            config.name = 'nb201';
            config.args.fidelity = 200;
            config.args.objs = ['err&params&flops&eyeriss_latency&' ...
                'eyeriss_energy&eyeriss_arithmetic_intensity'];
            config.args.dataset = 'cifar10';
            config.args.normalized_objectives = true;
            obj.Setting@EvoXBenchProblem(config);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function out = PlotPath(state, scalling, mapName)
% uses B-spline interpolation
    path = state.path';
    path(:,3) = 1;
    path(:,4) = 1;

    if nargin <3
        scalling = 1;
        imag = state.map;
    else
        tmp = LoadMap(  strcat(mapName, '.png'), 1);
        imag = tmp.map;
    end

    if scalling ~= 1
        path(:,1) = path(:,1)-2.5;
        path(:,2) = path(:,2)-2.5;
        pathInterpolated = BSpline(path*scalling);
    else
        pathInterpolated = path;
    end
    %%5
    for it = pathInterpolated'
        a = round(it);
        imag(a(1), a(2)) = 0.6;
    end
    figure(100)
    out.handle = imshow(imag);
    out.path = pathInterpolated;

end

author: LazyFalcon


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function bb = bb_shift_relative(bb,shift)
% Change

if isempty(bb)
    return;
end
bb(1,:) = bb(1,:) + bb_width(bb)*shift(1);
bb(2,:) = bb(2,:) + bb_height(bb)* shift(2);
bb(3,:) = bb(3,:) + bb_width(bb)*shift(1);
bb(4,:) = bb(4,:) + bb_height(bb)*shift(2);

author: yuxng


text:
function cvx_optval = norms_largest( varargin )

%NORMS_LARGEST Computation of multiple norm_largest() norms.
%   NORMS_LARGEST( X, K, DIM ) provides a means to compute the largest-k
%   norms of multiple vectors packed into a matrix or N-D vector. This is
%   useful for performing max-of-norms or sum-of-norms calculations.
%
%   If DIM is omitted, the norms are computed along the first non-singleton
%   dimension.
%
%   See NORM_LARGEST.
%
%   Disciplined convex programming information:
%       NORMS_LARGEST is convex and non-monotonic, so its input must be affine.

[ sx, x, k, dim ] = cvx_get_dimension( varargin, 3 );

if ~isnumeric( k ) || ~isreal( k ) || length( k ) ~= 1,
    cvx_throw( 'Second argument must be a real scalar.' );
end

cvx_optval = sum_largest( abs( x ), k, dim );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% Redundancy reduction principle
% This demo performs BMC on two nested models, one of which is fully
% redundant. The data is then varied on a pre-defined grid, to assess the
% impact of the complexity penalty term as a function of the data.

clear variables
close all

dim1.n_phi = 1;
dim1.n_theta = 0;
dim1.n=0;
dim2.n_phi = 2;
dim2.n_theta = 0;
dim2.n=0;
opt1.inG.X = 1;
opt1.DisplayWin = 0;
opt1.verbose = 0;
opt2.inG.X = [1,1];
opt2.DisplayWin = 0;
opt2.verbose = 0;

gy = 2.^[-8:8];
ny = length(gy);
for i=1:ny
    i
    y = gy(i);
    [p1,o1] = VBA_NLStateSpaceModel(y,[],[],@g_GLM,dim1,opt1);
    set(gcf,'tag','1')
    [p2,o2] = VBA_NLStateSpaceModel(y,[],[],@g_GLM,dim2,opt2);
    dF(i) = o1.F(end) - o2.F(end);
    m1(i) = p1.muPhi;
    s1(i) = p1.SigmaPhi;
    m2(i) = ones(2,1)'*p2.muPhi;
    s2(i) = ones(2,1)'*p2.SigmaPhi*ones(2,1);
    c = VBA_cov2corr(p2.SigmaPhi);
    c2(i) = c(1,2);
end

hf = figure('color',[1 1 1]);
ha = subplot(2,2,1,'parent',hf);
plot(ha,gy,dF,'marker','.')
set(ha,'xscale','log')
ylabel(ha,'log p(y|m1) - log p(y|m2)')
xlabel(ha,'y')
title(ha,'Bayesian model comparison')
ha = subplot(2,2,2,'parent',hf);
plot(ha,gy,c2,'marker','.')
set(ha,'xscale','log')
ylabel(ha,'Corr[x1,x2|y,m2]')
xlabel(ha,'y')
title(ha,'parameter identifiability under m2')
ha = subplot(2,2,3,'parent',hf);
plotUncertainTimeSeries(m1,s1,gy,ha)
set(ha,'xscale','log')
ylabel(ha,'E[x1|y,m1]')
xlabel(ha,'y')
title(ha,'estimated effect under m1')
ha = subplot(2,2,4,'parent',hf);
plotUncertainTimeSeries(m2,s2,gy,ha)
set(ha,'xscale','log')
ylabel(ha,'E[x1+x2|y,m2]')
xlabel(ha,'y')
title(ha,'estimated effect under m2')

author: MBB-team

repo: VBA-toolbox


text:
function [fv, opt]= proc_subtractMean(fv, varargin)
% PROC_SUBTRACTMEAN - Subtract mean or median from data, given in fv.
%
%Synopsis:
% [fv, opt]= proc_normalize(fv, opt)
%
%Arguments:
%     fv    - struct of feature vectors
%     opt  struct and/or property/value list of properties
%      .Policy - one of 'mean' (default), 'median', 'nanmean', 'nanmedian'
%      .Dim    - dimension along which the mean should be computed.
%                1 computes along the first dimension,
%                2 computes along the second dimension (default).
%      .Bias   - vector which is subtracted from fv. If this option is
%                given, the 'Policy' option is ignored.
%                Typically the bias vector is computed in a first call to
%                proc_subtractMean. Thus, this field can be used to
%                apply the shift calculated from one data set (e.g.
%                training data) to another data set (e.g. test data)
%
%      fv   - struct of shifted feature vectors
%      opt  - a copy of the input options, with a new field .Bias that
%             contains the mean/median vector that has been
%             subtracted from the data.
%
% See also: nanmean

% bb 09/03, ida.first.fhg.de
% Anton Schwaighofer, Feb 2005
% Sven Daehne, Feb 2016, ported to git toolbox

fv = misc_history(fv);

props= { 'Policy'  'mean'    '!CHAR(mean median nanmean nanmedian)'
         'Dim'      2       'INT'
         };

if nargin==0,
  fv = props; return
end

misc_checkType(fv, 'STRUCT(x)');

% this allows the function to be called like
% "proc_subtractMean(fv, 'mean')", i.e. without explictly stating 'Policy'
if length(varargin)==1 && ischar(varargin{1}),
  opt= struct('Policy', varargin{1});
else
  opt= opt_proplistToStruct(varargin{:});
end

opt = opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

sz= size(fv.x);

if isfield(opt, 'Bias'),
  bsz= sz;
  bsz(opt.Dim)= 1;
  if ~isequal(size(opt.Bias), bsz),
    error('size of opt.Bias does not match fv');
  end
else
  switch(opt.Policy),
   case 'mean',
    opt.Bias= mean(fv.x, opt.Dim);
   case 'median',
    opt.Bias= median(fv.x, opt.Dim);
   case 'nanmean',
    opt.Bias= nanmean(fv.x, opt.Dim);
   case 'nanmedian',
    opt.Bias= nanmedian(fv.x, opt.Dim);
  end
end

rep_sz= ones(1, max(length(sz), opt.Dim));
rep_sz(opt.Dim)= sz(opt.Dim);

fv.x= fv.x - repmat(opt.Bias, rep_sz);

author: bbci


text:
function y=PSNR_RGB(X,Y)

% Y= PSNR_RGB(X,Y)
% Computes the Peak Signal to Noise Ratio for two RGB images
% Class input : double [0,1] ,
% july ,25, 2012
% KHMOU Youssef

if size(X)~=size(Y)
    error('The images must have the same size');
end

%if ~isa(X,'double')
%   X=double(X)./255.00;
%end
%if  ~isa(Y,'double')
%    Y=double(Y)./255.00;
%end

% begin
d1=max(X(:));
d2=max(Y(:));
d=max(d1,d2);
sigma=mean2((X-Y).^2);

y=10*log((d.^2)./sigma);

author: Sable

repo: mcbench-benchmarks


text:
function d=ltfatgetdefaults(fname)
%LTFATGETDEFAULTS  Get default parameters of function
%
%   `ltfatgetdefaults(fname)` returns the default parameters
%   of the function `fname` as a cell array.
%
%   `ltfatgetdefaults('all')` returns all the set defaults.
%
%   See also: ltfatsetdefaults, ltfatstart

if nargin<1
    error('%s: Too few input arguments',upper(mfilename));
end;

if strcmpi(fname,'all')
  d=ltfatarghelper('all');
else
  d=ltfatarghelper('get',fname);
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function a = p04_a ( m, n )

%*****************************************************************************80
%
%% P04_A returns the matrix A for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of equations.
%
%    Input, integer N, the number of variables.
%
%    Output, real A(M,N), the matrix.
%
  for j = 1 : n
    for i = 1 : m
      a(i,j) = j ^ ( i - 1 );
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%
% Call the Fit() method of the model to fit
%
function AMICO_Fit()

	global CONFIG niiSIGNAL niiMASK bMATRIX

	% Fit right model to each voxel
    % =============================
	if ~isempty(CONFIG.model)

        % fit the model to the data
        fprintf( '\n-> Fitting "%s" model to %d voxels:\n', CONFIG.model.name, nnz(niiMASK.img) );
        TIME = tic;
        progress = ProgressBar( nnz(niiMASK.img) );
        for iz = 1:niiSIGNAL.hdr.dime.dim(4)
        for iy = 1:niiSIGNAL.hdr.dime.dim(3)
        for ix = 1:niiSIGNAL.hdr.dime.dim(2)
            if niiMASK.img(ix,iy,iz)==0, continue, end
            progress.update();

            try
                % Read the signal
                b0 = mean( squeeze( niiSIGNAL.img(ix,iy,iz,CONFIG.scheme.b0_idx) ) );
                if ( b0 < 1e-3 ), continue, end
                y = double( squeeze( niiSIGNAL.img(ix,iy,iz,:) ) ./ ( b0 + eps ) );
                y( y < 0 ) = 0; % [NOTE] this should not happen!

                % Find the MAIN DIFFUSION DIRECTIONS
                if any(strcmp(properties(CONFIG.model), 'max_dirs'))==false || CONFIG.model.max_dirs>0
                    % using DTI
                    [ ~, ~, V ] = AMICO_FitTensor( y, bMATRIX );
                    vox_DIRs = V(:,1);
                    if ( vox_DIRs(2)<0 ), vox_DIRs = -vox_DIRs; end
                    [ i1, i2 ] = AMICO_Dir2idx( vox_DIRs );
                    DIRs(ix,iy,iz,:) = vox_DIRs;
                else
                    % not needed by the model
                    DIRs(ix,iy,iz,:) = 0;
                    i1 = 1;
                    i2 = 1;
                end

                % Dispatch to the right handler for each model
                vox_MAPs = CONFIG.model.Fit( y, i1, i2 );

                % Store results
                MAPs(ix,iy,iz,:) = vox_MAPs;
            catch exception
                % set output maps to NaN in case of problems
                MAPs(ix,iy,iz,:) = NaN;
            end
        end
        end
        end
        progress.close();

        TIME = toc(TIME);
        fprintf( '   [ %.0fh %.0fm %.0fs ]\n', floor(TIME/3600), floor(mod(TIME/60,60)), mod(TIME,60) )
		CONFIG.OPTIMIZATION.fit_time = TIME;
	else
		error( '[AMICO_Fit] Model not set' )
    end

	% Save CONFIGURATION and OUTPUT to file
    % =====================================
    fprintf( '\n-> Saving output to "AMICO/*":\n' );

    fprintf( '\t- CONFIG.mat' );
    save( fullfile(CONFIG.OUTPUT_path,'CONFIG.mat'), '-v6', 'CONFIG' )
    fprintf( ' [OK]\n' );

    fprintf( '\t- FIT_dir.nii' );
    niiMAP = niiMASK;
    niiMAP.hdr.dime.dim(1) = 4;
    niiMAP.hdr.dime.datatype = 16;
    niiMAP.hdr.dime.bitpix = 32;
    niiMAP.hdr.dime.glmin = -1;
    niiMAP.hdr.dime.glmax = 1;
    niiMAP.hdr.dime.calmin = niiMAP.hdr.dime.glmin;
    niiMAP.hdr.dime.calmax = niiMAP.hdr.dime.glmax;
    niiMAP.hdr.dime.scl_slope = 1;
    niiMAP.hdr.dime.scl_inter = 0;

    niiMAP.img = DIRs;
    niiMAP.hdr.dime.dim(5) = size(DIRs,4);
    save_untouch_nii( niiMAP, fullfile(CONFIG.OUTPUT_path,'FIT_dir.nii') );
    fprintf( ' [OK]\n' );

    niiMAP.hdr.dime.dim(5) = 1;
    for i = 1:numel(CONFIG.model.OUTPUT_names)
        fprintf( '\t- AMICO/FIT_%s.nii', CONFIG.model.OUTPUT_names{i} );

        niiMAP.img = MAPs(:,:,:,i);
        niiMAP.hdr.hist.descrip = CONFIG.model.OUTPUT_descriptions{i};
        niiMAP.hdr.dime.glmin = min(niiMAP.img(:));
        niiMAP.hdr.dime.glmax = max(niiMAP.img(:));
        niiMAP.hdr.dime.calmin = niiMAP.hdr.dime.glmin;
        niiMAP.hdr.dime.calmax = niiMAP.hdr.dime.glmax;
        save_untouch_nii( niiMAP, fullfile(CONFIG.OUTPUT_path,['FIT_' CONFIG.model.OUTPUT_names{i} '.nii']) );

        fprintf( ' [OK]\n' );
    end

    fprintf( '   [ DONE ]\n' )
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% pop_writeeeg - write EEGLAB dataset to disk in EDF/GDF or BDF format
%
% pop_writeeeg( EEG ) % pops up a window
% pop_writeeeg( EEG, filename, 'key', 'val' )
%
% Inputs:
%  EEG        - EEGLAB dataset
%  filename   - [string] filename
%
% Optional keys (same as writeeeg):
%  'TYPE'         - ['GDF'|'EDF'|'BDF'|'CFWB'|'CNT'] file format for writing
%                   default is 'EDF'.
%  See writeeeg for more information
%
% Author: Arnaud Delorme, SCCN, UCSD/CERCO, 2009
%         Based on BIOSIG, sopen and swrite

% Copyright (C) 22 March 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% 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 [command] = pop_writeeeg(EEG, filename, varargin);
command = '';

if nargin < 2
	% ask user
	[filename, filepath] = uiputfile('*.*', 'Enter a file name -- pop_writeeeg()');
	if filename == 0 return; end;
	filename = fullfile(filepath,filename);

    % file format
    % -----------
    fileformats = { 'EDF' 'GDF' 'BDF' };
    uilist = { { 'style' 'text' 'String' 'File format' } ...
               { 'style' 'listbox' 'string' strvcat(fileformats) 'value' 1 } };
    geom = [1 1];
    result = inputgui( 'geometry', geom, 'uilist', uilist, 'helpcom', 'pophelp(''pop_writeeeg'')', ...
                     'title', 'Write data using BIOSIG -- pop_writeeeg()', 'geomvert', [1 2.5]);
    if length(result) == 0 return; end;

    options = { 'TYPE' fileformats{result{1}} };
else
    options = varargin;
end;

warning('off', 'MATLAB:intConvertNonIntVal');
if ~isempty(EEG.chanlocs)
    tmpchanlocs = EEG.chanlocs;
    writeeeg(filename, EEG.data, EEG.srate, 'label', { tmpchanlocs.labels }, 'EVENT', EEG.event, options{:});
else
    writeeeg(filename, EEG.data, EEG.srate, 'EVENT', EEG.event, options{:});
end;
warning('on', 'MATLAB:intConvertNonIntVal');

command = sprintf('pop_writeeeg(EEG, ''%s'', %s);', filename, vararg2str(options));

return;

author: PatternRecognition

repo: OpenBMI


text:
function x = combin_eigen_right ( alpha, beta, n )

%*****************************************************************************80
%
%% COMBIN_EIGEN_RIGHT returns the right eigenvectors of the COMBIN matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, BETA, scalars that define A.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real X(N,N), the right eigenvectors.
%
  x(1:n,1:n) = 0.0;

  for j = 1 : n - 1
    x(  1,j) = +1.0;
    x(j+1,j) = -1.0;
  end

  j = n;
  x(1:n,j) = 1.0;

  return
end


text:
function [params, names] = linearExtractParam(model,dim);

% LINEAREXTRACTPARAM Extract weights from a linear model.
%
% MODIFICATIONS : Carl Henrik Ek, 2009
%
% MLTOOLS

if(nargin<2)
  params = [model.W(:)' model.b];
else
  params = model.W(:,dim);
  params = [params(:)' model.b(dim)];
end

if nargout > 1
  counter = 0;
  for j = 1:size(model.W, 2)
    for i = 1:size(model.W, 1)
      counter = counter + 1;
      names{counter} = ['Weight ' num2str(i) '-' num2str(j)];
    end
  end
    for j = 1:size(model.b, 2)
    counter = counter + 1;
    names{counter} = ['Bias ' num2str(j)];
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [I] = imgread(img, c)
%
% Copyright Aditya Khosla http://mit.edu/khosla
%
% Please cite this paper if you use this code in your publication:
%   A. Khosla, J. Xiao, A. Torralba, A. Oliva
%   Memorability of Image Regions
%   Advances in Neural Information Processing Systems (NIPS) 2012
%

if(ischar(img))
    I = imread(img);
else
    I=img;
end

scale = min(min(c.maxsize/size(I, 1), c.maxsize/size(I, 2)), 1);
if(scale~=1)
	I = imresize(I, scale);
end

if(size(I,3)==1)
	I=cat(3, I, I, I);
end

author: adikhosla

repo: feature-extraction


text:
function [cache_exists,cache_name] = find_cache()
  % FIND_CACHE Called at the beginning of a very expensive caller function,
  % this function writes the input parameters (whatever's currently in scope)
  % to a temporary file and gets that files md5 checksum as a "cachename". If a
  % file with path built from that cache name and caller name exists and is
  % younger than the caller file name or any of its dependencies, then this
  % returns true and the cache name.
  %
  % [cache_exists,cache_name] = find_cache()
  %
  % Outputs:
  %   cache_exists  true if the cache_name file exists
  %   cache_name  path to cache file for this caller function on these inputs.
  %
  % See also: create_cache

  %if ~ismac
  if ~isunix
    cache_exists = false;
    cache_name = '';
  end

  S = dbstack(1);
  caller_name = S.name;
  % get a list of current variables in this scope, this is the input "state"
  variables = evalin('caller','who');
  % get a temporary file's name
  tmpf = [tempname('.') '.mat'];

  variable_names = sprintf('^%s$|',variables{:});
  save_cmd = sprintf('save(''%s'',''-regexp'',''%s'');',tmpf,variable_names);
  % save the "state" of the caller to file, so we can get a md5 checksum
  evalin('caller',save_cmd);

  % get md5 checksum on input "state", we append .cache.mat to the check sum
  % because we'll use the checksum as the cache file name
  if ismac
    md5_cmd = 'md5 -r';
  elseif isunix
    md5_cmd = 'md5sum';
  else
    assert(false,'What OS is this?')
  end
  [s,md5_res] = system(sprintf('tail -c +117 %s | %s',tmpf,md5_cmd));
  if s~=0
    cache_exists = false;
    cache_name = '';
    return;
  end
  md5_res = sscanf(md5_res,'%s',1);
  cache_name = sprintf('.%s.cache.%s.mat',caller_name,md5_res);

  % clean up
  delete(tmpf);
  clear s tmpf variables;
  cache_exists = exist(cache_name,'file')~=0;
  % If cache exists, check that it's older than the caller file.
  if cache_exists
    cache_file = dir(cache_name);
    cache_date = datenum(cache_file.date);
    caller_dependencies = depends(caller_name);
    for si = 1:numel(caller_dependencies)
      caller_dep_name = caller_dependencies{si};
      caller_dep_file = dir(which(caller_dep_name));
      caller_dep_date = datenum(caller_dep_file.date);
      if cache_date < caller_dep_date
        cache_exists = false;
      end
    end
  end
end

author: alecjacobson

repo: gptoolbox


text:
function [V,compareFluxes] = exoMetLOOCV(model,measuredFluxes,param)
% model: COBRA model
%
% measuredFluxes:   table with the fluxes obtained from exometabolomics experiments with the following variables
%           * .rxns:    k x 1 cell array of reaction identifiers
%           * .mean:    k x 1 double of measured mean flux
%           * .SD:      k x 1 double standard deviation of the measured flux
%           * .labels:  k x 1 cell array of labels to display on y axis
%           * .Properties.Description: string describing the data, used in plot legend
%
%  param.alpha:              alpha in (alpha/2)(v-h)'*H*(v-h), default = 10000;
%  param.metabolomicWeights: String indicating the type of weights to be applied to penalise the difference
%                            between of predicted and experimentally measured fluxes by, where
%                            'SD'   weights = 1/(1+exoMet.SD^2)
%                            'mean' weights = 1/(1+exoMet.mean^2)
%                            'RSD'  weights = 1/((exoMet.SD./exoMet.mean)^2)
%  param.relaxBounds: True to relax bounds on reaction whose fluxes are to be fitted to exometabolomic data
%
%
% OUTPUTS:
%  V:     table with the fluxes predicted by leave one out cross validation
%           * .rxns:    n x 1 cell array of reaction identifiers, one for each in input model
%           * .v:       n x 1 double of predicted mean flux
%           * .lb:      n x 1 double lower bound on reaction flux
%           * .ub:      n x 1 double upper bound on reaction flux
%
% LIBkey:
%
% EXAMPLE:
%
% NOTE:
%
% Author(s):
if ~exist('param','var')
    param = struct;
end
if ~isfield(param,'alpha')
    param.alpha=10000;
end
if ~isfield(param,'relaxBounds')
    param.relaxBounds=1;
end
if ~isfield(param,'metabolomicWeights')
    param.metabolomicWeights = 'mean';
end
if ~isfield(param,'relaxBounds')
    param.relaxBounds=1;
end
if ~isfield(param,'approach')
    param.approach='QEFBA';
end

[~,nRxns] = size(model.S);

if param.relaxBounds
    bool=ismember(model.rxns,measuredFluxes.rxns);
    model.lb(bool) = model.lb_preconstrainRxns(bool);
    model.ub(bool) = model.ub_preconstrainRxns(bool);
end

modelOrig = model;

if 1
    [measuredFluxes,LIBkey,LOCAkey] = mapAontoB(measuredFluxes.rxns,model.rxns,measuredFluxes);

    %preallocate and add data from measurments
    V = zeros(nRxns,nnz(LIBkey)+4);
    varNames = cell(nnz(LIBkey)+4,1);
    V(:,1) = measuredFluxes.mean;
    varNames{1} = 'mean';
    V(:,2) = measuredFluxes.SD;
    varNames{2} = 'SD';

    %reduce to the measured fluxes part of the model
    measuredFluxes = measuredFluxes(LIBkey,:);
else
    nMeasuredRxns = size(measuredFluxes,1);

    %preallocate and add data from measurments
    V = zeros(nRxns,nMeasuredRxns+4);
    varNames = cell(nMeasuredRxns+4,1);
    [V(:,1),LIBkey] = mapAontoB(measuredFluxes.rxns,model.rxns,measuredFluxes.mean);
    varNames{1} = 'mean';
    V(:,2) = mapAontoB(measuredFluxes.rxns,model.rxns,measuredFluxes.SD);
    varNames{2} = 'SD';
end

nMeasuredRxns = size(measuredFluxes,1);
bool = true(nMeasuredRxns,1);

%loop through the measured fluxes, leaving one of them out each time
for i=1:nMeasuredRxns
    bool(i)=0;
    switch param.approach
        case 'QEFBA'
            if i==1
                %compare without leaving any measured flux out
                model = addExoMetToEFBA(modelOrig,measuredFluxes,param);
                efbaParam.method = 'fluxes';
                efbaParam.printLevel = 0;
                model.osenseStr = 'min';
                model.cf = 0;
                model.cr = 0;
                model.g = 2;
                model.u0 = 0;
                model.f = 1;
                [QEFBAsolution, ~] = entropicFluxBalanceAnalysis(model,efbaParam);
                V(:,3) = QEFBAsolution.v;
                varNames{3}='QEFBA';
            end

            model = changeRxnBounds(modelOrig, measuredFluxes.rxns(~bool), min(modelOrig.lb), 'l');
            model = changeRxnBounds(model, measuredFluxes.rxns(~bool), max(modelOrig.ub), 'u');
            model = addExoMetToEFBA(model,measuredFluxes(bool,:),param);
            efbaParam.method = 'fluxes';
            efbaParam.printLevel = 0;
            model.osenseStr = 'min';
            model.cf = 0;
            model.cr = 0;
            model.g = 2;
            model.u0 = 0;
            model.f = 1;
            [solution, modelOut] = entropicFluxBalanceAnalysis(model,efbaParam);
            V(:,i+4) = solution.v;
            varNames{i+4}=measuredFluxes.rxns{~bool};

            %generate a flux vector amalgamated from all of the leave one out predictions
            boolRxn = ismember(model.rxns,measuredFluxes.rxns(~bool));
            V(boolRxn,4) = solution.v(boolRxn);
            if i==1
                varNames{4}='LOOCV';
            end
    end
    bool(i)=1;
end

V = array2table(V,'VariableNames',varNames);
V = addvars(V,model.rxns,modelOrig.lb,modelOrig.ub,'NewVariableNames',{'rxns','lb','ub'},'Before','mean');
% V = V(LIBkey,:);

compareFluxes = mapAontoB('rxns','rxns',V,measuredFluxes);

author: opencobra

repo: cobratoolbox


text:
function [ xtab, weight ] = ncc_rule ( norder )

%*****************************************************************************80
%
%% NCC_RULE computes the coefficients of a Newton-Cotes closed quadrature rule.
%
%  Discussion:
%
%    For the interval [-1,1], the Newton-Cotes quadrature rule estimates
%
%      Integral ( -1 <= X <= 1 ) F(X) dX
%
%    using NORDER equally spaced abscissas XTAB(I) and a weight vector
%    WEIGHT(I):
%
%      Sum ( 1 <= I <= N ) WEIGHT(I) * F ( XTAB(I) ).
%
%    For the CLOSED rule, the abscissas include A and B.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NORDER, the order of the rule.
%
%    Output, real XTAB(NORDER), the abscissas of the rule.
%
%    Output, real WEIGHT(NORDER), the weights of the rule.
%

%
%  Compute a closed quadrature rule.
%
  a = -1.0;
  b =  1.0;

  for i = 1 : norder
    xtab(i) = ( ( norder - i     ) * a   ...
              + (          i - 1 ) * b ) ...
              / ( norder - 1     );
  end

  weight = nc_rule ( norder, a, b, xtab );

  return
end


text:
function C = boundaryConditions(X0, XF)
% C = boundaryConditions(X0, XF)
%
% This function computes the boundary conditions for periodic solutions.
%
% INPUTS:
%   X0 = Initial State
%   XF = Final State
%
% OUTPUTS:
%   C = defect constraint
%

C = X0 - XF;

end

author: MatthewPeterKelly

repo: dscTutorials


text:
%-Abstract
%
%   ZZMICE_DP converts an numeric input to double precision format.
%
%-Disclaimer
%
%   THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
%   CALIFORNIA  INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
%   GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
%   ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
%   PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED
%   "AS-IS" TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING
%   ANY WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR
%   A PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
%   SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
%   SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
%
%   IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY,
%   OR NASA BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING,
%   BUT NOT LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF
%   ANY KIND, INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY
%   AND LOST PROFITS, REGARDLESS OF WHETHER CALTECH, JPL, OR
%   NASA BE ADVISED, HAVE REASON TO KNOW, OR, IN FACT, SHALL
%   KNOW OF THE POSSIBILITY.
%
%   RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE
%   OF THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO
%   INDEMNIFY CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING
%   FROM THE ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
%
%-I/O
%
%   Given:
%
%      x   an input numeric to convert to double precision
%
%   the call:
%
%      y = zzmice_dp(x)
%
%   returns:
%
%      y   the double precision representation of 'x'
%
%-Examples
%
%   None.
%
%-Particulars
%
%   This routine exists to support the NAIF MATLAB-CSPICE interface.
%
%-Required Reading
%
%   MICE.REQ
%
%-Version
%
%   -Mice Version 1.1.0, 27-JUL-2009, EDW (JPL)
%
%      Added value check on 'nargin'. Incorrect input argument type/form
%      error tag changed from "MICE(BADVAL)" to "MICE(BADARG)."
%
%   -Mice Version 1.0.1, 30-DEC-2008, EDW (JPL)
%
%      Function ensures all input values as finite.
%
%      Corrected misspellings.
%
%   -Mice Version 1.0.0, 30-JAN-2006, EDW (JPL)
%
%-Index_Entries
%
%   None.
%
%-&

function [y] = zzmice_dp(x)

   if( ~isequal(nargin,1) )

      error( 'MICE(USAGE): _y_ = zzmice_dp( _x_ )' )

   end

   if( isnumeric(x) && all( isfinite( x(:) ) ) )

     y = double(x);

   else

      error( ['MICE(BADARG): Improper type of input ' ...
              'argument passed to function. Value ' ...
              'or values expected as finite double precision.'] )

   end

author: goGPS-Project


text:
function roi = roiGrayAll(seg, mr, roi);
%
% roi = roiGrayAll(seg, [mr], [roi]);
%
% Create an ROI containing all gray matter nodes for a
% segmentation, or append these nodes to an existing ROI.
%
% seg: segmentation object. See segCreate.
%
% roi: if an ROI structure is passed as a second argument,
%      will append the gray coords to the roi, rather than
%      creating a new one.
%
% ras, 01/2007.
if notDefined('seg'), error('Need a segmentation.'); end
if notDefined('roi'), roi = roiCreate('I|P|R');      end
if notDefined('mr'),  mr = [];                       end

if isempty(seg.nodes)
    [seg.nodes seg.edges] = segGet(seg, 'gray');
end

roi.coords = combineCoords(roi.coords, seg.nodes([2 1 3],:), 'union');
roi.name = sprintf('Gray Matter %s', seg.name);

if ~isempty(mr)
    roi = roiCheckCoords(roi, mr);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function U = uCheck(U,TolMesh,optimState,proj)
%UCHECK Remove duplicates, vectors outside bounds or previously evaluated

if nargin < 4 || isempty(proj); proj = 0; end

if proj
    % Project vectors outside bounds on search mesh points closest to bounds
    U = bsxfun(@max, bsxfun(@min, U, optimState.UBsearch), optimState.LBsearch);
else
    idx = any(bsxfun(@gt,U,optimState.UB) | bsxfun(@lt,U,optimState.LB),2);
    U(idx,:) = [];
end

% Remove duplicate vectors
U = unique(U,'rows');

% Remove previously evaluted vectors (within TolMesh)
if ~isempty(U)
    TolMesh = TolMesh/2;
    % TolMesh = exp(0.5*(log(TolMesh) + log(optimState.meshsize)))/2;
    % Convert everything to normalized grid coordinates for comparison
    u1 = round(U./TolMesh);
    U2 = optimState.U(1:optimState.Xmax,:);
    u2 = round(U2./TolMesh);
    [~,idx] = setdiff(u1,u2,'rows');
    U = U(idx,:);
end

% Evaluate non-bound constraints
if ~isempty(optimState.nonbcon)
    X = origunits(U,optimState);    % Convert back to original space
    C = optimState.nonbcon(X);      % Evaluate constraints
    idx = (C <= 0);                 % Keep points that satisfy constraints
    U = U(idx,:);
end

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [varargout] = read_plexon_nex(filename, varargin)

% READ_PLEXON_NEX reads header or data from a Plexon *.nex file, which
% is a file containing action-potential (spike) timestamps and waveforms
% (spike channels), event timestamps (event channels), and continuous
% variable data (continuous A/D channels).
%
% LFP and spike waveform data that is returned by this function is
% expressed in microVolt.
%
% Use as
%   [hdr] = read_plexon_nex(filename)
%   [dat] = read_plexon_nex(filename, ...)
%   [dat1, dat2, dat3, hdr] = read_plexon_nex(filename, ...)
%
% Optional arguments should be specified in key-value pairs and can be
%   header      structure with header information
%   feedback    0 or 1
%   tsonly      0 or 1, read only the timestamps and not the waveforms
%   channel     number, or list of numbers (that will result in multiple outputs)
%   begsample   number (for continuous only)
%   endsample   number (for continuous only)
%
% See also READ_PLEXON_PLX, READ_PLEXON_DDT

% Copyright (C) 2007, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% parse the optional input arguments
hdr       = ft_getopt(varargin, 'header');
channel   = ft_getopt(varargin, 'channel');
feedback  = ft_getopt(varargin, 'feedback', false);
tsonly    = ft_getopt(varargin, 'tsonly', false);
begsample = ft_getopt(varargin, 'begsample', 1);
endsample = ft_getopt(varargin, 'endsample', inf);

% start with empty return values and empty data
varargout = {};

% read header info from file, use Matlabs for automatic byte-ordering
fid = fopen_or_error(filename, 'r', 'ieee-le');
fseek(fid, 0, 'eof');
siz = ftell(fid);
fseek(fid, 0, 'bof');

if isempty(hdr)
  if feedback, fprintf('reading header from %s\n', filename); end
  % a NEX file consists of a file header, followed by a number of variable headers
  % sizeof(NexFileHeader) = 544
  % sizeof(NexVarHeader) = 208
  hdr.FileHeader = NexFileHeader(fid);
  if hdr.FileHeader.NumVars<1
    ft_error('no channels present in file');
  end
  hdr.VarHeader = NexVarHeader(fid, hdr.FileHeader.NumVars);
end

for i=1:length(channel)
  chan = channel(i);
  vh = hdr.VarHeader(chan);
  clear buf
  fseek(fid, vh.DataOffset, 'bof');
  switch vh.Type
    case 0
      % Neurons, only timestamps
      buf.ts = fread(fid, [1 vh.Count], 'int32=>int32');

    case 1
      % Events, only timestamps
      buf.ts = fread(fid, [1 vh.Count], 'int32=>int32');

    case 2
      % Interval variables
      buf.begs = fread(fid, [1 vh.Count], 'int32=>int32');
      buf.ends = fread(fid, [1 vh.Count], 'int32=>int32');

    case 3
      % Waveform variables
      buf.ts = fread(fid, [1 vh.Count], 'int32=>int32');
      if ~tsonly
        buf.dat = fread(fid, [vh.NPointsWave vh.Count], 'int16');
        % convert the AD values to miliVolt, subsequently convert from miliVolt to microVolt
        buf.dat = buf.dat * (vh.ADtoMV * 1000);
      end

    case 4
      % Population vector
      ft_error('population vectors are not supported');

    case 5
      % Continuously recorded variables
      buf.ts   = fread(fid, [1 vh.Count], 'int32=>int32');
      buf.indx = fread(fid, [1 vh.Count], 'int32=>int32');
      if vh.Count>1 && (begsample~=1 || endsample~=inf)
        ft_error('reading selected samples from multiple AD segments is not supported');
      end
      if ~tsonly
        numsample = min(endsample - begsample + 1, vh.NPointsWave);
        fseek(fid, (begsample-1)*2, 'cof');
        buf.dat  = fread(fid, [1 numsample], 'int16');
        % convert the AD values to miliVolt, subsequently convert from miliVolt to microVolt
        buf.dat = buf.dat * (vh.ADtoMV * 1000);
      end

    case 6
      % Markers
      buf.ts = fread(fid, [1 vh.Count], 'int32=>int32');
      for j=1:vh.NMarkers
        buf.MarkerNames{j,1} = fread(fid, [1 64], 'uint8=>char');
        for k=1:vh.Count
          buf.MarkerValues{j,k} = fread(fid, [1 vh.MarkerLength], 'uint8=>char');
        end
      end

    otherwise
      ft_error('incorrect channel type');
  end % switch channel type

  % return the data of this channel
  varargout{i} = buf;
end % for channel

% always return the header as last
varargout{end+1} = hdr;

fclose(fid);
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function hdr = NexFileHeader(fid)
hdr.NexFileHeader  = fread(fid,4,'uint8=>char')';    % string NEX1
hdr.Version        = fread(fid,1,'int32');
hdr.Comment        = fread(fid,256,'uint8=>char')';
hdr.Frequency      = fread(fid,1,'double');         % timestamped freq. - tics per second
hdr.Beg            = fread(fid,1,'int32');          % usually 0
hdr.End            = fread(fid,1,'int32');          % maximum timestamp + 1
hdr.NumVars        = fread(fid,1,'int32');          % number of variables in the first batch
hdr.NextFileHeader = fread(fid,1,'int32');          % position of the next file header in the file, not implemented yet
Padding = fread(fid,256,'uint8=>char')';             % future expansion

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function hdr = NexVarHeader(fid, numvar)
for varlop=1:numvar
  hdr(varlop).Type         = fread(fid,1,'int32');        % 0 - neuron, 1 event, 2- interval, 3 - waveform, 4 - pop. vector, 5 - continuously recorded
  hdr(varlop).Version      = fread(fid,1,'int32');        % 100
  hdr(varlop).Name         = fread(fid,64,'uint8=>char')'; % variable name
  hdr(varlop).DataOffset   = fread(fid,1,'int32');        % where the data array for this variable is located in the file
  hdr(varlop).Count        = fread(fid,1,'int32');        % number of events, intervals, waveforms or weights
  hdr(varlop).WireNumber   = fread(fid,1,'int32');        % neuron only, not used now
  hdr(varlop).UnitNumber   = fread(fid,1,'int32');        % neuron only, not used now
  hdr(varlop).Gain         = fread(fid,1,'int32');        % neuron only, not used now
  hdr(varlop).Filter       = fread(fid,1,'int32');        % neuron only, not used now
  hdr(varlop).XPos         = fread(fid,1,'double');       % neuron only, electrode position in (0,100) range, used in 3D
  hdr(varlop).YPos         = fread(fid,1,'double');       % neuron only, electrode position in (0,100) range, used in 3D
  hdr(varlop).WFrequency   = fread(fid,1,'double');       % waveform and continuous vars only, w/f sampling frequency
  hdr(varlop).ADtoMV       = fread(fid,1,'double');       % waveform continuous vars only, coeff. to convert from A/D values to Millivolts
  hdr(varlop).NPointsWave  = fread(fid,1,'int32');        % waveform only, number of points in each wave
  hdr(varlop).NMarkers     = fread(fid,1,'int32');        % how many values are associated with each marker
  hdr(varlop).MarkerLength = fread(fid,1,'int32');        % how many characters are in each marker value
  Padding = fread(fid,68,'uint8=>char')';
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = remez(f, varargin)
%REMEZ   Best polynomial or rational approximation for real valued CHEBFUNs.
%   This code is obsolete as of May 2017.  Use MINIMAX instead.

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

warning('CHEBFUN:CHEBFUN:remez', ...
        [' This command is deprecated.', ...
         ' Use minimax instead.']);
polyOutput = detectType(f,varargin{:});

if ( polyOutput )
    [p,err,status] = minimax(f,varargin{:});
    varargout = {p, err, status};
else
    [p,q,r,err,status] = minimax(f,varargin{:});
    varargout = {p, q, r, err, status};
end

end

function polyOutput = detectType(~, varargin)

isSilent = 0;
for k = 1:length(varargin)
    if ( ischar(varargin{k}) && strcmpi('silent', varargin{k}) )
        isSilent = 1;
    end
end

% Detect polynomial / rational approximation type.
polyOutput = true;
if ( mod(nargin - isSilent, 2) ) % Odd number of inputs --> rational case.
    polyOutput = false;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function test_bug1483

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_selectdata ft_selectdata_new ft_selectdata_old

%% first confirm the reported bug, i.e. the incapability of ft_selectdata_new
% to select 'rpt'

% create some dummy data
for k = 1:5
  data.trial{k} = ones(2,3).*k;
  data.time{k}  = -1:1;
end
data.label = {'chan1';'chan2'};

freq.powspctrm = rand(5,2,3);
freq.dimord    = 'rpt_chan_freq';
freq.label     = {'chan1';'chan2'};
freq.freq      = [1 2 3];
freq.cumtapcnt = ones(5,1);

% old style -> this does not work anymore because ft_selectdata in its old
% implementation has been moved to compat/obsolete, and currently the input
% data checking is so strict, that a non cfg-like data structure causes a
% crash
%dataold1 = ft_selectdata(data, 'rpt', 2:4);
%dataold2 = ft_selectdata(data, 'channel', data.label(1));

%freqold1 = ft_selectdata(freq, 'rpt', 2:4);
%freqold2 = ft_selectdata(freq, 'channel', freq.label(1));

% new style
cfg = [];
cfg.trials = 2:4;
datanew1 = ft_selectdata(cfg, data);
freqnew1 = ft_selectdata(cfg, freq);

cfg = [];
cfg.channel = data.label(1);
datanew2 = ft_selectdata(cfg, data);
freqnew2 = ft_selectdata(cfg, freq);

% BUG CONFIRMED

%%


text:
function [perimeterEdges,eulerCondition] = findLegalPerimeters2(mesh,perimDist)
% Calculate edges that on distinct perimeters perimDist mm from startVertex
%
%  [perimeterEdges,eulerCondition] = findLegalPerimeters2(mesh,perimDist)
%
% Given a mesh structure (nodes, edges, distances from start point) and a
% threshold distance return a list of edges that constitute separate
% perimeters at a distance of perimDist from the start point.
%
% Note that because of intrinsic curvature, there can be more than one
% perimeter at the required distance. This routine makes sure that it
% returns separate perimeters (ones with no common nodes).
%

% Find perims with simple threshold
insideNodes = find(mesh.dist<=perimDist);
insideNodes = insideNodes(:);

% Some triangles have one or two vertices outside the perimeter.  These
% triangles will not be included; any nodes that are within the perimeter
% but only part of these excluded triangles are removed here.
insideNodes = removeHangingNodes2(mesh,insideNodes);

% We could write a short piece of code that verifies these new nodes have
% no hanging nodes ...

numBadNodes = 1;

while (numBadNodes>0)

    [perimeterEdges,eulerCondition] = findGroupPerimeter2(mesh,insideNodes);

    length(perimeterEdges);
    length(unique(perimeterEdges,'rows'));
    fprintf('Euler number=%d\n',eulerCondition);

    badPerimNodes = findBadPerimNodes2(mesh,perimeterEdges);
    numBadNodes   = length(badPerimNodes);
    fprintf('There are %d bad perim nodes.\n',numBadNodes);

    if(numBadNodes)
        % Splits up joined perimeters
        [insideNodes] = correctBadNodes2(mesh,insideNodes,badPerimNodes);

        % Cleans up mesh again
        insideNodes = removeHangingNodes2(mesh,insideNodes);
    end
end

return;


text:
%% (Internal) GUI for inspecting and correcting timeseries
%
% Description:
% This function implements a graphical user interface (GUI) to correct annotations
% possibly previous performed by an automatic algorithm. The idea is to
% easily visualize, compare and correct annotations. For this purpose the
% GUI presents several representations of the time evolution of the events
% (possibly but not necessarily heartbeats).
%
% Arguments:
%
%     +ECG: [numeric or cell]
%
%           [numeric]: signal matrix of dimension [sig_length sig_size
%           repetitions_size] where:
%             - sig_length: time length in samples
%             - sig_size: number of ECG leads or number of signals.
%             - repetitions_size: number of repetitions of the same
%             signals. Typically used when time-synchronized events, like
%             heartbeats.
%
%           [cell]: cell array of length repetitions_size, where each cell
%           is (probably a time alligned event) a signal of dimension
%           [sig_length sig_size]
%
%     +QRS_locations: [numeric] OPTIONAL. Default values enclosed in ()
%                 Synchronization sample. In ECG context, this values are
%                 the QRS fiducial point. (empty)
%
%
%     +ECG_header: [struct] OPTIONAL.
%
%             Description of the ECG typically available in the
%             header. Structure with fields:
%
%               -freq: Sampling rate in Hz. (1)
%
%               -nsig: Number of ECG leads. (size(ECG,2))
%
%               -nsamp: Number of ECG samples. (size(ECG,1))
%
%               -adczero: ADC offset (e.g. 2^(adc_res-1) when
%               using unsigned integers). ( repmat(0, ECG_header.nsig , 1) )
%
%               -adcgain: ADC gain in units/adc_sample
%               (typically uV/adc_unit). ( repmat(1, ECG_header.nsig , 1) )
%
%               -units: Measurement units. ( repmat('uV', ECG_header.nsig , 1) )
%
%               -desc: Signal description. ( num2str(colvec(1:ECG_header.nsig)) )
%
%     +QRS_annotations: [cell] OPTIONAL. Default values enclosed in ()
%               Annotations to be included in the mosaic. The funcion
%               accepts 2 type of annotations: points and lines.
%
% Example:
%
%    https://www.youtube.com/watch?v=qgWjvsvafVg&list=PLlD2eDv5CIe9sA2atmnb-DX48FIRG46z7&index=3
%
% See also ECGtask_QRS_corrector
%
% Limits and Known bugs:
%   Probably a lot :( ... but dont panic! send me feedback if you need help.
%
% Author: Mariano Llamedo Soria (llamedom at frba.utn.edu.ar)
% Version: 0.1 beta
% Birthdate  : 16/08/2017
% Last update: 16/08/2017
% Copyright 2008-2017
%
function ts_out = timeseries_explorer(tsc_in)

    tsc_out

    %% Global variables

    fig_hdl = figure();

    title_efimero_hdl = [];

    %%

    set(fig_hdl, 'WindowButtonDownFcn', @WindowButtonDownCallback);
    set(fig_hdl, 'WindowButtonUpFcn',   @WindowButtonUpCallback);
    set(fig_hdl, 'WindowKeyPressFcn',   @WindowKeyPress);

    maximize(fig_hdl);
    maximized_size = get(fig_hdl, 'Position');

    set(fig_hdl, 'Position', [ maximized_size(3:4) maximized_size(3:4) ] .* [ 0.05 0.13 0.95 0.9] );

    set(fig_hdl,'CloseRequestFcn',@my_closefcn)

    user_data.RRserie = colvec(diff(user_data.anns_under_edition));
    user_data.RRserie = [user_data.RRserie(1); user_data.RRserie] * 1/user_data.ECG_struct.header.freq;

    limits = prctile(user_data.RRserie, [5 95]);

    if(limits(1) > 0)
        limits(1) = limits(1) * 0.9;
    else
        limits(1) = limits(1) * 1.1;
    end
    limits(2) = limits(2) * 1.1;

    if( user_data.bLockScatter )
        x_lims_scatter = get(user_data.Scatter_axes_hdl, 'Xlim');
        y_lims_scatter = get(user_data.Scatter_axes_hdl, 'Ylim');
    end

    if( user_data.bLockRRserie )
        x_lims_RRserie = get(user_data.RRserie_axes_hdl, 'Xlim');
        y_lims_RRserie = get(user_data.RRserie_axes_hdl, 'Ylim');
    end

    user_data.fig_hdl = figure(1);
%     set(user_data.fig_hdl, 'Toolbar', 'none');
    clf

    %% Axis de Poincar�

    user_data.Scatter_axes_hdl = axes('Position', [0.55 0.41 0.355 0.52 ]);

    RRscatter_hdl = scatter(user_data.Scatter_axes_hdl, user_data.RRserie, [ user_data.RRserie(2:end); user_data.RRserie(end) ] );

    if( ~isempty(user_data.selected_hb_idx) )
        hold(user_data.Scatter_axes_hdl, 'on')
        RRnext = [user_data.RRserie(2:end); user_data.RRserie(end)];
        user_data.RRscatter_selection_hdl = scatter(user_data.Scatter_axes_hdl, user_data.RRserie(user_data.selected_hb_idx), RRnext(user_data.selected_hb_idx), 'xg' );
        hold(user_data.Scatter_axes_hdl, 'off')
    end

    if( user_data.bLockScatter )
        xlim(user_data.Scatter_axes_hdl, x_lims_scatter);
        ylim(user_data.Scatter_axes_hdl, y_lims_scatter);
    else
        xlim(user_data.Scatter_axes_hdl, limits);
        ylim(user_data.Scatter_axes_hdl, limits);
        zoom reset
    end

    title(user_data.Scatter_axes_hdl, ['Poincar� plot - Registro ' user_data.ECG_struct.header.recname ]);
    xlabel(user_data.Scatter_axes_hdl, 'RR current')
    ylabel(user_data.Scatter_axes_hdl, 'RR next')

    %% Axis de RR serie
    user_data.RRserie_axes_hdl = axes('Position', [0.13 0.709 0.375 0.216 ]);

    RRserie_hdl = plot(user_data.RRserie_axes_hdl, user_data.RRserie, '-xb' );

    if( ~isempty(user_data.selected_hb_idx) )
        hold(user_data.RRserie_axes_hdl, 'on')
        user_data.RRserie_selection_hdl = plot(user_data.RRserie_axes_hdl, user_data.selected_hb_idx, user_data.RRserie(user_data.selected_hb_idx), 'og' );
        hold(user_data.RRserie_axes_hdl, 'off')
    end

    if( user_data.bLockRRserie )
        xlim(user_data.RRserie_axes_hdl, x_lims_RRserie);
        ylim(user_data.RRserie_axes_hdl, y_lims_RRserie);
    else
        ylim(user_data.RRserie_axes_hdl, limits);
        zoom reset
    end

    cant_ticks = 5;
    aux_idx = round(linspace(1, length(user_data.RRserie), cant_ticks));
    set(user_data.RRserie_axes_hdl, 'XTick', rowvec(aux_idx) );
    aux_str = strcat(cellstr(num2str(colvec(aux_idx))), cellstr(repmat(' (',cant_ticks,1)), cellstr(Seconds2HMS(user_data.anns_under_edition(aux_idx)*1/user_data.ECG_struct.header.freq)), cellstr(repmat(')',cant_ticks,1)) );    set(user_data.RRserie_axes_hdl, 'XTickLabel', char(aux_str) );
    set(user_data.RRserie_axes_hdl, 'XTickLabel', char(aux_str) );
    xlabel(user_data.RRserie_axes_hdl, 'heartbeat #')
    ylabel(user_data.RRserie_axes_hdl, 'RR interval')

%% global functions

    function update_title_efimero( strTitle, delay )

        if( ~isempty(title_efimero_hdl) && ~isempty(ancestor(title_efimero_hdl,'figure')) && gcf == ancestor(title_efimero_hdl,'figure') && ishandle(title_efimero_hdl) )
            set(title_efimero_hdl, 'String',strTitle )
            set(title_efimero_hdl, 'Visible', 'on');
        else

            title_efimero_hdl = annotation( 'textbox', [0 0.98 0.3 0.01 ], ...
                                  'String', strTitle, ...
                                  'Tag', 'title_efimero', ...
                                  'FontSize', 8, ...
                                  'Interpreter', 'none', ...
                                  'FitBoxToText', 'on', ...
                                  'HorizontalAlignment', 'left', ...
                                  'EdgeColor', 'none', ...
                                  'BackgroundColor', 'r' );

        end

        if( ~isinf(delay) && strcmpi(my_timer.Running, 'off') )
            my_timer.StartDelay = delay;
            start(my_timer)
        end

    end

    function timer_fcn(obj,event_obj)

        set(title_efimero_hdl, 'Visible', 'off');

    end

%% Window callbacks

    function WindowButtonDownCallback(obj,event_obj)

        if (strcmp(get(fig_hdl,'SelectionType'),'alt'))

            prev_u = get(fig_hdl, 'units');
            set(fig_hdl, 'units','normalized');
            crd = get(fig_hdl, 'CurrentPoint');
            xp = crd(1);
            yp = crd(2);
            set(fig_hdl, 'units', prev_u);

            if( ~isempty([min_y_drag max_x_drag min_y_drag]) && xp <= max_x_drag && yp <= max_y_drag && yp >= min_y_drag)
                DragMouseBegin()
            end
        end

    end

    function WindowButtonUpCallback(obj,event_obj)

        DragMouseEnd()

    end

    function KeyPress(obj,event_obj)

        if (strcmp(event_obj.Key,'c') && strcmp(event_obj.Modifier,'control'))
            % copy selection
            if( ~isempty(selected_hb_idx) )
                copy_paste_buffer = anns_under_edition(selected_hb_idx);
                update_title_efimero('Selection copied', 5 );

            end

        elseif (strcmp(event_obj.Key,'v') && strcmp(event_obj.Modifier,'control'))
            % paste selection
            if( ~isempty(copy_paste_buffer) )
                bAnnsEdited = true;
                bRecEdited = true;
                PushUndoAction();
                anns_under_edition = sort( unique([anns_under_edition; colvec(copy_paste_buffer) ]) );
                selected_hb_idx = [];
                RRserie_filt = {[]};

                Redraw();
            end

        elseif (strcmp(event_obj.Key,'delete'))
            % delete selection
            if( ~isempty(selected_hb_idx) )
                bAnnsEdited = true;
                bRecEdited = true;
                PushUndoAction();

                aux_idx = [find(anns_under_edition < anns_under_edition(anns_under_edition_idx(selected_hb_idx(1))),1, 'last') find( anns_under_edition > anns_under_edition(anns_under_edition_idx(selected_hb_idx(end))), 1, 'first') ];

                if( isempty(aux_idx) )
                    aux_idx = 1;
                else
                    aux_idx = aux_idx(1);
                end

                aux_val = anns_under_edition(aux_idx);

                if( bSeries )
                    anns_under_edition(anns_under_edition_idx(selected_hb_idx)) = nan;
                else
                    anns_under_edition(anns_under_edition_idx(selected_hb_idx)) = [];
                end

                anns_under_edition_idx = find( anns_under_edition >= start_idx &  anns_under_edition <= end_idx );
                hb_idx = find(anns_under_edition(anns_under_edition_idx) == aux_val);
                selected_hb_idx = [];

                aux_rr_filt = RRserie_filt{1};
                aux_rr_filt(anns_under_edition_idx(selected_hb_idx)) = [];
                RRserie_filt{1} = aux_rr_filt;

                Redraw();
            end

        elseif (strcmp(event_obj.Key,'y') && ~isempty(event_obj.Modifier) && strcmp(event_obj.Modifier,'control'))
            % redo last change
            if( undo_buffer_idx < length(undo_buffer) )
                undo_buffer_idx = undo_buffer_idx + 1;
                anns_under_edition = undo_buffer{undo_buffer_idx};
                selected_hb_idx = [];
                RRserie_filt = {[]};
                Redraw();
            end

        elseif (strcmp(event_obj.Key,'z') && strcmp(event_obj.Modifier,'control'))
            % undo last change
            if( undo_buffer_idx > 1 )
                undo_buffer_idx = undo_buffer_idx - 1;
                anns_under_edition = undo_buffer{undo_buffer_idx};
                selected_hb_idx = [];
                RRserie_filt = {[]};

                Redraw();
            end

        elseif (strcmp(event_obj.Key,'rightarrow') && ~isempty(event_obj.Modifier) && strcmp(event_obj.Modifier,'control'))

            if(bRecEdited)

                update_annotations();

                if( bLoadECG )
                    update_title_efimero('Saving data ...', 5 );
                    save(rec_path, '-struct', 'ECG_struct');
                    update_title_efimero(['Saved ' rec_path], 5 );

                    if( rec_idx >= 1 && rec_idx < length(recording_indexes) )
                        rec_idx = rec_idx + 1;
                    else
                        rec_idx = 1;
                    end

                    set(recordings_control, 'Value', rec_idx );

                    DoRecording();

                else

                    if( isfield(ECG_struct, 'series_quality' ) )
                        anns_struct.series_quality = ECG_struct.series_quality;
                    end
                    for ii = 1:size(AnnNames,1)
                        anns_struct.(AnnNames{ii,1}) = ECG_struct.(AnnNames{ii,1});
                    end
                    assignin( 'caller', OutputVarName, anns_struct );
                    update_title_efimero(sprintf('Saving ''%s'' variable in caller workspace', OutputVarName), 5 );

                    bAnnsEdited = false;
                    bRecEdited = false;

                end
            end

        elseif (strcmp(event_obj.Key,'leftarrow') && ~isempty(event_obj.Modifier) && strcmp(event_obj.Modifier,'control'))

            if(bRecEdited)

                update_annotations();

                if( bLoadECG )
                    update_title_efimero('Saving data ...', 5 );
                    save(rec_path, '-struct', 'ECG_struct');
                    update_title_efimero(['Saved ' rec_path], 5 );

                    if( rec_idx > 1 && rec_idx <= length(recording_indexes) )
                        rec_idx = rec_idx - 1;
                    else
                        rec_idx = length(recording_indexes);
                    end

                    set(recordings_control, 'Value', rec_idx );

                    DoRecording();

                else
                    if( isfield(ECG_struct, 'series_quality' ) )
                        anns_struct.series_quality = ECG_struct.series_quality;
                    end
                    for ii = 1:size(AnnNames,1)
                        anns_struct.(AnnNames{ii,1}) = ECG_struct.(AnnNames{ii,1});
                    end
                    assignin( 'caller', OutputVarName, anns_struct );
                    update_title_efimero(sprintf('Saving ''%s'' variable in caller workspace', OutputVarName), 5 );

                    bAnnsEdited = false;
                    bRecEdited = false;

                end
            end

        elseif (strcmp(event_obj.Key, 'l') )
            % lock ECG Y axis

            bLockECG = ~bLockECG;

            if( bLockECG )
                update_title_efimero('ECG Y axis locked', 5 );
            else
                update_title_efimero('ECG Y axis unlocked', 5 );
            end

        elseif (strcmp(event_obj.Key,'p') )

            update_title_efimero('Click and drag the time interval where the pattern is ...', 10 );

            set(fig_hdl, 'WindowButtonDownFcn', []);
            set(ECG_axes_hdl,'ButtonDownFcn', []);

            set(obj,'CurrentAxes', ECG_axes_hdl);
            waitforbuttonpress;

            point1 = get(ECG_axes_hdl,'CurrentPoint');    % button down detected
            rbbox;
            point2 = get(ECG_axes_hdl,'CurrentPoint');    % button up detected

            pattern_match_xlims = round(sort([ point1(1,1) point2(1,1) ]));
            aux_seq = pattern_match_xlims(1):pattern_match_xlims(2);

            update_title_efimero('Filtering ECG ...', 5 );

            if( isempty(filtro) )
                ECG = ECG_struct.signal(:,lead_idx);
            else
                ECG = filter(filtro, flipud(ECG_struct.signal(:,lead_idx)) );
                ECG = filter(filtro, flipud(ECG) );
            end

            if( ishandle(Pattern_hdl) )
                delete(Pattern_hdl)
            end

            hold(ECG_axes_hdl, 'on')
            Pattern_hdl = plot(ECG_axes_hdl, aux_seq, ECG(aux_seq,:), '.g' );
            hold(ECG_axes_hdl, 'off')
            drawnow;

            try
                SearchPattern();
                set(fig_hdl, 'WindowButtonDownFcn', @WindowButtonDownCallback2D);
                set(ECG_axes_hdl,'ButtonDownFcn',@inspect_ECG);

            catch ME

                set(fig_hdl, 'WindowButtonDownFcn', @WindowButtonDownCallback2D);
                set(ECG_axes_hdl,'ButtonDownFcn',@inspect_ECG);
                rethrow(ME);
            end

            % store the axes_hdl created in SearchPattern();
            axes_hdl( RR_global_PM_k, 1:2) = [RR_global_PM_hdl, figPatternMatch_hdl];
            axes_hdl( proximity_k, 1:2) = [proximity_hdl, figPatternMatch_hdl];
            axes_hdl( similarity_hdl_k, 1:2) = [similarity_hdl, figPatternMatch_hdl];
            axes_hdl( pattMatch_hdl_k, 1:2) = [pattMatch_hdl, figPatternMatch_hdl];

        elseif ( ~isempty(event_obj.Modifier) && strcmp(event_obj.Key,'g') && strcmp(event_obj.Modifier,'control'))

            if(bRecEdited)

                update_annotations();

                if( bLoadECG )
                    update_title_efimero('Saving data ...', 5 );
                    save(rec_path, '-struct', 'ECG_struct');
                    update_title_efimero(['Saved ' rec_path], 5 );
                else

                    if( isfield(ECG_struct, 'series_quality' ) )
                        anns_struct.series_quality = ECG_struct.series_quality;
                    end
                    for ii = 1:size(AnnNames,1)
                        anns_struct.(AnnNames{ii,1}) = ECG_struct.(AnnNames{ii,1});
                    end
                    assignin( 'caller', OutputVarName, anns_struct );
                    update_title_efimero(sprintf('Saving ''%s'' variable in caller workspace', OutputVarName), 5 );

                    bAnnsEdited = false;
                    bRecEdited = false;
                end

            end

        elseif (strcmp(event_obj.Key,'t'))
            %% Toggle ECG lead
            if( length(lead_idx) > 1 )
                lead_idx = 1;
            else
                if( lead_idx < ECG_struct.header.nsig )
                    lead_idx = lead_idx + 1;
                else
                    lead_idx = 1:ECG_struct.header.nsig;
                end
            end
            cla(ECG_axes_hdl);

            if( bSeries )
                this_all_anns = all_annotations_selected_serie_location;
                aux_val = this_all_anns{1};
                aux_val(serie_location_mask) = nan;
                this_all_anns{1} = aux_val;
            else
                this_all_anns = all_annotations_selected;
            end

            if( ~bLockECG )
                ECG_limits = [];
            end

            [ECG_hdl, ECG_limits ] = plot_ecg_heartbeat(ECG_struct.signal, lead_idx, this_all_anns, start_idx, anns_under_edition_idx(hb_idx) , hb_detail_window, ECG_struct.header, filtro, ECG_axes_hdl, ECG_limits);

            if( length(lead_idx) > 1 )
                aux_str = rowvec(colvec([repmat(',', length(lead_idx), 1) ECG_struct.header.desc(lead_idx,:) ]'));
                title(ECG_axes_hdl, ['Heartbeat ' num2str(hb_idx) ' : Leads ' aux_str(2:end) ] )
            else
                title(ECG_axes_hdl, ['Heartbeat ' num2str(hb_idx) ' : Lead ' ECG_struct.header.desc(lead_idx,:)] )
            end

    %         zoom reset
            cellfun(@(a)( set(a,'ButtonDownFcn',@inspect_ECG)), ECG_hdl);
            set(ECG_axes_hdl,'ButtonDownFcn',@inspect_ECG);

        elseif (strcmp(event_obj.Key,'rightarrow'))
            %% right arrow

            if(bAnnsEdited)
                update_annotations();
            end

            AnnNames_idx = AnnNames_idx + 1;
            if( AnnNames_idx > length(AnnNames) )
                AnnNames_idx = 1;
            end

            update_title_efimero(['Using ' AnnNames{AnnNames_idx,1} ' annotations (' num2str(ratios(AnnNames_idx)) ')'], 5 );

            undo_buffer_idx = 1;

            anns_under_edition = unique(round(colvec( ECG_struct.(AnnNames{AnnNames_idx,1}).(AnnNames{AnnNames_idx,2}) )));

            bAnnsEdited = false;

            selected_hb_idx = [];

            Redraw();

        elseif (strcmp(event_obj.Key,'leftarrow'))
            %% left arrow

            if(bAnnsEdited)
                update_annotations();
            end

            AnnNames_idx = AnnNames_idx - 1;
            if( AnnNames_idx < 1  )
                AnnNames_idx = length(AnnNames);
            end

            update_title_efimero(['Using ' AnnNames{AnnNames_idx,1} ' annotations (' num2str(ratios(AnnNames_idx)) ')'], 5 );

            undo_buffer_idx = 1;

            anns_under_edition = unique(round(colvec( ECG_struct.(AnnNames{AnnNames_idx,1}).(AnnNames{AnnNames_idx,2}) )));

            bAnnsEdited = false;

            selected_hb_idx = [];

            Redraw();

        else
%             bLockRRserie = false;
%             bLockScatter = false;
%             lead_idx = 1:ECG_struct.header.nsig;
        end

    end

end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function per_init(num_pkts, snr_vec)
% per_graph.m - draws PER graph
global numPkts PER_snr;

% PER graph settings (actual values are user-selected (GUI))
numPkts = num_pkts;
PER_snr = snr_vec;


text:
function [L,U,p,q,inform] = luSOL(A)

%        [L,U,p,q,inform] = luSOL(A);
%
% This is an example script that uses lusolFactor
% to find LU factors of sparse matrix A
% using predefined options (hardwired here).

% In general, it is better to call lusolSet yourself
% and then reset a few options if necessary as shown.

% 28 Apr 2004: First version of luSOL based in splu.m.
%              Michael Saunders, SOL, Stanford University.

options = lusolSet;
options.Pivoting   = 'TRP';
options.FactorTol  = 4.0;

[L,U,p,q,options]  = lusolFactor(A,options);

inform = options.Inform;


text:
function ndx = sub2indFast(siz, varargin)
% ndx = sub2indFast(siz, varargin)
%
% Like sub2ind but without time-intensive checks.
%
% Copyright by Mathworks
% Modified by Holger Caesar, 2014

assert(all(numel(siz) == numel(varargin)));

%Compute linear indices
k = [1, cumprod(siz(1:end-1))];
ndx = 1;
numOfIndInput = nargin-1;
for i = 1:numOfIndInput
    v = varargin{i};
    ndx = ndx + (double(v)-1) * k(i);
end

author: nightrome

repo: matconvnet-calvin


text:
function vizLineRpts(ax,data,fld1,fld2,line)
% vizLineRpts(ax,data,fld1,fld2,line)
% Visualize experimental repeats for given stat/line.
%
% data: data struct
% fld1: field for x-axis of scatterplot. Defaults to day_num.
% fld2: field for y-axis of scatterplot
% line: line to visualize

if isempty(ax)
    figure;
    ax = axes;
end

if isempty(fld1)
    fld1 = 'day_num';
end

tfLine = strcmp({data.line_name}',line);
tmpGF = {data.bdpOrNot}';
%tmpGF = repmat({'not'},numel(tfLine),1);
tmpGF(tfLine) = {line};

axes(ax);
x = OlyDat.Analysis.getNumericScalarFieldSafe(data,fld1);
y = OlyDat.Analysis.getNumericScalarFieldSafe(data,fld2);
h = gscatter(x,y,tmpGF);
set(h(1),'MarkerEdgeColor',[0.75 0.75 0.75]);
set(h(2),'MarkerEdgeColor',[0.3 0.3 0.8]);
set(h(3),'MarkerEdgeColor',[1 0 0],'MarkerSize',14);
hLeg = legend;
set(hLeg,'interpreter','none');
xlabel(fld1,'interpreter','none','FontWeight','bold');
ylabel(fld2,'interpreter','none','FontWeight','bold');
titlestr = sprintf('%s: %d exps',line,nnz(tfLine));
title(titlestr,'interpreter','none','FontWeight','bold');
grid on;

% plot ctl mean, std in x- and y- dirs
hold on;
tfCtl = [data.tfBDP]';
xCtl = [data(tfCtl).(fld1)];
yCtl = [data(tfCtl).(fld2)];
mux = mean(xCtl);
sdx = std(xCtl);
lox = mux-2*sdx;
hix = mux+2*sdx;
muy = mean(yCtl);
sdy = std(yCtl);
loy = muy-2*sdy;
hiy = muy+2*sdy;
%plot(mux,muy,'ob','MarkerSize',14,'MarkerFaceColor',[0 0 1]);
% plot([lox hix],[loy loy],'b');
% plot([lox hix],[hiy hiy],'b');
% plot([lox lox],[loy hiy],'b');
% plot([hix hix],[loy hiy],'b');
hold off;

author: kristinbranson

repo: JAABA


text:
function imin = i4vec_imin ( n, a )

%*****************************************************************************80
%
%% I4VEC_IMIN computes the index of the minimum element of an I4VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input, integer A(N), the array.
%
%    Output, integer IMIN, the index of the smallest entry.
%
  if ( n <= 0 )

    imin = 0;

  else

    amin = a(1);
    imin = 1;

    for i = 2 : n

      if ( a(i) < amin )
        amin = a(i);
        imin = i;
      end

    end

  end

  return
end


text:
function value = p32_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P32_F evaluates the integrand for problem 32.
%
%  Dimension:
%
%    N arbitrary.
%
%  Region:
%
%    0 <= X(1:DIM_NUM) <= 1
%
%  Integral Parameters:
%
%    The integral depends on vectors C(1:N) and Z(1:N).
%
%    The reference suggests choosing C at random in [0,1]
%    and then multiplying by the normalizing factor sqrt(140/N**(3/2)).
%
%    The default value of C(1:N) is (1/2)^(1/N).
%
%    The default value of Z(1:N) is (1/2)^(1/N).
%
%  Integrand:
%
%    exp ( c(1:n)*x(1:n) ) if all x(1:n) <= z(1:n)
%    0                        otherwise
%
%  Exact Integral:
%
%    product ( g(1:n)(x,z,a,b,c) )
%
%    where g(i)(x,z,a,b,c) =
%
%      0                                         if z(i) <= a(i)
%      ( e^(c(i)*z(i) ) - e^(c(i)*a(i)) ) / c(i) if a(i) <= z(i) <= b(i)
%      ( e^(c(i)*b(i) ) - e^(c(i)*a(i)) ) / c(i) if b(i) <= z(i)
%
%    with obvious modifications when c(i) = 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Alan Genz,
%    [Integral #6]
%    A Package for Testing Multiple Integration Subroutines,
%    in Numerical Integration: Recent Developments, Software
%    and Applications,
%    edited by Patrick Keast and Graeme Fairweather,
%    D Reidel, 1987, pages 337-340,
%    LC: QA299.3.N38.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  c = [];
  c = p32_r8vec ( 'G', 'C', dim_num, c );

  z = [];
  z = p32_r8vec ( 'G', 'Z', dim_num, z );

  value(1:point_num) = 0.0;

  for point = 1 : point_num

    if ( all ( x(1:dim_num,point) <= z(1:dim_num)' ) )
      value(point) = exp ( c(1:dim_num) * x(1:dim_num,point) );
    else
      value(point) = 0.0;
    end

  end

  p32_i4 ( 'I', '#', point_num );

  return
end


text:
function x = norm_ball( x, varargin ) %#ok

%NORM_BALL   Norm ball.
%   NORM_BALL( sz, ... ) returns a variable of size sz, say 'x', that is
%   constrained to satisfy NORM( x, ... ) <= 1. Any syntactically valid
%   and _convex_ use of the NORM() function has a direct analog in
%   NORM_BALL. The convex requirement specifically excludes, then, all
%   instances of NORM( x, p ) where p < 1.
%
%   See NORM for more detaills.
%
%   Disciplined convex programming information:
%       NORM_BALL is a cvx set specification. See the user guide for
%       details on how to use sets.

cvx_begin set
    variable x( size(x) )
    norm( x, varargin{:} ) <= 1; %#ok
cvx_end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook


function varths = Get_VarthsPSI_swing(zetaj,ntilj,PSIdj,PSIdjm1)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here

zetaj2     = zetaj^2; zetaj3 = zetaj*zetaj; zetaj4 = zetaj*zetaj3;

varths = zeros(1,4);
gamms  = zeros(1,3);

gamms(1)  = (1 - zetaj4)/8 + zetaj3/3 - zetaj2/4;
gamms(2)  = 5/12 + zetaj4/4 - zetaj3/3 - zetaj2/2 + zetaj;
gamms(3)  = 1/12 - (1 + zetaj4)/8 + zetaj2/4;

varths(3) = PSIdjm1(ntilj(1)-1)*gamms(1) + PSIdjm1(ntilj(1))*gamms(2) + PSIdjm1(ntilj(1)+1)*gamms(3);

gamms(1) = 1/12 - gamms(1);
gamms(2) = 5/6  - gamms(2);
gamms(3) = 1/12 - gamms(3);

varths(1) = PSIdj(ntilj(1)-1)*gamms(1) + PSIdj(ntilj(1))*gamms(2) + PSIdj(ntilj(1)+1)*gamms(3);

%%%----------------------

gamms(1)  = zetaj4/8 - zetaj3/2 + zetaj2/2;
gamms(2)  = -zetaj4/4 + 2*zetaj3/3;
gamms(3)  = zetaj4/8 - zetaj3/6;

varths(4) = PSIdjm1(ntilj(1))*gamms(1) + PSIdjm1(ntilj(1)+1)*gamms(2) + PSIdjm1(ntilj(1)+2)*gamms(3);

gamms(1) = 1/12 - gamms(1);
gamms(2) = 5/6  - gamms(2);
gamms(3) = 1/12 - gamms(3);

varths(2) = PSIdj(ntilj(1))*gamms(1) + PSIdj(ntilj(1)+1)*gamms(2) + PSIdj(ntilj(1)+2)*gamms(3);

end


text:
function d = vgg_H_sampson_distance_sqr(H,X1,X2)

% d = vgg_H_sampson_distance_sqr(H,X1,X2)
%
% Returns an approximation of the squared geometric distance from the
% 4d joint-space point [X1;X2] to the H manifold. See Torr CVIU 97.

if (size(X1) ~= size(X2))
  error('Point sets not same size!');
end

p1 = X1 ./ repmat(X1(3,:),3,1);
p2 = X2 ./ repmat(X2(3,:),3,1);

alg = vgg_H_algebraic_distance(H,p1,p2);

N = size(X1,2);

h = reshape(H',9,1);

G1 = [ H(1,1) - p2(1,:) * H(3,1) ; ...
  H(1,2) - p2(1,:) * H(3,2) ; ...
  -p1(1,:) * H(3,1) - p1(2,:) * H(3,2) - H(3,3) ; ...
  zeros(1,N) ];

G2 = [ H(2,1) - p2(2,:) * H(3,1) ; ...
  H(2,2) - p2(2,:) * H(3,2) ; ...
  zeros(1,N) ; ...
  -p1(1,:) * H(3,1) - p1(2,:) * H(3,2) - H(3,3) ];

magG1 = sqrt(sum(G1 .* G1));
magG2 = sqrt(sum(G2 .* G2));
magG1G2 = sum(G1 .*  G2);

alpha = acos( magG1G2 ./ (magG1 .* magG2) );

D1 = alg(1,:) ./ magG1;
D2 = alg(2,:) ./ magG2;

d = (D1.*D1 + D2.*D2 - 2 * D1 .* D2 .* cos(alpha)) ./ sin(alpha);

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function pz=v_lpcpp2pz(pp)
%V_LPCPP2PZ LPC: Convert power spectrum polynomial in cos(w) to power spectrum zeros PZ=(RP)
% pp is a polynomial such that |polyval(ra,e^jw)| = polyval(pp,cos(w))

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: v_lpcpp2pz.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pz=roots(pp);

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function Y = LSA(X, nLowVec)

k = nLowVec;

[Y,~,~] = svds(X,k);

end

author: jacoxu


text:
function Yout = prtUtilEvalCAPtree(tree,X,nClasses)
%Yout = evalCAPtree(tree,X)
%   Evaluate a CAP tree on a 1xN data point X
% Internal
% xxx Need Help xxx

if nargin < 3 || isempty(nClasses)
    nClasses = 2;
end

index = 1;
voted = false;
while ~voted
    if any(isfinite(tree.W(:,index)))
        %disp(((tree.W(:,index)'*X(:,tree.featureIndices(:,index))')') - tree.threshold(:,index))
        Yout = double(((tree.W(:,index)'*X(:,tree.featureIndices(:,index))')' - tree.threshold(:,index)) >= 0);
        if Yout == 0
            index = find(tree.treeIndices(:) == index,1,'first');
        elseif Yout > 0
            index = find(tree.treeIndices(:) == index,1,'last');
        end
    else
        Yout = zeros(1,nClasses);
        Yout(tree.terminalVote(index)) = 1;
        voted = true;
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function pde = fracLapdata1
%% FRACLAPDATA1 data for fractional Laplacian problem
%
% s a parameter in (0,1)
% k is given integer
% u = 2^(1-s)/gamma(s)*(k*pi*y)^s*besselk(s,k*pi*y).*sin(k*pi*x);
% u|_(y==0) =  sin(k*pi*x);
% g_N = ds*(k*pi)^(2*s)*sin(k*pi*x);
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

pde = struct('f',0,'exactu',@exactu,'g_D',0,'g_N',@g_N,'d',@d,'Du',@Du);

    % exact solution
    function u =  exactu(p)
        global s
        u = zeros(size(p,1),1);
        k = 3;
        x = p(:,1); y = p(:,2);
        idx = (y>eps);
        sqrtlambda = k*pi;
        z = sqrtlambda*y(idx);
        C = 2^(1-s)/gamma(s);
        u(idx) =  C*z.^s.*besselk(s,z).*sin(sqrtlambda*x(idx));
        u(~idx) = sin(sqrtlambda*x(~idx));  % at y = 0.  %% NOS 2013 p.31
    end
    function K =  d(p)  % Diffusion constant
        global s
        alpha = 1-2*s;
        %x = p(:,1);
        y = p(:,2);
        K = y.^alpha;
    end
    function z = g_N(p)
        global s
        k = 3;
        x = p(:,1); %y = p(:,2);
        ds = 2^(1-2*s)*gamma(1-s)/gamma(s);  % (2.23) in NOS 2013
        sqrtlambda = k*pi;
        z = ds*sqrtlambda^(2*s)*sin(sqrtlambda*x);
    end
    function uprime =  Du(p)
        global s
        k = 3;
        x = p(:,1); y = p(:,2);
        sqrtlambda = k*pi;
        z = sqrtlambda*y;
        uprime(:,1) = 2^(1-s)/gamma(s)*sqrtlambda*z.^s.*besselk(s,z).*cos(sqrtlambda*x);
        uprime(:,2) = -2^(1-s)/gamma(s)*sqrtlambda*z.^s.*besselk(1-s,z).*sin(sqrtlambda*x);
    end
end

repo: Aero-matlab


text:
function [F, Srez, Scel]=statindexes(XData, YData, a, b)
%
% [F, Srez, Scel]=statindexes(xdata, ydata, a, b)
%
% Computation of statistical indicators for linear model
%
% Input parameters:
%  - xdata: input data block -- x: axis
%  - ydata: input data block -- y: axis
%  - a: parameter of linear model \  y = bx + a
%  - b: parameter of linear model /
%
% Return parameters:
%  - F: value for F-test
%  - Srez: residual dispersion
%  - Scel: total dispersion
%
% Authors:
% Ivo Petras (ivo.petras@tuke.sk)
% Dagmar Bednarova (dagmar.bednarova@tuke.sk)
%
% Date: 25/05/2009
%
kx=length(XData);
ky=length(YData);;
if kx ~= ky
   disp('Incompatible X and Y data.');
   close all;
end
n=size(YData,2);
Srez = sum((YData-(a+b*XData)/(sqrt(1+b^2))).^2)/(n-2);
Scel=sum((YData-mean(YData)).^2+(XData-mean(XData)).^2)/(n-1);
F=Scel/Srez;
%

author: Sable

repo: mcbench-benchmarks


text:
function [f]= spm_fx_adem_pursuit(x,v,a,P)
% returns the flow for occulomotor pursuit (with action)
% FORMAT [f]= spm_fx_adem_pursuit(x,v,a,P)
%
% x    - hidden states:
%   x.o(1) - oculomotor angle
%   x.o(2) - oculomotor angle
%   x.x(1) - target location (visual) - extrinsic coordinates (Cartesian)
%   x.x(2) - target location (visual) - extrinsic coordinates (Cartesian)
%   x.a(:) - attractor (SHC) states
%
% v    - hidden cause (speed)
% P    - parameters
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_fx_adem_pursuit.m 4625 2012-01-24 20:53:10Z karl $

% intisaise flow (to ensure fields are aligned)
%--------------------------------------------------------------------------
f    = x;

% motion of attractor states
%==========================================================================
f.a  = spm_lotka_volterra(x.a,v);

% motion of target states
%==========================================================================

% target location is determined by the attractor state x.a
%--------------------------------------------------------------------------
t    = P*spm_softmax(x.a,1/2);
f.x  = (t - x.x)/2;

% motion of oculomotor angles (driven by bounded action)
%==========================================================================
f.o  = tanh(a) - x.o/16;

author: spm


text:
function hsvMean = selectPixelsAndGetHSV(RGB, Area)

%
% function hsvMean = selectPixelsAndGetHSV(RGB, Area)
%
% Use this function in order to select multiple points from an image (use
% right click to stop process). The selected points are used to calculate
% the average HSV values.

% ARGUMENTS:
% RGB: the RGB image
% Area: the area size used to calulate the HSV values of each point
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Theodoros Giannakopoulos - January 2008
% www.di.uoa.gr/~tyiannak
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clf;
warning off;
imshow(RGB); hold on;
HSV = rgb2hsv(RGB);
HSV2 = HSV;
numOfSelectedPixels = 0;
right_not_pressed = 1;
BUTTON = 1;
while (BUTTON~=3)
    numOfSelectedPixels = numOfSelectedPixels + 1;
    [X,Y,BUTTON] = GINPUT(1);

    hsvTemp2 = HSV(Y-(Area-1)/2:Y+(Area-1)/2, X-(Area-1)/2:X+(Area-1)/2, :);

    HSV2(Y-round((Area-1)/2):Y+round((Area-1)/2), X-round((Area-1)/2):X+round((Area-1)/2), :) = 0;

    hsvTemp = zeros(3,1);
    [K,L,M] = size(hsvTemp2);
    for (i=1:K)
        for (j=1:L)
            hsvTemp(1) = hsvTemp(1) + hsvTemp2(i,j,1);
            hsvTemp(2) = hsvTemp(2) + hsvTemp2(i,j,2);
            hsvTemp(3) = hsvTemp(3) + hsvTemp2(i,j,3);
        end
    end

    hsvTemp = hsvTemp / (K*L);

    hsv(numOfSelectedPixels,:) = hsvTemp;
    hsvMean = median(hsv,1);

    line([X-(Area-1)/2 X+(Area-1)/2] , [Y-(Area-1)/2 Y-(Area-1)/2]);
    line([X+(Area-1)/2 X+(Area-1)/2] , [Y-(Area-1)/2 Y+(Area-1)/2]);
    line([X+(Area-1)/2 X-(Area-1)/2] , [Y+(Area-1)/2 Y+(Area-1)/2]);
    line([X-(Area-1)/2 X-(Area-1)/2] , [Y+(Area-1)/2 Y-(Area-1)/2]);

    %rgbTemp = hsv2rgb(hsvTemp);
    %fprintf('Cur RGV Values:  %.3f %.3f %.3f\n', rgbTemp(1), rgbTemp(2), rgbTemp(3));
    %fprintf('Cur HSV Values:  %.3f %.3f %.3f\n', hsvTemp(1), hsvTemp(2), hsvTemp(3));
    %fprintf('Mean HSV Values: %.3f %.3f %.3f\n', hsvMean(1), hsvMean(2), hsvMean(3));
end

[N, t] = size(hsv);

% for (i=1:N) hsvM1(i) = median(hsv(1:i,1)); end
% for (i=1:N) hsvM2(i) = median(hsv(1:i,2)); end
% for (i=1:N) hsvM3(i) = median(hsv(1:i,3)); end

hsvMean = median(hsv);

% figure;
% subplot(3,1,1); plot(hsvM1); title(sprintf('H-->%.4f',hsvM1(end)));
% subplot(3,1,2); plot(hsvM2); title(sprintf('S-->%.4f',hsvM2(end)));
% subplot(3,1,3); plot(hsvM3); title(sprintf('V-->%.4f',hsvM3(end)));
%
% figure;
% RGB2 = hsv2rgb(HSV2);
% imshow(RGB2);


text:
classdef PTKROIDividedIntoLeftAndRight < PTKPlugin
    % PTKROIDividedIntoLeftAndRight. Plugin to divide the roi region into left and right
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Left and <br>Right ROI'
        ToolTip = 'Separate and label left and right lungs'
        Category = 'Lungs'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
    end

    methods (Static)
        function results = RunPlugin(dataset, reporting)
            results = dataset.GetResult('PTKLeftAndRightLungs');
            [~, nn_index] = bwdist(results.RawImage > 0);
            results_raw = results.RawImage(nn_index);
            results.ChangeRawImage(results_raw);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function short_sample()
% Short syntax highlighting sample
global gbl_value;
gbl_value = 20
hghlghtvar = 'a_string_example_here';

%% Title of cell
badstring = 'unterminated_string;
unusedvar = hghlghtvar;
!echo "Example `!` system command"
for iterval = 1:2:gbl_value
  % This is a very informative comment
  fprintf('%s:%d\n',hghlghtvar,iterval);
end
end

author: scottclowe

repo: matlab-schemer


text:
%% Coarsening for Algebraic Multigrid
%
% Given a SPD matrix A, we describe an algebraic coarsening of a graph of A
% based on the concept of strong connectness. The measure of strong
% connectness is slightly different with the standard definition. The
% parameter theta is used to define strong connectness and the default
% value is 0.025.

%% Usage of the function
clear all
help coarsenAMGc

%% Generate a test matrix
[node,elem] = squaremesh([0,1,0,1],1/8);
% [node,elem] = uniformrefine(node,elem);
% load lakemesh
[A,M] = assemblematrix(node,elem);
% A = M;  % test mass matrix. No coarsening is needed.

%% Parameters
theta = 0.025;
N = size(A,1);
N0 = min(sqrt(N),50);       % number of the coarest nodes

%% Generate strong connectness matrix
D = spdiags(1./sqrt(diag(A)),0,N,N);
Am = D*A*D;  % normalize diagonal of A
[im,jm,sm] = find(Am);
idx = (-sm > theta);   % delete weakly connect off-diagonal and diagonal
As = sparse(im(idx),jm(idx),sm(idx),N,N); % matrix for strong connectness
% The diagonal of Am is 1. The negative off-diagonal measures the
% diffusivity. The positive off-diagonal is filtered.

%% Compute degree of vertex
deg = sum(spones(As)); % number of strongly connected neighbors
deg = full(deg');
% deg = deg + rand(N,1); % break the equal degree case but deteriorate performance
if sum(deg>0) < 0.1*sqrt(N)   % too few connected nodes e.g. A is mass matrix
    isC(round(rand(N0,1)*N)) = true; % randomly chose N0 nodes
    return                    % smoother is a good preconditioner
end
idx = (deg>0);
deg(idx) = deg(idx) + 0.1*rand(sum(idx),1); % break the equal degree case

%% Find an approximate maximal independent set and put to C set
isC = false(N,1);       % C: coarse node
isF = false(N,1);       % F: fine node
isU = true(N,1);        % S: selected set
isF(deg == 0) = true;   % isolated nodes are added into F set
% debug
close all;
%%
% * magneta dots: indepedent nodes in U
% * yellow dots: F (fine) nodes
% * red dots: C (coarse) nodes
% * black dots: U (undecided) nodes
set(gcf,'Units','normal');
set(gcf,'Position',[0.5,0.5,0.5,0.5]);
showmesh(node,elem);
findnode(node,isU,'noindex','Color','k','MarkerSize',32)
m = 1;
while sum(isC) < N/2 && sum(isU) >N0
    % Mark all undecided nodes
    isS = false(N,1);  % S: selected set, changing in the coarsening
    isS(deg>0) = true;
    S = find(isS);
    % debug
    fprintf('Coarsening ... \n');

    % Find marked nodes with local maximum degree
    [i,j] = find(triu(As(S,S),1));    % i,j and i<j: edges of subgraph S
    idx = deg(S(i)) >= deg(S(j));     % compare degree of vertices
    isS(S(j(idx))) = false;  % remove vertices with smaller degree
    isS(S(i(~idx))) = false;
    isC(isS) = true;
    findnode(node,isS,'noindex','Color','m');
    fprintf('Number of chosen points: %6.0u\n',sum(isS));
    snapnow

    % Remove coarse nodes and neighboring nodes from undecided set
    [i,j] = find(As(:,isC)); %#ok<*NASGU>
    isF(i) = true;        % neighbor of C nodes are F nodes
    isU = ~(isF | isC);   % U: undecided set
    deg(~isU) = 0;        % remove current C and F from the graph
    % -- No improvement by adding weight to nodes connected to F nodes --
%     degFin = sum(spones(As(isF,isU)));
%     degFin = full(degFin'); % degrer of strong connected with F points
%     deg(isU) = deg(isU) + degFin; % add weight to nodes connected to F

    if sum(isU) <= N0   % add small undecided nodes into C nodes
        isC(isU) = true;
        isU = [];       % to exit the while loop;
    end
    % debug
    showmesh(node,elem);
    findnode(node,isU,'noindex','Color','k','MarkerSize',32)
    findnode(node,isF,'noindex','Color','y','MarkerSize',32)
    findnode(node,isC,'noindex','Color','r','MarkerSize',36);
    snapnow
    m = m + 1;
end
fprintf('Apply %2.0u times and Number of Coarse Nodes: %6.0u\n',m,sum(isC));

%%
% Note that the red nodes are connected in the grid but not connected in
% the graph of the 5-point stencil.


text:
function HDR = mat2sel(FileName,fnout);
% MAT2SEL is useful for converting visual artifact selection
%   with gBSanalyze into BKR/EPS compliant *.SEL files.
%
%   HDR = mat2sel(infile, [selfile]);
%      infile   input filename (*.BKR or *.MAT are supported)
%      selfile  output file with info artifact selection [OPTIONAL]
%               if no selfile is provided, the information will be
%               stored with extension '.SEL' and the same filename and
%               path.
%
%   The artifact selection is available in HDR.ArtifactSelection and
%   can be also loaded with
%       [s,HDR]=sload(filename);
%   Currently, artifact selection of  EPS/BKR-Software and
%   gBSanalyze is supported.
%
% see also: SLOAD

%	$Revision: 1.1 $
%	$Id: mat2sel.m 2205 2009-10-27 12:18:15Z schloegl $
%	Copyright (C) 2004 by Alois Schloegl
%	a.schloegl@ieee.org

% This library is free software; you can redistribute it and/or
% modify it under the terms of the GNU Library General Public
% License as published by the Free Software Foundation; either
% Version 2 of the License, or (at your option) any later version.
%
% This library 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
% Library General Public License for more details.
%
% You should have received a copy of the GNU Library General Public
% License along with this library; if not, write to the
% Free Software Foundation, Inc., 59 Temple Place - Suite 330,
% Boston, MA  02111-1307, USA.

[s,HDR] = sload(FileName);

if isfield(HDR,'ArtifactSelection'),
        artifact = HDR.ArtifactSelection;

        save([HDR.FILE.Name, '_artifact'], 'artifact');

        if nargin<2,
                fnout = [HDR.FILE.Name, '.sel']
                fnout = fullfile(HDR.FILE.Path,[HDR.FILE.Name, '.sel']);
        end
        if 0, exist(fnout)==2,
                fprintf(1,'File %s exists! Do you want to overwrite ',fnout);
                answer = input('[Y/N] ?');
        else
                answer = 'Y';
        end;
        if any(answer=='jJyYzZ'),
                fid = fopen(fnout,'w');
                fprintf(fid, '%i\r\n', artifact);
                fclose(fid);
        end;
end;

author: PatternRecognition

repo: OpenBMI


text:
function test_suite = test_survival_weibull

%   Run specific demo and save values for comparison.
%
%   See also
%     TEST_ALL, DEMO_SURVIVAL_WEIBULL

% Copyright (c) 2011-2012 Ville Tolvanen

initTestSuite;

function testDemo
% Set random number stream so that failing isn't because randomness. Run
% demo & save test values.
prevstream=setrandstream(0);

disp('Running: demo_survival_weibull')
demo_survival_weibull;
path = which('test_survival_weibull.m');
path = strrep(path,'test_survival_weibull.m', 'testValues');
if ~(exist(path, 'dir') == 7)
    mkdir(path)
end
path = strcat(path, '/testSurvival_weibull');
save(path, 'Ef1', 'Ef2', 'Varf1', 'Varf2');

% Set back initial random stream
setrandstream(prevstream);
drawnow;clear;close all

% Compare test values to real values.

function testPredictionsWeibull
values.real = load('realValuesSurvival_weibull', 'Ef1', 'Varf1', 'Ef2', 'Varf2');
values.test = load(strrep(which('test_survival_weibull.m'), 'test_survival_weibull.m', 'testValues/testSurvival_weibull'), 'Ef1', 'Varf1', 'Ef2', 'Varf2');
assertElementsAlmostEqual(values.real.Ef1, values.test.Ef1, 'relative', 0.10);
assertElementsAlmostEqual(values.real.Ef2, values.test.Ef2, 'relative', 0.10);
assertElementsAlmostEqual(values.real.Varf1, values.test.Varf1, 'relative', 0.10);
assertElementsAlmostEqual(values.real.Varf2, values.test.Varf2, 'relative', 0.10);