text:
function fn_handle = GiveMeFunctionHandle(cfnParams)
% GiveMeFunctionHandle  Returns a function handle for computing classification accuracy
%
%---INPUTS:
% cfnParams, parameters for the classification, e.g., set with
%           GiveMeDefaultClassificationParams.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

% Set empty function handle:
fn_handle = [];

%-------------------------------------------------------------------------------
% Set the function handle to compute the accuracy/loss measure:
%-------------------------------------------------------------------------------
if strcmp(cfnParams.whatClassifier,'fast-linear')
    fn_loss = @(yTest,yPredict) BF_LossFunction(yTest,yPredict,cfnParams.whatLoss,cfnParams.classLabels);
    fn_handle = @(XTrain,yTrain,XTest,yTest) fn_loss(yTest,classify(XTest,XTrain,yTrain,'linear'));
    return
end

%-------------------------------------------------------------------------------
% Binary model (easier), we can do it in one line:
%-------------------------------------------------------------------------------
if cfnParams.numClasses==2
    % Set the loss function:
    fn_loss = @(yTest,yPredict) BF_LossFunction(yTest,yPredict,cfnParams.whatLoss,cfnParams.classLabels);

    switch cfnParams.whatClassifier
    case 'knn'
        fn_handle = @(XTrain,yTrain,XTest,yTest) fn_loss(yTest,predict(fitcknn(XTrain,yTrain),XTest));
    case 'tree'
        fn_handle = @(XTrain,yTrain,XTest,yTest) fn_loss(yTest,predict(fitctree(XTrain,yTrain),XTest));
    case {'linear','linclass'}
        if cfnParams.doReweight
            fn_handle = @(XTrain,yTrain,XTest,yTest) ...
                        fn_loss(yTest,predict(fitcdiscr(XTrain,yTrain,'FillCoeffs','off',...
                                    'SaveMemory','on','Weights',InverseProbWeight(yTrain)),XTest));
        else
            fn_handle = @(XTrain,yTrain,XTest,yTest) ...
                        fn_loss(yTest,predict(fitcdiscr(XTrain,yTrain,'FillCoeffs','off',...
                                    'SaveMemory','on'),XTest));
        end
    case {'svm','svm-linear','svm-linear-highdim'}
        % This one uses the fitclinear function (fast; suited for high-dimensional datasets)
        if cfnParams.doReweight
            fn_handle = @(XTrain,yTrain,XTest,yTest) fn_loss(yTest,predict(fitclinear(XTrain,yTrain,...
                                    'KernelFunction','linear','Weights',InverseProbWeight(yTrain)),XTest));
        else
            fn_handle = @(XTrain,yTrain,XTest,yTest) fn_loss(yTest,predict(fitclinear(XTrain,yTrain,...
                                    'KernelFunction','linear'),XTest));
        end
    case {'svm-linear-lowdim'}
        % This one uses the fitcsvm function (slower; suited for low-dimensional datasets)
        if cfnParams.doReweight
            fn_handle = @(XTrain,yTrain,XTest,yTest) fn_loss(yTest,predict(fitcsvm(XTrain,yTrain,...
                                    'KernelFunction','linear','Weights',InverseProbWeight(yTrain)),XTest));
        else
            fn_handle = @(XTrain,yTrain,XTest,yTest) fn_loss(yTest,predict(fitcsvm(XTrain,yTrain,...
                                    'KernelFunction','linear'),XTest));
        end
    case 'svm-rbf'
        if cfnParams.doReweight
            fn_handle = @(XTrain,yTrain,XTest,yTest) fn_loss(yTest,predict(fitcsvm(XTrain,yTrain,...
                                    'KernelFunction','rbf','Weights',InverseProbWeight(yTrain)),XTest));
        else
            fn_handle = @(XTrain,yTrain,XTest,yTest) fn_loss(yTest,predict(fitcsvm(XTrain,yTrain,...
                                    'KernelFunction','rbf'),XTest));
        end
    case 'diaglinear'
        if cfnParams.doReweight
            fn_handle = @(XTrain,yTrain,XTest,yTest) fn_loss(yTest,predict(fitcnb(XTrain,yTrain,...
                                                    'Weights',InverseProbWeight(yTrain)),XTest));
        else
            fn_handle = @(XTrain,yTrain,XTest,yTest) fn_loss(yTest,predict(fitcnb(XTrain,yTrain),XTest));
        end
    otherwise
        warning('Unknown classifier: ''%s''',cfnParams.whatClassifier);
    end
end

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

if isempty(fn_handle)
    % Multiclass, have to use the Cfn function (which is slower to use as a function handle)
    % OR: unknown classifier set above
    fn_handle = @(XTrain,yTrain,XTest,yTest) ...
        GiveMeCfn(XTrain,yTrain,XTest,yTest,cfnParams,false);
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% This demo efficiently computes emd_hat between two random SIFTs
% histograms (3d histograms, see Lowe's paper: "Distinctive image
% features from scale-invariant keypoints" for more detail) where the
% ground distance between the bins is the thresholded sum of orientation and
% spatial distanced.
% emd_hat was described in the paper:
% A Linear Time Histogram Metric for Improved SIFT Matching
% Ofir Pele, Michael Werman
%  ECCV 2008
% The efficient algorithm is described in the paper:
%  Fast and Robust Earth Mover's Distances
%  Ofir Pele, Michael Werman
%  ICCV 2009

clc; close all; clear all;
rand('state',sum(100*clock));

XNBP= 4; % SIFT's X-dimension
YNBP= 4; % SIFT's Y-dimension
NBO= 16; % SIFT's Orientation-dimension
N= XNBP*YNBP*NBO;
thresh= 2;
extra_mass_penalty= -1; % Default of maximum distance
flowType= 3; % Regular flows

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SIFT ground distance matrix computation.
% Note: loops in Matlab are costly.
% However, this should be done only once.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D= zeros(N,N);
i= 1;
for y1=1:YNBP
    for x1=1:XNBP
        for nbo1=1:NBO

            j= 1;
            for y2=1:YNBP
                for x2=1:XNBP
                    for nbo2=1:NBO
                        D(i,j)= (sqrt((y1-y2)^2 + (x1-x2)^2) + ...
                                 min( [abs(nbo1-nbo2) NBO-abs(nbo1-nbo2)] ));

                        j=j+1;
                    end
                end
            end
            i= i+1;

        end
    end
end
maxDist= max(D(:));
D= min(D,thresh);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

P= rand(N,1);
Q= rand(N,1);

% The demo includes several ways to call emd_hat_mex and emd_hat_gd_metric_mex
demo_FastEMD_compute(P,Q,D,extra_mass_penalty,flowType);

% Copyright (c) 2009-2012, Ofir Pele
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%    * Redistributions of source code must retain the above copyright
%    notice, this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright
%    notice, this list of conditions and the following disclaimer in the
%    documentation and/or other materials provided with the distribution.
%    * Neither the name of the The Hebrew University of Jerusalem nor the
%    names of its contributors may be used to endorse or promote products
%    derived from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
function epoch = mydatesowi (sow, epoch2)
    bow2 = mydatesow_aux (epoch2);
    epoch = sow + bow2;
end

%!test
%! % on a Sunday, sow == sod.
%! vec = [2008 11 30  13 14 15];
%! num = mydatenum(vec);
%! sod = mydatesod(num);
%! sow = mydatesow(num);
%! num2 = mydatesodi(sod, num);
%! num3 = mydatesowi(sow, num);
%! %sod, sow  % DEBUG
%! %num, num2, num3  % DEBUG
%! myassert (sow, sod, -sqrt(eps));
%! myassert (num2, num, -sqrt(eps));
%! myassert (num3, num, -sqrt(eps));

%!test
%! % random sow (within one weeek):
%! sow = randint(0, 7*24*3600);
%! bow_vec = [2008 11 30  0 0 0];
%! bow_num = mydatenum(bow_vec);
%! num = bow_num + sow;
%! sow2 = mydatesow(num);
%! num2 = mydatesowi(sow2, num);
%! %sow2, sow, sow2-sow  % DEBUG
%! %num2, num, num2-num  % DEBUG
%! tol = 1e-7;
%! myassert (sow2, sow, -tol);
%! myassert (num2, num, -tol);

%!test
%! % epoch2 <> epoch.
%! sow = randint(0, 7*24*3600);
%! bow_vec = [2008 11 30  0 0 0];
%! bow_num = mydatenum(bow_vec);
%! num = bow_num + sow;
%! bow_vec2 = [2008 11 30-7  0 0 0];
%! bow_num2 = mydatenum(bow_vec2);
%! sow2a = mydatesow(num, bow_num);
%! sow2b = mydatesow(num, bow_num2);
%! num2a = mydatesowi(sow2a, bow_num);
%! num2b = mydatesowi(sow2b, bow_num2);
%! %num, num2a, num2b  % DEBUG
%! tol = 1e-7;
%! myassert (sow2a, sow, -tol);
%! myassert (sow2b, sow+7*24*3600, -tol);
%! myassert (num2a, num, -tol);
%! myassert (num2b, num, -tol);

author: Sable

repo: mcbench-benchmarks


text:
function [tagName, hexString] = tag2strings(tag, dcmobj)
%"tag2strings"
%   Convert a decimal representation of a Dicom tag into two strings, the
%   tagName which is a decription of the field, and a hexString which gives
%   the tag in hex form '(xxxx,yyyy)'.
%
%   A dcmobj must be passed in to provide a link to the dicomdictionary.
%
%JRA 06/13/06
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Usage:
%   [tagName, hexString] = tag2strings(tag, dcmobj)

hexString    = char(org.dcm4che2.util.TagUtils.toString(tag));
tagName      = char(dcmobj.nameOf(tag));

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [opts, imdb] = model_setup(varargin)

% Copyright (C) 2015 Tsung-Yu Lin, Aruni RoyChowdhury, Subhransu Maji.
% All rights reserved.
%
% This file is part of the BCNN and is made available under
% the terms of the BSD license (see the COPYING file).

% setup ;

opts.seed = 1 ;
opts.batchSize = 128 ;
opts.numEpochs = 100;
opts.momentum = 0.9;
opts.learningRate = 0.001;
opts.numSubBatches = 1;
opts.keepAspect = true;
opts.useVal = false;
opts.fromScratch = false;
opts.useGpu = 1 ;
opts.regionBorder = 0.05 ;
opts.numDCNNWords = 64 ;
opts.numDSIFTWords = 256 ;
opts.numSamplesPerWord = 1000 ;
opts.printDatasetInfo = false ;
opts.excludeDifficult = true ;
opts.datasetSize = inf;
% opts.encoders = {struct('name', 'rcnn', 'opts', {})} ;
opts.encoders = {} ;
opts.dataset = 'cub' ;
opts.csiqDir = 'dataset\CSIQ database';
opts.liveDir = 'dataset\databaserelease2';
opts.tidDir = 'dataset\TID2013';
opts.cliveDir = 'dataset\ChallengeDB_release';
opts.mliveDir = 'dataset\LIVEmultidistortiondatabase\To_Release';
opts.ilsvrcDir = '/home/tsungyulin/dataset/ILSVRC2014/CLS-LOC/';
opts.ilsvrcDir_224 = '/home/tsungyulin/dataset/ILSVRC2014/CLS-LOC-224/';
opts.suffix = 'baseline' ;
opts.prefix = 'v1' ;
opts.model  = 'imagenet-vgg-m.mat';
opts.modela = 'imagenet-vgg-m.mat';
opts.modelb = [];
opts.layer  = 14;
opts.layera = [];
opts.layerb = [];
opts.imgScale = 1;
opts.bcnnLRinit = false;
opts.bcnnLayer = 14;
opts.rgbJitter = false;
opts.dataAugmentation = {'none', 'none', 'none'};
opts.cudnn = true;
opts.nonftbcnnDir = 'nonftbcnn';
opts.batchNormalization = false;
opts.cudnnWorkspaceLimit = 1024*1024*1204;
opts.dataset_path = 'dataset\databaserelease2';

[opts, varargin] = vl_argparse(opts,varargin) ;

opts.expDir = sprintf('data/%s/%s-seed-%02d', opts.prefix, opts.dataset, opts.seed) ;
opts.nonftbcnnDir = fullfile(opts.expDir, opts.nonftbcnnDir);
opts.imdbDir = fullfile(opts.expDir, 'imdb') ;
opts.resultPath = fullfile(opts.expDir, sprintf('result-%s.mat', opts.suffix)) ;

opts = vl_argparse(opts,varargin) ;

if nargout <= 1, return ; end

% % Setup GPU if needed
% if opts.useGpu
%   gpuDevice(opts.useGpu) ;
% end

% -------------------------------------------------------------------------
%                                                            Setup encoders
% -------------------------------------------------------------------------

models = {} ;
modelPath = {};
for i = 1:numel(opts.encoders)
  if isstruct(opts.encoders{i})
    name = opts.encoders{i}.name ;
    opts.encoders{i}.path = fullfile(opts.expDir, [name '-encoder.mat']) ;
    opts.encoders{i}.codePath = fullfile(opts.expDir, [name '-codes.mat']) ;
    [md, mdpath] = get_cnn_model_from_encoder_opts(opts.encoders{i});
    models = horzcat(models, md) ;
    modelPath = horzcat(modelPath, mdpath);
%     models = horzcat(models, get_cnn_model_from_encoder_opts(opts.encoders{i})) ;
  else
    for j = 1:numel(opts.encoders{i})
      name = opts.encoders{i}{j}.name ;
      opts.encoders{i}{j}.path = fullfile(opts.expDir, [name '-encoder.mat']) ;
      opts.encoders{i}{j}.codePath = fullfile(opts.expDir, [name '-codes.mat']) ;
      [md, mdpath] = get_cnn_model_from_encoder_opts(opts.encoders{i}{j});
      models = horzcat(models, md) ;
      modelPath = horzcat(modelPath, mdpath);
%       models = horzcat(models, get_cnn_model_from_encoder_opts(opts.encoders{i}{j})) ;
    end
  end
end

% -------------------------------------------------------------------------
%                                                       Download CNN models
% -------------------------------------------------------------------------

for i = 1:numel(models)
    if ~exist(modelPath{i})
        error(['cannot find model ', models{i}]) ;
    end
end

% -------------------------------------------------------------------------
%                                                              Load dataset
% -------------------------------------------------------------------------

vl_xmkdir(opts.expDir) ;
vl_xmkdir(opts.imdbDir) ;

switch opts.dataset
    case 'live'
        opts.dataset_path = opts.liveDir;
    case 'clive'
        opts.dataset_path = opts.cliveDir;
    case 'csiq'
        opts.dataset_path = opts.csiqDir;
    case 'tid'
        opts.dataset_path = opts.tidDir;
    case 'mlive'
        opts.dataset_path = opts.mliveDir;
    otherwise
        error('Unknown dataset %s', opts.dataset) ;
end

imdbPath = fullfile(opts.imdbDir, sprintf('imdb-seed-%d.mat', opts.seed)) ;
if exist(imdbPath)
  imdb = load(imdbPath) ;
  if(opts.rgbJitter)
      opts.pca = imdb_compute_pca(imdb, opts.expDir);
  end
  return ;
end

switch opts.dataset
    case 'live'
        imdb = getLiveDatabase_train(opts.liveDir);
    case 'clive'
        imdb = getCliveDatabase_train(opts.cliveDir);
    case 'csiq'
        imdb = getCsiqDatabase_train(opts.csiqDir);
    case 'tid'
        imdb = getTIDDatabase_train(opts.tidDir);
    case 'mlive'
        imdb = getMLIVEDatabase_train(opts.mliveDir);
    otherwise
        error('Unknown dataset %s', opts.dataset) ;
end

save(imdbPath, '-struct', 'imdb') ;

if(opts.rgbJitter)
   opts.pca = imdb_compute_pca(imdb, opts.expDir);
end

if opts.printDatasetInfo
  print_dataset_info(imdb) ;
end

% -------------------------------------------------------------------------
function [model, modelPath] = get_cnn_model_from_encoder_opts(encoder)
% -------------------------------------------------------------------------
p = find(strcmp('model', encoder.opts)) ;
if ~isempty(p)
  [~,m,e] = fileparts(encoder.opts{p+1}) ;
  model = {[m e]} ;
  modelPath = encoder.opts{p+1};
else
  model = {} ;
  modelPath = {};
end

% bilinear cnn models
p = find(strcmp('modela', encoder.opts)) ;
if ~isempty(p)
  [~,m,e] = fileparts(encoder.opts{p+1}) ;
  model = horzcat(model,{[m e]}) ;
  modelPath = horzcat(modelPath, encoder.opts{p+1});
end
p = find(strcmp('modelb', encoder.opts)) ;
if ~isempty(p) && ~isempty(encoder.opts{p+1})
  [~,m,e] = fileparts(encoder.opts{p+1}) ;
  model = horzcat(model,{[m e]}) ;
  modelPath = horzcat(modelPath, encoder.opts{p+1});
end

repo: DBCNN


text:
function varargout = drawVector3d(pos, vect, varargin)
%DRAWVECTOR3D Draw vector at a given position
%
%   drawVector3d(POS, VECT)
%   Draws the vector VECT starting at the position POS. Both VECT and POS
%   are N-by-3 arrays.
%
%   drawVector3d(..., PNAME, PVALUE)
%   Specifies additional optional parameters that will be given to the
%   quiver3 function.
%
%   Example
%     figure; hold on;
%     drawVector3d([2 3 4], [1 0 0]);
%     drawVector3d([2 3 4], [0 1 0]);
%     drawVector3d([2 3 4], [0 0 1]);
%     view(3);
%
%   See also
%   vectors3d, quiver3
%

% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2011-12-19,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

h = quiver3(pos(:, 1), pos(:, 2), pos(:, 3), ...
    vect(:, 1), vect(:, 2), vect(:, 3), 0, varargin{:});

% format output
if nargout > 0
    varargout{1} = h;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% determine searching area
active_pixel = (sum(neuron_bk.A(:,IDs), 2)>0);

% update spatial/temporal components of the merged neuron
data = neuron_bk.A(active_pixel, IDs)*neuron_bk.C_raw(IDs, :);
ci = neuron_bk.C_raw(IDs(1), :);
for miter=1:10
    ai = data*ci'/(ci*ci');
    ci = ai'*data/(ai'*ai);
end

sn = GetSn(ci);
neuron.A(active_pixel, IDs(1)) = ai*sn;
neuron.C_raw(IDs(1), :) = ci/sn;
neuron.S(IDs(1), :) = ci/sn;
newIDs(IDs(1)) = IDs(1);
% remove merged elements
ind_del(IDs(2:end)) = true;

%% show merged shape and trace
axes(ax_merged);
neuron.image(neuron.A(:, IDs(1)));
axis equal off tight;
set(gca, 'xlim', get(ax_selected, 'xlim'));
set(gca, 'ylim', get(ax_selected, 'ylim'));

axes(ax_merged_trace);
plot(ci, 'k', 'linewidth', 2);

%% clear IDs;
IDs_last = IDs;
IDs = [];

author: zhoupc


text:
classdef SimplePeakDetector < Algorithm

    properties (Access = public)
        minPeakHeight = single(0.8);
        minPeakDistance  = int32(100);
    end

    methods  (Access = public)

        function obj = SimplePeakDetector(minPeakHeight, minPeakDistance)
            if nargin > 0
                obj.minPeakHeight = minPeakHeight;
                obj.minPeakDistance = minPeakDistance;
            end
            obj.name = 'simplePeakDet';
            obj.inputPort = DataType.kSignal;
            obj.outputPort = DataType.kEvent;
        end

        function str = toString(obj)
            str = sprintf('%s_%.3f_%d',obj.name,obj.minPeakHeight,obj.minPeakDistance);
        end

        function editableProperties = getEditableProperties(obj)
            minPeakHeightProperty = Property('minPeakHeight',obj.minPeakHeight,0,1);
            minPeakDistanceProperty = Property('minPeakDistance',obj.minPeakDistance,50,200);
            editableProperties = [minPeakHeightProperty,minPeakDistanceProperty];
        end

        function events = compute(obj,signal)
            eventLocations = obj.detectEvents(signal);
            events = Event.EventsArrayWithEventLocations(eventLocations);
        end

        function metrics = computeMetrics(~,input)
            flops = size(input,1) * 4;
            memory = 1;
            outputSize = 0;
            metrics = Metric(flops,memory,outputSize);
        end
    end

    methods (Access = private)

        function peakLocations = detectEvents(obj,signal)

            signal = single(signal);
            lastPeakLocation = -1000;
            lastPeak = 0;

            [peakLocations, ~,~] = obj.findPeakLocations(signal,obj.minPeakHeight,obj.minPeakDistance,lastPeakLocation,lastPeak);
            peakLocations = uint32(peakLocations);
        end

        function [peakLocations, lastPeakLocation, lastPeakValue] = findPeakLocations(~,magnitude,...
                minPeakHeight, minPeakDistance, lastPeakLocation, lastPeakValue)
            numPeakLocations = uint8(0);
            peakLocations = int32(zeros(1,floor(length(magnitude) / minPeakDistance)));

            for i = 1 : length(magnitude)
                sample = magnitude(i);

                if sample >= minPeakHeight
                    if sample > lastPeakValue || i >= (lastPeakLocation + minPeakDistance)
                        lastPeakValue = sample;
                        lastPeakLocation = int32(i);
                    end
                end

                if lastPeakValue > single(0.00001) && i >= lastPeakLocation + minPeakDistance
                    lastPeakValue = single(0);
                    numPeakLocations = numPeakLocations + 1;
                    peakLocations(numPeakLocations) = lastPeakLocation;
                end
            end
            peakLocations = peakLocations(1:numPeakLocations);
        end
    end

end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function [ y, symm ] = cvx_s_symmetric_ut( m, n, symm )
%CVX_S_SYMMETRIC_UT Symmetric matrices (upper triangle storage).
if m ~= n,
    error( 'Symmetric structure requires square matrices.' );
end
symm = false;
nsq = n * n;
ntr = 0.5 * ( nsq + n );
c   = 0 : n - 1;
c   = c( ones( 1, n ), : );
r   = c';
mn  = min( r, c );
mx  = max( r, c );
y   = mn + 0.5 * mx .* ( mx + 1 ) + 1;
y   = sparse( y( : ), 1 : nsq, 1, ntr, nsq );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function con = spm_dcm_connectivity_ui(DCM,D,title_text,defaults,enabled)
% GUI for manually specifying connection values in a DCM
% FORMAT [con] = spm_dcm_connectivity_ui(DCM,D,title_text,defaults,enabled)
%
% DCM        - DCM structure
% D          - 'A','B' or 'C' i.e. connectivity matrix of interest
% title_text - Text to display above the matrix, e.g. 'Enter contrast for '
% defaults   - (optional) structure of default values containing
%              defaults.A, defaults.B and defaults.C
% enabled    - (optional) structure of inputs to enable with binary
%              matrices enabled.A, enabled.B and enabled.C
%
% Returns:
% con        - structure with con.A, con.B and con.C of user-entered values
%__________________________________________________________________________
% Copyright (C) 2014-2016 Wellcome Trust Centre for Neuroimaging

% Will Penny & Peter Zeidman
% $Id: spm_dcm_connectivity_ui.m 6808 2016-06-13 16:48:30Z guillaume $

% Set-up data
%--------------------------------------------------------------------------
n       = DCM.n;              % number of regions
m       = length(DCM.U.name); % number of inputs

region_names = DCM.Y.name;
input_names  = DCM.U.name;

try, A = defaults.A; catch, A = zeros(n,n); end
try, B = defaults.B; catch, B = zeros(n,n,m); end
try, C = defaults.C; catch, C = zeros(n,m); end

if nargin < 5
    enabled = struct('A',ones(size(A)),'B',ones(size(B)),'C',ones(size(C)));
end

% Set-up GUI
%--------------------------------------------------------------------------
Finter  = spm_figure('GetWin','Interactive');
header  = get(Finter,'Name');
set(Finter,'Name','Dynamic Causal Modelling')
WS      = spm('WinScale');

% prompt for contrast if necessary
%--------------------------------------------------------------------------
if (nargin < 2) || isempty(D)
    str = 'contrast for';
    D   = spm_input(str,1,'b',{'A','B','C'});
end

% set object sizes
%--------------------------------------------------------------------------
dx    = 35;
wx    = 30;
wy    = 20;
uic   = uicontrol(Finter,'String','done','Position',[300 50 060 020].*WS,...
        'Callback', 'uiresume(gcbf)');

% Define left edge dialogue
% itext_left = 080;
% inum_left  = 180;
%--------------------------------------------------------------------------
BackgroundColor = get(get(uic,'Parent'),'Color');
itext_left  = 030;
inum_left   = 80;
name_length = 8; % number of characters in short names
for i=1:n
    if length(region_names{i}) > name_length
        short_name(i).str = region_names{i}(1:name_length);
    else
        short_name(i).str = region_names{i};
    end
end
text_top = 336;

switch D

    % Get contrast weights
    %======================================================================
    case 'A' % intrinsic

        str   = sprintf('%s A: ',title_text);
        spm_input(str,1,'d');

        % Warn if no matrix elements enabled
        any_enabled = any(enabled.A(:));
        if any_enabled
            h0 = [];
        else
            h0 = add_none_enabled_warning();
        end

        % Print names and numbers of regions
        %------------------------------------------------------------------
        for i = 1:n
            str    = sprintf('%s   %i',short_name(i).str,i);

            h1(i) = add_text([itext_left text_top-dx*i 080 020].*WS,...
                        str, any_enabled, 'right');

            h2(i) = add_text([inum_left+dx*i text_top 020 020].*WS,...
                        sprintf('%i',i), any_enabled);
        end

        % Set contrast values and display
        %------------------------------------------------------------------
        for i = 1:n
            for j = 1:n
                cc=ceil([inum_left+dx*j text_top+4-dx*i wx wy].*WS);
                h3(i,j) = add_matrix_input(cc, num2str(A(i,j)), enabled.A(i,j)==1);
            end
        end
        drawnow

        % wait for 'done'
        %------------------------------------------------------------------
        uiwait(Finter);
        for i = 1:n
            for j = 1:n
                A(i,j) = get_matrix_input(h3(i,j));
            end
        end

        % clean up
        %------------------------------------------------------------------
        delete([h0; h1(:); h2(:); h3(:)]);
        spm_input(' ',1,'d');

    case 'B' % modulatory
        %------------------------------------------------------------------
        for k = 1:m,
            str   = sprintf(...
                '%s B: effects of input  %-12s',...
                title_text, input_names{k});
            spm_input(str,1,'d')

            % Warn if no matrix elements enabled
            any_enabled = any(any(enabled.B(:,:,k)));
            if any_enabled
                h0 = [];
            else
                h0 = add_none_enabled_warning();
            end

            % Print names and numbers of regions
            %--------------------------------------------------------------
            for i = 1:n
                str    = sprintf('%s   %i',short_name(i).str,i);

                h1(i) = add_text([itext_left text_top-dx*i 080 020].*WS,...
                        str, any_enabled);

                h2(i) = add_text([inum_left+dx*i text_top 020 020].*WS,...
                        sprintf('%i',i), any_enabled);
            end

            % Set contrast values and display
            %--------------------------------------------------------------
            for i = 1:n
                for j = 1:n
                    cc=ceil([inum_left+dx*j text_top+4-dx*i wx wy].*WS);
                    h3(i,j) = add_matrix_input(cc, num2str(B(i,j,k)), enabled.B(i,j,k)==1);
                end
            end
            drawnow

            % wait for 'done'
            %--------------------------------------------------------------
            set(gcf,'CurrentObject',h3(1));
            uiwait(Finter);
            for i = 1:n
                for j = 1:n
                    B(i,j,k) = get_matrix_input(h3(i,j));
                end
            end

            % clean up
            %--------------------------------------------------------------
            delete([h0; h1(:); h2(:); h3(:)]);
            spm_input(' ',1,'d');

        end

    case 'C' % input
        %------------------------------------------------------------------
        for k = 1:m
            str   = sprintf(...
                '%s C: Effects of input %-12s',...
                title_text, input_names{k});
            spm_input(str,1,'d');

            % Warn if no matrix elements enabled
            any_enabled = any(enabled.C(:,k));
            if any_enabled
                h0 = [];
            else
                h0 = add_none_enabled_warning();
            end

            % Print names and numbers of regions
            %--------------------------------------------------------------
            for i = 1:n
                str    = sprintf('%s   %i',short_name(i).str,i);

                h1(i) = add_text([itext_left text_top-dx*i 080 020].*WS,...
                        str, any_enabled);

                cc = [inum_left+dx text_top+4-dx*i wx wy].*WS;
                h2(i) = add_matrix_input(cc, num2str(C(i,k)), enabled.C(i,k)==1);

            end
            drawnow

            % wait for 'done'
            %--------------------------------------------------------------
            set(gcf,'CurrentObject',h2(1))
            uiwait(Finter);

            for i = 1:n
                C(i,k)   = get_matrix_input(h2(i));
            end

            % clean up
            %--------------------------------------------------------------
            delete([h0; h1(:); h2(:)]);
            spm_input(' ',1,'d');
        end
    otherwise
        disp('Error in spm_dcm_connectivity_ui: matrix must be A, B or C');
        close
        return
end

% Bundle outputs
%--------------------------------------------------------------------------
con = struct('A',A,'B',B,'C',C);

delete(uic);

%--------------------------------------------------------------------------
function h = add_matrix_input(position, default_txt, enabled)
    % Adds an editor for a single element of a matrix
    if enabled == 1
        isvisible = 'On';
    else
        isvisible = 'Off';
    end

    h = uicontrol(Finter,...
        'Position',position,...
        'BackgroundColor',BackgroundColor,...
        'Style','edit',...
        'Visible',isvisible,...
        'String',default_txt);
end
%--------------------------------------------------------------------------
function num = get_matrix_input(h)
    % Reads a number from a single element of a matrix
    entry=get(h,'string');
    entry=strtrim(entry);

    if isempty(entry)
        entry = '0';
    end

    try
        num = str2double(entry);
    catch
        error('Could not read matrix input');
    end
end
%--------------------------------------------------------------------------
function h = add_text(position, default_txt, visible, horizontalAlignment)
    % Adds a text label to the GUI
    if nargin < 3
        visible = 1;
    end
    if nargin < 4
        horizontalAlignment = 'left';
    end

    if visible
        visible = 'On';
    else
        visible = 'Off';
    end

    h  = uicontrol(Finter,'String',default_txt,...
        'Style','text',...
        'HorizontalAlignment',horizontalAlignment,...
        'BackgroundColor',BackgroundColor,...
        'Position',position,...
        'Visible',visible);
end
%--------------------------------------------------------------------------
function h = add_none_enabled_warning()
    % Adds a warning for when there are no enabled matrix entries
    position = [itext_left text_top-dx*1 250 020].*WS;
    nonestring = 'Nothing to enter - press done to continue';
    h = add_text(position, nonestring);
end

end

author: spm


text:
classdef dme_branch < mp.dm_element
%MP.DME_BRANCH  MATPOWER data model class for branch data

%   MATPOWER
%   Copyright (c) 2020-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

    properties
        fbus    %% bus index vector for "from" port (port 1) (all branches)
        tbus    %% bus index vector for "to" port (port 2) (all branches)
        r       %% series resistance (p.u.) for branches that are on
        x       %% series reactance (p.u.) for branches that are on
        g_fr    %% shunt conductance (p.u.) at "from" end for branches that are on
        g_to    %% shunt conductance (p.u.) at "to" end for branches that are on
        b_fr    %% shunt susceptance (p.u.) at "from" end for branches that are on
        b_to    %% shunt susceptance (p.u.) at "to" end for branches that are on
        tm      %% transformer off-nominal turns ratio for branches that are on
        ta      %% xformer phase-shift angle (radians) for branches that are on
        rate_a  %% long term flow limit (p.u.) for branches that are on
    end     %% properties

    methods
        function name = name(obj)
            name = 'branch';
        end

        function label = label(obj)
            label = 'Branch';
        end

        function label = labels(obj)
            label = 'Branches';
        end

        function name = cxn_type(obj)
            name = 'bus';
        end

        function name = cxn_idx_prop(obj)
            name = {'fbus', 'tbus'};
        end

        function names = main_table_var_names(obj)
            names = horzcat( main_table_var_names@mp.dm_element(obj), ...
                {'bus_fr', 'bus_to', 'r', 'x', 'g_fr', 'b_fr', ...
                'g_to', 'b_to', 'sm_ub_a', 'sm_ub_b', 'sm_ub_c', ...
                'cm_ub_a', 'cm_ub_b', 'cm_ub_c', 'vad_lb', 'vad_ub', ...
                'tm', 'ta', ... %% remove these when we separate out xformers
                'pl_fr', 'ql_fr', 'pl_to', 'ql_to'} );
        end

        function vars = export_vars(obj)
            vars = {'pl_fr', 'ql_fr', 'pl_to', 'ql_to'};
        end

        function obj = initialize(obj, dm)
            initialize@mp.dm_element(obj, dm);  %% call parent

            %% get bus mapping info
            b2i = dm.elements.bus.ID2i;     %% bus num to idx mapping

            %% set bus index vectors for port connectivity
            obj.fbus = b2i(obj.tab.bus_fr);
            obj.tbus = b2i(obj.tab.bus_to);
        end

        function obj = update_status(obj, dm)
            %% get bus status info
            bs = dm.elements.bus.tab.status;    %% bus status

            %% update status of branches connected to isolated/offline buses
            obj.tab.status = obj.tab.status & bs(obj.fbus) & ...
                                              bs(obj.tbus);

            %% call parent to fill in on/off
            update_status@mp.dm_element(obj, dm);
        end

        function obj = build_params(obj, dm)
            obj.r  = obj.tab.r(obj.on);
            obj.x  = obj.tab.x(obj.on);
            obj.g_fr  = obj.tab.g_fr(obj.on);
            obj.b_fr  = obj.tab.b_fr(obj.on);
            obj.g_to  = obj.tab.g_to(obj.on);
            obj.b_to  = obj.tab.b_to(obj.on);
            obj.tm = obj.tab.tm(obj.on);
            obj.ta = obj.tab.ta(obj.on) * pi/180;
            obj.rate_a = obj.tab.sm_ub_a(obj.on) / dm.base_mva;
        end

        function obj = pp_data_cnt(obj, dm, rows, out_e, mpopt, fd, pp_args)
            %% call parent
            pp_data_cnt@mp.dm_element(obj, dm, rows, out_e, mpopt, fd, pp_args);

            num_xf = length(find(obj.tab.tm));
            num_xf_on = length(find(obj.tab.tm(obj.on)));
            num_ln = obj.nr - num_xf;
            num_ln_on = obj.n - num_xf_on;
            ln = sprintf('%d', num_ln);
            if num_ln == num_ln_on
                ln_on = ln;
                ln_off = '-';
            else
                ln_on = sprintf('%d', num_ln_on);;
                ln_off = sprintf('%d', num_ln - num_ln_on);
            end
            xf = sprintf('%d', num_xf);
            if num_xf == num_xf_on
                xf_on = xf;
                xf_off = '-';
            else
                xf_on = sprintf('%d', num_xf_on);;
                xf_off = sprintf('%d', num_xf - num_xf_on);
            end

            %% print line, transformer counts
            fprintf(fd, '  %-20s%7s %7s %7s\n', '  Lines', ln_on, ln_off, ln);
            fprintf(fd, '  %-20s%7s %7s %7s\n', '  Transformers', xf_on, xf_off, xf);
        end

        function TorF = pp_have_section_sum(obj, mpopt, pp_args)
            TorF = true;
        end

        function obj = pp_data_sum(obj, dm, rows, out_e, mpopt, fd, pp_args)
            %% call parent
            pp_data_sum@mp.dm_element(obj, dm, rows, out_e, mpopt, fd, pp_args);

            %% print branch summary
            fprintf(fd, '  %-29s  %12.2f MW', 'Total branch losses', ...
                sum(obj.tab.pl_fr(obj.on)) + sum(obj.tab.pl_to(obj.on)) );
            if mpopt.model(1) ~= 'D'    %% AC model
                fprintf(fd, ' %12.2f MVAr', ...
                    sum(obj.tab.ql_fr(obj.on)) + sum(obj.tab.ql_to(obj.on)) );
            end
            fprintf(fd, '\n');
        end

        function h = pp_get_headers_det(obj, dm, out_e, mpopt, pp_args)
            h = [ pp_get_headers_det@mp.dm_element(obj, dm, out_e, mpopt, pp_args) ...
                {   ' Branch     From       To             From Bus Injection   To Bus Injection', ...
                    '   ID      Bus ID    Bus ID   Status   P (MW)   Q (MVAr)   P (MW)   Q (MVAr)', ...
                    '--------  --------  --------  ------  --------  --------  --------  --------' } ];
            %%       1234567 123456789 123456789 -----1 1234567.90 123456.89 123456.89 123456.89
        end

        function TorF = pp_have_section_det(obj, mpopt, pp_args)
            TorF = true;
        end

        function str = pp_data_row_det(obj, dm, k, out_e, mpopt, fd, pp_args)
            str = sprintf('%7d %9d %9d %6d %10.2f %9.2f %9.2f %9.2f', ...
                obj.tab.uid(k), obj.tab.bus_fr(k), obj.tab.bus_to(k), ...
                obj.tab.status(k), ...
                obj.tab.pl_fr(k), obj.tab.ql_fr(k), ...
                obj.tab.pl_to(k), obj.tab.ql_to(k) );
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% Generate score matrix for the 2D-2D cue
% Author: Junaid Ahmed Ansari @ Robotics Research Center, IIIT Hyderabad
% Email : junaid.ansari@research.iiit.ac.in
function [scoreMatrix_2D2D] = generate2D2DScoreMatrix(featureQ, featureT)

    scoreMatrix_2D2D = zeros(size(featureQ,1), size(featureT,1));

    % 2D 2D score
    for i = 1:size(featureQ,1)
        for j = 1:size(featureT,1)

            fQ = featureQ(i,2:4097) - mean(featureQ(i,2:4097));
            fT = featureT(j,2:4097) - mean(featureT(j,2:4097));
            fQ = fQ/norm(fQ);
            fT = fT/norm(fT);

            score = sum(fQ .* fT);
            scoreMatrix_2D2D(i,j) = score;

        end
    end

end

repo: MOTBeyondPixels


text:
function M = elliptopefactory(n, k)
% Manifold of n-by-n psd matrices of rank k with unit diagonal elements.
%
% function M = elliptopefactory(n, k)
%
% A point X on the manifold is parameterized as YY^T where Y is a matrix of
% size nxk. As such, X is symmetric, positive semidefinite. We restrict to
% full-rank Y's, such that X has rank exactly k. The point X is numerically
% represented by Y (this is more efficient than working with X, which may
% be big). Tangent vectors are represented as matrices of the same size as
% Y, call them Ydot, so that Xdot = Y Ydot' + Ydot Y and diag(Xdot) == 0.
% The metric is the canonical Euclidean metric on Y.
%
% The diagonal constraints on X (X(i, i) == 1 for all i) translate to
% unit-norm constraints on the rows of Y: norm(Y(i, :)) == 1 for all i.
% The set of such Y's forms the oblique manifold. But because for any
% orthogonal Q of size k it holds that (YQ)(YQ)' = YY', we "group" all
% matrices of the form YQ in an equivalence class. The set of equivalence
% classes is a Riemannian quotient manifold, implemented here.
%
% Note that this geometry formally breaks down at rank-deficient Y's.
% This does not appear to be a major issue in practice when optimization
% algorithms converge to rank-deficient Y's, but convergence theorems no
% longer hold. As an alternative, you may use the oblique manifold (it has
% larger dimension, but does not break down at rank drop.)
%
% The geometry is taken from the 2010 paper:
% M. Journee, P.-A. Absil, F. Bach and R. Sepulchre,
% "Low-Rank Optimization on the Cone of Positive Semidefinite Matrices".
% Paper link: http://www.di.ens.fr/~fbach/journee2010_sdp.pdf
%
%
% Please cite the Manopt paper as well as the research paper:
%     @Article{journee2010low,
%       Title   = {Low-rank optimization on the cone of positive semidefinite matrices},
%       Author  = {Journ{\'e}e, M. and Bach, F. and Absil, P.-A. and Sepulchre, R.},
%       Journal = {SIAM Journal on Optimization},
%       Year    = {2010},
%       Number  = {5},
%       Pages   = {2327--2351},
%       Volume  = {20},
%       Doi     = {10.1137/080731359}
%     }
%
%
% See also: obliquefactory symfixedrankYYfactory spectrahedronfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, July 12, 2013.
% Contributors:
% Change log:
%   July 18, 2013 (NB):
%       Fixed projection operator for rank-deficient Y'Y.
%
%   Aug.  8, 2013 (NB):
%       No longer using nested functions, to aim at Octave compatibility.
%       Sign error in right hand side of the call to minres corrected.
%
%   June 24, 2014 (NB):
%       Used code snippets from obliquefactory to speed up projection,
%       retraction, egrad2rgrad and rand: the code now uses bsxfun for this.
%
%   Apr. 3, 2015 (NB):
%       Replaced trace(A'*B) by A(:)'*B(:) : equivalent but faster.
%
%   Apr. 17, 2018 (NB):
%       Removed dependency on lyap entirely.
%
%   Sep.  6, 2018 (NB):
%       Removed M.exp() as it was not implemented.

% TODO: modify normalize_rows and project_rows to work without transposes.
% TODO: enhance ehess2rhess to also use bsxfun.

    if k < 2
        warning('manopt:elliptopefactory:lowk', ...
                'k should be an integer >= 2. At k = 1, the set is discrete.');
    end

    M.name = @() sprintf('YY'' quotient manifold of %dx%d psd matrices of rank %d with diagonal elements being 1', n, k);

    M.dim = @() n*(k-1) - k*(k-1)/2; % Extra -1 is because of the diagonal constraint that

    % Euclidean metric on the total space
    M.inner = @(Y, eta, zeta) eta(:)'*zeta(:);

    M.norm = @(Y, eta) sqrt(M.inner(Y, eta, eta));

    M.dist = @(Y, Z) error('elliptopefactory.dist not implemented yet.');

    M.typicaldist = @() 10*k;

    M.proj = @projection;

    M.tangent = M.proj;
    M.tangent2ambient = @(Y, eta) eta;

    M.retr = @retraction;

    M.egrad2rgrad = @egrad2rgrad;

    M.ehess2rhess = @ehess2rhess;

    % Notice that the hash of two equivalent points will be different...
    M.hash = @(Y) ['z' hashmd5(Y(:))];

    M.rand = @() random(n, k);

    M.randvec = @randomvec;

    M.lincomb = @matrixlincomb;

    M.zerovec = @(Y) zeros(n, k);

    M.transp = @(Y1, Y2, d) projection(Y2, d);

    M.vec = @(Y, u_mat) u_mat(:);
    M.mat = @(Y, u_vec) reshape(u_vec, [n, k]);
    M.vecmatareisometries = @() true;

end

% Given a matrix X, returns the same matrix but with each column scaled so
% that they have unit 2-norm.
% See obliquefactory.
function X = normalize_rows(X)
    X = X';
    norms = sqrt(sum(X.^2, 1));
    X = bsxfun(@times, X, 1./norms);
    X = X';
end

% Orthogonal projection of each row of H to the tangent space at the
% corresponding row of X, seen as a point on a sphere.
% See obliquefactory.
function PXH = project_rows(X, H)
    X = X';
    H = H';
    % Compute the inner product between each vector H(:, i) with its root
    % point X(:, i), that is, X(:, i).' * H(:, i). Returns a row vector.
    inners = sum(X.*H, 1);
    % Subtract from H the components of the H(:, i)'s that are parallel to
    % the root points X(:, i).
    PXH = H - bsxfun(@times, X, inners);
    PXH = PXH';
end

% Projection onto the tangent space, i.e., on the tangent space of
% ||Y(i, :)|| = 1
function etaproj = projection(Y, eta)
    eta = project_rows(Y, eta);

    % Projection onto the horizontal space
    YtY = Y'*Y;
    SS = YtY;
    AS = Y'*eta - eta'*Y;
    Omega = lyapunov_symmetric(SS, AS);

    % It does not seem necessary to enforce skew-symmetry numerically.
    % Omega = (Omega-Omega')/2;

    etaproj = eta - Y*Omega;
end

% Retraction
function Ynew = retraction(Y, eta, t)
    if nargin < 3
        t = 1.0;
    end
    Ynew = Y + t*eta;
    Ynew = normalize_rows(Ynew);
end

% Euclidean gradient to Riemannian gradient conversion.
% We only need the ambient space projection: the remainder of the
% projection function is not necessary because the Euclidean gradient must
% already be orthogonal to the vertical space.
function rgrad = egrad2rgrad(Y, egrad)
    rgrad = project_rows(Y, egrad);
end

% Euclidean Hessian to Riemannian Hessian conversion.
% TODO: speed this function up using bsxfun.
function Hess = ehess2rhess(Y, egrad, ehess, eta)
    k = size(Y, 2);

    % Directional derivative of the Riemannian gradient
    scaling_grad = sum((egrad.*Y), 2); % column vector of size n
    scaling_grad_repeat = scaling_grad*ones(1, k);

    Hess = ehess - scaling_grad_repeat.*eta;

    scaling_hess = sum((eta.*egrad) + (Y.*ehess), 2);
    scaling_hess_repeat = scaling_hess*ones(1, k);
    % directional derivative of scaling_grad_repeat
    Hess = Hess - scaling_hess_repeat.*Y;

    % Project on the horizontal space
    Hess = projection(Y, Hess);
end

% Random point generation on the manifold
function Y = random(n, k)
    Y = randn(n, k);
    Y = normalize_rows(Y);
end

% Random vector generation at Y
function eta = randomvec(Y)
    eta = randn(size(Y));
    eta = projection(Y, eta);
    nrm = norm(eta, 'fro');
    eta = eta / nrm;
end

author: NicolasBoumal

repo: manopt


text:
function outsig = signal_from_spectrum(amplitude,phase,f,dim,conf)
%SIGNAL_FROM_SPECTRUM time signal from single-sided spectrum
%
%   Usage: outsig = signal_from_spectrum(amplitude,phase,f,[dim],conf)
%
%   Input parameters:
%       amplitude   - the single-sided amplitude spectrum
%       phase       - the single-sided phase spectrum / rad
%       f           - the corresponding frequency vector
%       dim         - dimension along which the ifft is performed
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       outsig      - one channel audio (time) signal
%
%   SIGNAL_FROM_SPECTRUM(amplitude,phase,f,conf) generates the time signal from
%   single-sided amplitude and phase spectra using ifft. It is the counterpart
%   of SPECTRUM_FROM_SIGNAL and not interchangeable with calling ifft.
%
%   See also: spectrum_from_signal, ifft

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking input arguments ========================================
nargmin = 4;
nargmax = 5;
narginchk(nargmin,nargmax);
if nargin == nargmin
   conf = dim;
   dim = find(size(amplitude)~=1, 1);  % find first non-singleton dimension
else
   isargpositivescalar(dim);
end
isargstruct(conf);

Ndims = ndims(amplitude);
dim = min(dim,Ndims);
amplitude = permute(amplitude, [dim:Ndims, 1:dim-1]);  % move dim to first dimension
phase = permute(phase, [dim:Ndims, 1:dim-1]);  % move dim to first dimension
s = size(amplitude);
Nx = s(1);
amplitude = reshape(amplitude, Nx, []);  % squeeze all other dimensions
phase = reshape(phase, Nx, []);  % squeeze all other dimensions

%% ===== Configuration ===================================================
fs = conf.fs;

%% ===== Regenerating wave form from spectrum ============================
% Provided number of frequency bins
bins = length(f);
if bins ~= Nx
   error('%s: size of spectrum does not match length of frequency vector',...
       upper(mfilename));
end

if f(end) == fs/2  % -> even time signal length
    % Length of the signal to generate
    samples = 2 * (bins-1);
    % Rescaling (see spectrum_from_signal())
    amplitude = [amplitude(1,:); amplitude(2:end-1,:)/2; amplitude(end,:)] ...
        * samples;
    % Mirror the amplitude spectrum ( 2*pi periodic [0, fs[ )
    amplitude = [amplitude; amplitude(end-1:-1:2,:)];
    % Mirror the phase spectrum and build the inverse (complex conjugate)
    phase = [phase; -1 * phase(end-1:-1:2,:)];

else  % -> odd time signal length
    % Length of the signal to generate
    samples = 2*bins - 1;
    % Rescaling (see signal_from_spectrum)
    amplitude = [amplitude(1,:); amplitude(2:end,:)/2] * samples;
    % Mirror the amplitude spectrum ( 2*pi periodic [0, fs-bin] )
    amplitude = [amplitude; amplitude(end:-1:2,:)];
    % Mirror the phase spectrum and build the inverse (complex conjugate)
    phase = [phase; -1*phase(end:-1:2,:)];
end

% Convert to complex spectrum
compspec = amplitude .* exp(1i*phase);

% Build the inverse fft and assume spectrum is conjugate symmetric
outsig = real(ifft(compspec));

%% ===== Output ==========================================================
% undo reshape and permute
outsig = reshape(outsig, [samples, s(2:end)]);
outsig = permute(outsig, [Ndims-dim+2:Ndims, 1:Ndims-dim+1]);

author: sfstoolbox

repo: sfs-matlab


text:
function file_name_sequence ( )

%*****************************************************************************80
%
%% FILE_NAME_SEQUENCE demonstrates ways of generating a sequence of filenames.
%
%  Discussion:
%
%    There are situations such as animations or parallel processing in which
%    it is necessary to generate a sequence of file names which include
%    an embedded index that increases.  A simple example might be
%
%      'fred0.txt', 'fred1.txt', 'fred2.txt'
%
%    A side issue arises when the number of files is large enough that the
%    number of digits in the index will vary.  Thus, if we are going to have
%    15 files, do we want to number them as
%
%      'fred00.txt' through 'fred14.txt'
%
%    which means, for one thing, that they will alphabetize properly, or
%    will we be satisfied with
%
%      'fred0.txt' through 'fred14.txt' ?
%
%    Schemes for generating such a sequence in MATLAB can involve the
%    NUM2STR function, the SPRINTF function, or a more elaborate function
%    called FILE_NAME_INC which searches a string for embedded numeric
%    data and increments it.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 March 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FILE_NAME_SEQUENCE:\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Demonstrate ways of generating a numeric sequence of file names.\n' );

  test01 ( 'frad', '.txt', 0, 12 );
  test02 ( 'fred', '.txt', 0, 12 );
  test03 ( 'frid', '.txt', 99, 105 );
  test04 ( 'frodo_01345_lives.txt', 10 );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FILE_NAME_SEQUENCE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

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

  return
end
function test01 ( prefix, suffix, first, last )

%*****************************************************************************80
%
%% TEST01 uses NUM2STR and concatenation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 March 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01:\n' );
  fprintf ( 1, '  FILENAME(I) = PREFIX + NUM2STR(I) + SUFFIX\n' );
  fprintf ( 1, '  PREFIX = "%s"\n', prefix );
  fprintf ( 1, '  SUFFIX = "%s"\n', suffix );
  fprintf ( 1, '  %d <= I <= %d\n', first, last );
  fprintf ( 1, '  Numbers do NOT include leading zeros.\n' );
  fprintf ( 1, '\n' );

  for i = first : last
    filename = strcat ( prefix, num2str ( i ), suffix );
    fprintf ( 1, '  %4d:  "%s"\n', i, filename );
  end

  return
end
function test02 ( prefix, suffix, first, last )

%*****************************************************************************80
%
%% TEST02 uses SPRINTF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 March 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02:\n' );
  fprintf ( 1, '  FILENAME(I) = SPRINTF ( PREFIX, I, SUFFIX )\n' );
  fprintf ( 1, '  PREFIX = "%s"\n', prefix );
  fprintf ( 1, '  SUFFIX = "%s"\n', suffix );
  fprintf ( 1, '  %d <= I <= %d\n', first, last );
  fprintf ( 1, '  Numbers do NOT include leading zeros.\n' );
  fprintf ( 1, '\n' );

  for i = first : last
    filename = sprintf ( '%s%d%s', prefix, i, suffix );
    fprintf ( 1, '  %4d:  "%s"\n', i, filename );
  end

  return
end
function test03 ( prefix, suffix, first, last )

%*****************************************************************************80
%
%% TEST03 uses SPRINTF and leading 0's.
%
%  Discussion:
%
%    The format in SPRINTF must include a leading "0" for the integer,
%    plus a desired width.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 March 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03:\n' );
  fprintf ( 1, '  FILENAME(I) = SPRINTF ( PREFIX, I, SUFFIX )\n' );
  fprintf ( 1, '  PREFIX = "%s"\n', prefix );
  fprintf ( 1, '  SUFFIX = "%s"\n', suffix );
  fprintf ( 1, '  %d <= I <= %d\n', first, last );
  fprintf ( 1, '  Numbers DO include leading zeros.\n' );
  fprintf ( 1, '\n' );

  for i = first : last
    filename = sprintf ( '%s%04d%s', prefix, i, suffix );
    fprintf ( 1, '  %4d:  "%s"\n', i, filename );
  end

  return
end
function test04 ( filename, filename_num )

%*****************************************************************************80
%
%% TEST04 uses FILENAME_INC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 March 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04:\n' );
  fprintf ( 1, '  FILENAME(I+1) = FILENAME_INC ( FILENAME(I) )\n' );
  fprintf ( 1, '  First FILENAME = "%s"\n', filename );
  fprintf ( 1, '  Number of filenames = %d\n', filename_num );
  fprintf ( 1, '  Numbers may include leading zeros.\n' );
  fprintf ( 1, '\n' );

  for i = 1 : filename_num
    fprintf ( 1, '  %4d:  "%s"\n', i, filename );
    filename = filename_inc ( filename );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [eeg] = read_eep_cnt(fn);

% READ_EEP_CNT reads continuous EEG data from an EEProbe *.cnt file
% and returns a structure containing the header and data information.
%
% eeg = read_eep_cnt(filename, sample1, sample2)
%
% where sample1 and sample2 are the begin and end sample of the data
% to be read.
%
% eeg.label    ... labels of EEG channels
% eeg.rate     ... sampling rate
% eeg.npnt     ... number of sample in data segment
% eeg.nchan    ... number of channels
% eeg.nsample
% eeg.time     ... array [1 x npnt] of time points (ms)
% eeg.data     ... array [nchan x npnt] containing eeg data (uV)
% eeg.trigger  ... array [n, 4] containing trigger information
%
% Author: Robert Oostenveld, Aalborg University, Denmark, 11 March 2003
%
% See also READ_EEP_TRG, READ_EEP_REJ, READ_EEP_AVR
%

% Copyright (C) 2002, Robert Oostenveld
%                     Aalborg University, Denmark
%                     http://www.smi.auc.dk/~roberto/
%
% 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

% $Log: not supported by cvs2svn $
% Revision 1.1  2004/11/19 14:55:38  jwiskerke
% Added files for use with matlab
%
% Revision 1.1  2004/11/19 14:39:29  jwiskerke
% Initial input into cvs. These files should work.
%
% Revision 1.2  2003/10/24 13:34:41  Maarten-Jan Hoeve
% Added GNU Licence and updated revision history
%
% Revision 1.1.1.2  2003/10/17 09:55:20  mvelde
% updated: consistent copyrights, arguments/data labels, fixed some typos
%
% Revision 1.1.1.1  2003/03/11 15:24:51  roberto
% updated help and copyrights
%
% ANT Software BV, The Netherlands, www.ant-neuro.com / info@ant-neuro.com
%

error('could not locate mex file');


text:
% For transforming the sparse adjacency matrices in 'DD' and 'NCI1' to dense
% in order to be read by torch.matio
% Note torch.matio does not support dense .mat

dataset = strvcat('DD', 'NCI1');
for ith_data = 1: size(dataset, 1)
    clear save_struct
    data = strcat(dataset(ith_data, :))
    dataname = strcat('../data/raw_data/', data, '.mat');
    load(dataname);
    X = eval(data);
    for i = 1: size(X, 2)
        am = full(X(i).am);
        l = X(i).nl.values;
        X(i).am = am;
    end
    save_struct.(data) = X;
    save(dataname, '-struct', 'save_struct', '-append');
end

author: muhanzhang

repo: DGCNN

save_path: github-repos/MATLAB/muhanzhang-DGCNN


text:
function f = initialize_variables(N, M, V, min_range, max_range)

%% function f = initialize_variables(N, M, V, min_range, max_range)
% This function initializes the chromosomes. Each chromosome has the
% following at this stage
%       * set of decision variables
%       * objective function values
%
% where,
% N - Population size
% M - Number of objective functions
% V - Number of decision variables
% min_range - A vector of decimal values which indicate the minimum value
% for each decision variable.
% max_range - Vector of maximum possible values for decision variables.

%  Copyright (c) 2009, Aravind Seshadri
%  All rights reserved.
%

min = min_range;
max = max_range;

% K is the total number of array elements. For ease of computation decision
% variables and objective functions are concatenated to form a single
% array. For crossover and mutation only the decision variables are used
% while for selection, only the objective variable are utilized.

K = M + V;

%% Initialize each chromosome
% For each chromosome perform the following (N is the population size)
f = zeros(N,K); % modified by zzb
for i = 1 : N
    % Initialize the decision variables based on the minimum and maximum
    % possible values. V is the number of decision variable. A random
    % number is picked between the minimum and maximum possible values for
    % the each decision variable.
    for j = 1 : V
        f(i,j) = min(j) + (max(j) - min(j))*rand(1);
    end
    % For ease of computation and handling data the chromosome also has the
    % vlaue of the objective function concatenated at the end. The elements
    % V + 1 to K has the objective function valued.
    % The function evaluate_objective takes one chromosome at a time,
    % infact only the decision variables are passed to the function along
    % with information about the number of objective functions which are
    % processed and returns the value for the objective functions. These
    % values are now stored at the end of the chromosome itself.
    f(i,V + 1: K) = evaluate_objective(f(i,1:V), M);
end

author: sfvsfv

repo: Mathematical-modeling


text:
function nc = wavepaste(type,c,s,n,x)
%WAVEPASTE Puts coefficients in a wavelet decomposition structure.
%   NC = WAVEPASTE(TYPE,C,S,N,X) returns the new decomposition structure
%   after pasting X into it based on TYPE and N.
%
%   INPUTS:
%     TYPE      Coefficient category
%     -------------------------------------
%     'a'       Approximation coefficients
%     'h'       Horizontal details
%     'v'       Vertical details
%     'd'       Diagonal details
%
%     [C, S] is a wavelet data structure.
%     N specifies a decomposition level (Ignored if TYPE = 'a').
%     X is a 2- or 3-D approximation or detail coefficient matrix whose
%       dimensions are appropriate for decomposition level N.
%
%   See also WAVEWORK, WAVECUT, and WAVECOPY.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

narginchk(5,5)
nc = wavework('paste',type,c,s,n,x);

author: dipum

repo: dipum-toolbox

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


text:
% EEG_GLOBAL - declare global EEGLAB variables. These variables are
%                used only as global by the main function EEGLAB,
%                the function POP_REJMENU and the history EEGH function.
%
% Global variables:
%  EEG        - structure containing the current dataset
%  ALLEEG     - array of structures containing all the loaded datasets
%  CURRENTSET - index of the current dataset in the ALLEEG array
%  LASTCOM    - most recent command run by EEGLAB
%  ALLCOM     - cell array containing the EEGLAB session command history
%               See the history (h) function for more details.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: EEGLAB

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% global variables
% ----------------
globalvars = who('global');
if ~isempty(strmatch('ALLCOM', globalvars, 'exact')) || exist('ALLCOM') ~= 1
    tmpEEG = whos('EEG');
    % avoid overwritting variable in the current space
    % in case it is not global
    if ~isempty(tmpEEG) && ~tmpEEG.global
        TMPEEG = EEG;
        clear EEG;
        global EEG;
        EEG = TMPEEG;
    else
        global EEG;		% current dataset
    end
    global ALLEEG;		% all datasets
    global CURRENTSET;	% current set index

    %global W_MAIN;		% main window
    %global H_MAIN;		% main window
    %global EEGMENU;		% main menu
    global ALLCOM;		% all commands (history)
    global LASTCOM;		% last command
    global STUDY;
    global CURRENTSTUDY;
    global PLUGINLIST;

    if exist('eegplugin_erplab.m')
        global ALLERP; % Javier Lopez-Calderon for ERPLAB
    end
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function etaButtonCallback(blockHandle)
% Callback function for ETAButton
% Allows user to specify what .mat file to use for elevationImport wave type

    [filename,filepath] = uigetfile('.mat');               % Return filename, filepath of user-chosen file
    % Don't set value if no file is chosen, or prompt canceled.
    if ~isequal(filename,0) && ~isequal(filepath,0)
        mask = Simulink.Mask.get(blockHandle);
        fileParam = mask.getParameter('elevationFile');
        fileParam.Value = [filepath,filename];             % Update elevationFile with new filename
    end
end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function [A_hat E_hat iter svp elapsed] = ialm_rpca(D, lambda, tol, maxIter, blk)

% Oct 2009
% This matlab code implements the inexact augmented Lagrange multiplier
% method for Robust PCA.
%
% D - m x n matrix of observations/data (required input)
%
% lambda - weight on sparse error term in the cost function
%
% tol - tolerance for stopping criterion.
%     - DEFAULT 1e-7 if omitted or -1.
%
% maxIter - maximum number of iterations
%         - DEFAULT 1000, if omitted or -1.
% blk - indicate whether to use BLWS
%
% Initialize A,E,Y,u
% while ~converged
%   minimize (inexactly, update A and E only once)
%     L(A,E,Y,u) = |A|_* + lambda * |E|_1 + <Y,D-A-E> + mu/2 * |D-A-E|_F^2;
%   Y = Y + \mu * (D - A - E);
%   \mu = \rho * \mu;
% end

%addpath PROPACK;

%elapsed = tic;

[m n] = size(D);

if nargin < 2
    lambda = 1 / sqrt(max(m,n));
elseif lambda == -1
    lambda = 1 / sqrt(max(m,n));
end

if nargin < 3
    tol = 1e-7;
elseif tol == -1
    tol = 1e-7;
end

if nargin < 4
    maxIter = 1000;
elseif maxIter == -1
    maxIter = 1000;
end

if nargin < 5
    blk = 0;
elseif blk == -1
    blk = 0;
end

% initialize
Y = D;
norm_two = lansvd(Y, 1, 'L');
norm_inf = norm( Y(:), inf) / lambda;
dual_norm = max(norm_two, norm_inf);
Y = Y / dual_norm;

A_hat = zeros( m, n);
E_hat = zeros( m, n);
mu = 1.25/norm_two % this one can be tuned
mu_bar = mu * 1e7
rho = 1.2         % this one can be tuned
d_norm = norm(D, 'fro');

iter = 0;
total_svd = 0;
converged = false;
stopCriterion = 1;
sv = 10;
d = min(m, n);

mark = false;
Block_mark = false ;

while ~converged
    iter = iter + 1;

    temp_T = D - A_hat + (1/mu)*Y;
    E_hat = max(temp_T - lambda/mu, 0);
    E_hat = E_hat+min(temp_T + lambda/mu, 0);

%     if Block_mark==true
%         uv_temp=[U_temp; V_temp]/sqrt(2);
%         [U S V]=BL_SVD(D - E_hat + (1/mu)*Y, U_temp, V_temp, 1);
%     else
          [U S V] = lansvd(D - E_hat + (1/mu)*Y, sv, 'L');
%     end

    if mark==true
        U_temp = U;
        V_temp = V;
        Block_mark=true;
    end

    diagS = diag(S);
    svp = length(find(diagS > 1/mu));

    len = length(diagS);
    ratio = diagS(1:len-1) ./ diagS(2:len);
%     ind = find(ratio > 2);
%     if blk && length(ind) > 0
%         svp = min(svp, ind(1));
%         mark = true;
%     end

    if svp < sv
        sv = min(svp + 1, d);
    else
        sv = min(svp + round(0.05*d), d);
    end

    A_hat = U(:, 1:svp) * diag(diagS(1:svp) - 1/mu) * V(:, 1:svp)';

    total_svd = total_svd + 1;

    Z = D - A_hat - E_hat;

    Y = Y + mu*Z;
    mu = min(mu*rho, mu_bar);

    %% stop Criterion
    stopCriterion = norm(Z, 'fro') / d_norm;
    if stopCriterion < tol
        converged = true;
    end

    if mod( total_svd, 10) == 0
        disp(['#svd ' num2str(total_svd) ' r(A) ' num2str(svp)...
            ' |E|_0 ' num2str(length(find(abs(E_hat)>0)))...
            ' stopCriterion ' num2str(stopCriterion)]);
    end

    if ~converged && iter >= maxIter
        disp('Maximum iterations reached') ;
        converged = 1 ;
    end
end

%elapsed = toc(elapsed);

author: andrewssobral

repo: lrslibrary


text:
function f = fracDiff(f, mu, type)
%FRACDIFF  Fractional derivative of a CHEBFUN.
%   FRACINT(F, MU) gives the order MU Riemann-Liouville fractional derivative of
%   a CHEBFUN object F.
%
%   FRACINT(F, MU, 'Caputo') instead uses the Caputo definition of the
%   fractional derivative.
%
%   Currently this only supports the situation where F is smooth (i.e., it has
%   no breakpoints or endpoint singularities) and on a finite domain.

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

% Default to Riemann-Liouville:
if ( nargin < 3  )
    type = 'RL';
end

% Extract the fractional part:
n = ceil(mu);

if ( n == mu )
    f = diff(f, n);
    return
end

if ( numel(f) > 0 )
    f = quasimatrix(f);
    for k = 1:numel(f)
        f(k) = fracDiffcol(f(k), mu, type);
    end
else
    f = fracDiffcol(f, mu, type);
end

end

function f = fracDiffcol(f, mu, type)

% No piecewise support yet:
if ( numel(f(1).funs) > 1 )
    error('CHEBFUN:CHEBFUN:fracDiff:breakpoints', ...
        'FRACDIFF does not currently support piecewise functions.');
end

% Extract the fractional part:
n = ceil(mu);

if ( strcmpi(type, 'Caputo') )
    % Caputo:
    f = diff(f, n);
    f = fracInt(f, n - mu);

else
    % Riemann-Liouville:
    f = fracInt(f, n - mu);
    f = diff(f, n);

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function el = export_structure_set_ROI_sequence(args)
%"export_structure_set_ROI_sequence"
%   Subfunction to handle structure_set_ROI sequences within the
%   structure_set module.  Uses the same layout and principle as the parent
%   function.
%
%   This function takes a CERR structure element and the index of that
%   element within the planC.structures array.
%
%JRA 06/23/06
%NAV 07/19/16 updated to dcm4che3
%   replaced ml2dcm_Element to data2dcmElement
%
%Usage:
%   @export_structure_set_ROI_sequence(args)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Init output element to empty.
el = [];

%Unpack input data.
tag         = args.tag;
structS     = args.data{1};
index       = args.data{2};
template    = args.template;

switch tag
    case 805699618  %3006,0020  ROI Number
        data = uint32(index);
       % el = template.get(tag);
        el = data2dcmElement(data, tag);

    case 805699620  %3006,0024  Referenced Frame of Reference UID
        UID = structS.Frame_Of_Reference_UID;
        el = data2dcmElement(UID, tag);

    case 805699622  %3006,0026  ROI Name
        data = structS.structureName;
        el = data2dcmElement(data, tag);

    case 805699624  %3006,0028  ROI Description
        %Currently unsupported.

    case 805699628  %3006,002C  ROI Volume
        %Currently unsupported.

    case 805699638  %3006,0036  ROI Generation Algorithm
       % el = template.get(tag);
        %el = org.dcm4che3.data.Attributes;
        %el.setString(tag, template.getVR(tag), template.getString(tag));
        data = structS.roiGenerationAlgorithm;
        el = data2dcmElement(data, tag);

    case 805699640  %3006,0038  ROI Generation Description
        data = structS.roiGenerationDescription;
        el = data2dcmElement(data, tag);

    otherwise
        warning(['No methods exist to populate DICOM structure_set module''s structure_set_ROI sequence field: ' dec2hex(tag,8) '.']);
end

author: cerr

repo: CERR


text:
function feature_map = extract_features(image, pos, scales, features, gparams, extract_info)

% Sample image patches at given position and scales. Then extract features
% from these patches.
% Requires that cell size and image sample size is set for each feature.

if ~iscell(features)
    error('Wrong input');
end;

num_features = length(features);
num_scales = length(scales);
num_sizes = length(extract_info.img_sample_sizes);

% Extract image patches
img_samples = cell(2,1);
for sz_ind = 1:num_sizes
    img_sample_sz = extract_info.img_sample_sizes{sz_ind};
    img_input_sz = extract_info.img_input_sizes{sz_ind};
    img_samples{sz_ind} = zeros(img_input_sz(1), img_input_sz(2), size(image,3), num_scales, 'uint8');
    for scale_ind = 1:num_scales
        img_samples{sz_ind}(:,:,:,scale_ind) = sample_patch(image, pos, round(img_sample_sz*scales(scale_ind)), img_input_sz, gparams);
    end
end

% Find the number of feature blocks and total dimensionality
num_feature_blocks = 0;
total_dim = 0;
for feat_ind = 1:num_features
    num_feature_blocks = num_feature_blocks + length(features{feat_ind}.fparams.nDim);
    total_dim = total_dim + sum(features{feat_ind}.fparams.nDim);
end

feature_map = cell(1, 1, num_feature_blocks);

% Extract feature maps for each feature in the list
ind = 1;
for feat_ind = 1:num_features
    feat = features{feat_ind};
    gparams.cell_size = feat.fparams.cell_size;

    % get the image patch index
    img_sample_ind = cellfun(@(sz) isequal(feat.img_sample_sz, sz), extract_info.img_sample_sizes);

    % do feature computation
    if feat.is_cell
        num_blocks = length(feat.fparams.nDim);
        feature_map(ind:ind+num_blocks-1) = feat.getFeature(img_samples{img_sample_ind}, feat.fparams, gparams);
    else
        num_blocks = 1;
        feature_map{ind} = feat.getFeature(img_samples{img_sample_ind}, feat.fparams, gparams);
    end

    ind = ind + num_blocks;
end

% Do feature normalization
if ~isempty(gparams.normalize_power) && gparams.normalize_power > 0
    if gparams.normalize_power == 2
        feature_map = cellfun(@(x) bsxfun(@times, x, ...
            sqrt((size(x,1)*size(x,2))^gparams.normalize_size * size(x,3)^gparams.normalize_dim ./ ...
            (sum(reshape(x, [], 1, 1, size(x,4)).^2, 1) + eps))), ...
            feature_map, 'uniformoutput', false);
    else
        feature_map = cellfun(@(x) bsxfun(@times, x, ...
            ((size(x,1)*size(x,2))^gparams.normalize_size * size(x,3)^gparams.normalize_dim ./ ...
            (sum(abs(reshape(x, [], 1, 1, size(x,4))).^gparams.normalize_power, 1) + eps)).^(1/gparams.normalize_power)), ...
            feature_map, 'uniformoutput', false);
    end
end
if gparams.square_root_normalization
    feature_map = cellfun(@(x) sign(x) .* sqrt(abs(x)), feature_map, 'uniformoutput', false);
end
end

repo: CFWCR


text:
function s=norm(a,P)
%norm: Matrix or vector norm of vector/matrix of mp-type

%     For matrices...
%       NORM(X) is the largest singular value of X, max(svd(X)).
%       NORM(X,2) is the same as NORM(X).
%       NORM(X,1) is the 1-norm of X, the largest column sum,
%                       = max(sum(abs(X))).
%       NORM(X,inf) is the infinity norm of X, the largest row sum,
%                       = max(sum(abs(X'))).
%       NORM(X,'fro') is the Frobenius norm, sqrt(sum(diag(X'*X))).
%       NORM(X,P) is available for matrix X only if P is 1, 2, inf or 'fro'.
%
%     For vectors...
%       NORM(V,P) = sum(abs(V).^P)^(1/P).
%       NORM(V) = norm(V,2).
%       NORM(V,inf) = max(abs(V)).
%       NORM(V,-inf) = min(abs(V)).
%
%     See also COND, RCOND, CONDEST, NORMEST.

if nargin<2,P=2;end
if isa(P,'mp')
 error(['Second parameter should not be of mp-type'])
end
if min(size(a))==1
 %is a vector
 if isinf(P) | strcmp(lower(P),'inf')
  if P>0
   s=max(abs(a));
  else
   s=min(abs(a));
  end
  return
 end
 A=max(abs(a(:)));a=a/A;
 if P==2
  %naive implementation
  s=A*sqrt(sum(a.*a));
 else
  s=A*(sum(a.^P))^(1/P);
 end
else
 %is a matrix
 if isinf(P) | strcmp(lower(P),'inf')
  s=max(sum(abs(a')));
  return
 end
 if ischar(P)
  if strcmp(lower(P),'fro')
   s=sqrt(sum(diag(a'*a)));
  else
   error(['Second parameter is not recognized as an option: ' P])
  end
  return
 end
 switch P
  case 1
   s=max(sum(abs(a)));
  case 2
   s=max(svd(a));
  otherwise
   error(['NORM(X,P) is available for matrix X only if P is 1, 2, inf or ' setstr(39) 'fro' setstr(39) '.'])
 end
end

author: opencobra

repo: cobratoolbox


text:
% twa_demo.m, The following script demonstrates how to use the TWA code in
% the MVtoolbox to estimate TWAs using the modified moving average (MMA)
% method coupled with a non parametric surrogate statistical reshuffling
% method to reduce false positive TWA detections. The algorithm is
% referenced as follows,
%
%   Nemati S, Abdala O, Monasterio V, Yim-Yeh S, Malhotra A, Clifford GD.
%   A nonparametric surrogate-based test of significance for T-wave alternans detection.
%   IEEE Trans Biomed Eng. 2011;58(5):1356–1364. doi:10.1109/TBME.2010.2047859
%
% The main function is ComputeTWA. The output is stored in the TWAResult
% structure with the following fields with the expected output listed as follows,
%
%   Field                   Expected values
%   TWAResult.HR            [80.32,80.21,79.79,80.65,80.65,80.43,80,80,80.21,80.32,80,79.47]
%   TWAResult.VAlt          [38.84,40.10,37.81,37.38,37.70,38.14,39.32,39.74,38.31,40.39,39.07,42.20]
%   TWAResult.VAlt_Sig      [38.84,40.10,37.81,37.38,37.70,38.14,39.32,39.74,38.31,40.39,39.07,42.20]
%   TWAResult.Noise_Median  [4.28,3.92,4.40,3.75,3.68,4.14,3.00,3.21,4.15,3.95,3.78,7.20]
%   TWAResult.Noise_95      [11.19,10.75,9.61,12.67,11.69,9.72,9.15,9.84,10.29,11.03,9.89,19.40]
%   TWAResult.VAltPt        [132,116,115,125,116,122,130,117,91,111,101,159]
%   TWAResult.successful    1
%
% REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Ismail Sadiq
%	COPYRIGHT (C) 2019
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information. The license may be found in
%       the Documents folder of the Physionet-Cardiovascular-Signal-Toolbox.

% add dependencies within cardiovascular toolbox
addpath(genpath('../Tools/'));
addpath(genpath('../../../../../PhysioNet-Cardiovascular-Signal-Toolbox-master/'));
% load data and annotations
sig = load('./testdata/twa/sample_ecg_twam');
siginfo = readheader('./testdata/twa/sample_ecg_twam.hea');
Fs = siginfo.freq;
ecg = (sig.val - siginfo.adczero)./siginfo.gain; % Adjust signal according to gain and dc offset
clear sig siginfo
if (size(ecg,2) > size(ecg,1))
   ecg = ecg';
end

% load annotations
ann.QRSon = read_ann('./testdata/twa/sample_ecg_twa','qrson')'; ann.Q = read_ann('./testdata/twa/sample_ecg_twa','q')'; ann.R = read_ann('./testdata/twa/sample_ecg_twa','r')'; ann.S = read_ann('./testdata/twa/sample_ecg_twa','s')'; ann.QRSoff = read_ann('./testdata/twa/sample_ecg_twa','qrsoff')'; ann.Toff = read_ann('./testdata/twa/sample_ecg_twa','toff')';

% TWA evaluation
TWAResult = ComputeTWA(ecg,ann,Fs);

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function testFinalModels_HN_OldExp(pathExperiments,nExp,fSetNames,paramCells,nameCells)
% - paramCells: One cell of parameters for each fSet
% - nameCells: One cell of names for each fSet

startpath = pwd;
nFset = numel(fSetNames);
for exp = 1:nExp
    cd(fullfile(pathExperiments,['Experiment',num2str(exp)])), pathExperiment = pwd;
    load('training'), load('testing')
    nameOutcomes = fieldnames(training); nOutcomes = numel(nameOutcomes);
    cd('FINAL_MODELS'), pathFinalModels = pwd;
    for o = 1:nOutcomes
        for f = 1:nFset
            results = []; results = struct;
            cd(fullfile(pathFinalModels,nameOutcomes{o},fSetNames{f}))
            load('finalModel'), load('coeff'), load('response'), model = struct;

            % OBTAINING MODEL CHARACTERISTIC
            order = numel(finalModel.Name);
            paramCell = paramCells{f}; nameCell = nameCells{f};
            nParam = numel(paramCell); scans = cell(1,order);
            for n = 1:order
                name = finalModel.Name{n};
                indA = strfind(name,'--'); indB = strfind(name,'-');
                if ~isempty(indA) % This is a texture
                    model.(['Feature',num2str(n)]).isTexture = true;
                    model.(['Feature',num2str(n)]).textType = name((indA(1)+2):(indB(3)-1));
                    model.(['Feature',num2str(n)]).textName = name((indB(3)+1):end);
                    indScan = strfind(name,'(');
                    scans{n} = name(1:indScan-1);
                    model.(['Feature',num2str(n)]).scan = scans{n};
                    indMat = zeros(1,nParam); sizeParam = zeros(1,nParam);
                    for p = 1:nParam
                        paramAll = paramCell{p}; nSpecParam = numel(paramAll); sizeParam(p) = nSpecParam;
                        nameParam = nameCell{p};
                        indP = strfind(name,nameParam);
                        for start = indP:numel(name)
                            if strcmp(name(start),'=')
                                break
                            end
                        end
                        for final = indP:numel(name)
                            if strcmp(name(final),',') | strcmp(name(final),')')
                                break
                            end
                        end
                        if iscell(paramCell{p})
                            param = name((start+1):(final-1));
                            for nsp = 1:nSpecParam
                                if strcmp(param,paramAll{nsp})
                                    indMat(p) = nsp;
                                    break
                                end
                            end
                        else
                            param = str2num(name((start+1):(final-1)));
                            for nsp = 1:nSpecParam
                                if param > (paramAll(nsp)-0.05) & param < (paramAll(nsp)+0.05)
                                    indMat(p) = nsp;
                                    break
                                end
                            end
                        end
                    end
                    if numel(sizeParam) == 2
                        model.(['Feature',num2str(n)]).indText = sub2ind(sizeParam,indMat(1),indMat(2));
                    elseif numel(sizeParam) == 3
                        model.(['Feature',num2str(n)]).indText = sub2ind(sizeParam,indMat(1),indMat(2),indMat(3));
                    elseif numel(sizeParam) == 4
                        model.(['Feature',num2str(n)]).indText = sub2ind(sizeParam,indMat(1),indMat(2),indMat(3),indMat(4));
                    elseif numel(sizeParam) == 5
                        model.(['Feature',num2str(n)]).indText = sub2ind(sizeParam,indMat(1),indMat(2),indMat(3),indMat(4),indMat(5));
                    elseif numel(sizeParam) == 6
                        model.(['Feature',num2str(n)]).indText = sub2ind(sizeParam,indMat(1),indMat(2),indMat(3),indMat(4),indMat(5),indMat(6));
                    end
                    model.(['Feature',num2str(n)]).FullName = name;
                else % This is a non-texture
                    scans{n} = 'nonText';
                    model.(['Feature',num2str(n)]).isTexture = false;
                    model.(['Feature',num2str(n)]).nonText = name;
                end
            end
            model.coeff = coeff;
            model.order = order;
            results.model = model;
            results.trainData.data = finalModel.Data; results.trainData.response = response; results.trainData.outcome = training.(nameOutcomes{o}).outcome;

            % TESTING MODEL
            text = testing.(nameOutcomes{o}).text; nonText = testing.(nameOutcomes{o}).nonText;
            outcome = testing.(nameOutcomes{o}).outcome;
            resp = zeros(numel(outcome),1);
            data = zeros(numel(outcome),order);
            for n = 1:order
                coeff = model.coeff(n);
                feat = ['Feature',num2str(n)];
                scan = scans{n};
                if model.(feat).isTexture
                    data(:,n) = text.(fSetNames{f}).(scan){model.(feat).indText}.(model.(feat).textType).(model.(feat).textName).Data;
                else
                    data(:,n) = nonText.(model.(feat).nonText).Data;
                end
                resp = resp + coeff*data(:,n);
            end
            resp = resp + model.coeff(end);
            results.testData.data = data; results.testData.response = resp; results.testData.outcome = outcome;
            [AUC,sensitivity,specificity,accuracy] = calcPerformMetrics(resp,outcome,0);
            results.AUC = AUC;
            results.Sensitivity = sensitivity;
            results.Specificity = specificity;
            results.Accuracy = accuracy;
            cd(pathExperiment), save(['testResultsModels_',fSetNames{f},'_',nameOutcomes{o}],'results')
        end
    end
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
% Predict the label for a trained one-vs-all classifier. The labels
% are in the range 1..K, where K = size(all_theta, 1)
function p = one_vs_all_predict(all_theta, X)
    m = size(X, 1);
    num_labels = size(all_theta, 1);

    % We need to return the following variables correctly.
    p = zeros(size(X, 1), 1);

    % Add ones to the X data matrix
    X = [ones(m, 1) X];

    % Calculate probabilities of each number for each input example.
    % Each row relates to the input image and each column is a probability that this example is 1 or 2 or 3 etc.
    h = sigmoid(X * all_theta');

    % Now let's find the highest predicted probability for each row.
    % Also find out the row index with highest probability since the index is the number we're trying to predict.
    [p_vals, p] = max(h, [], 2);
end

author: trekhleb

repo: machine-learning-octave


text:
%kkbitorinv 'Output = NOT(Input 1) OR (Input 2 or Constant)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kbitorinv.pane file
%
% Parameters:
% InputFile: i1 'Input 1', required: 'First input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
% InputFile: i2 'Input 2', optional: 'Second input data object'
% Integer: real 'Constant', default: 0: 'Constant value'
%
% Example: o = kkbitorinv({i1, i2}, {'i1','';'o','';'i2','';'real',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kbitorinv - Output = NOT(Input 1)  OR  (Input 2 or Constant)
%
%  DESCRIPTION
% The bitwise OR Inverted operator performs a bitwise NOT on \fBInput 1"
% and then a bitwise OR between the result of the NOT operation and
% either \fBInput 2" or the \fBConstant\fP value, which
% ever is specified by the user.
%
%  "Data Type - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "Data Type - Two Input Objects" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_2input
% The bitwise operators will not operate on float or complex data types.
%
%  "Map Data" 5
% The bitwise operations have not been written to be fully polymorphic yet.
% They does not
% check for map data, and will therefore always operate on the value data,
% even if a map exists.  This will most likely corrupt indexing into the map.
% In the case of a single input, it is recommended to use the "Copy to
% Value" (kcptoval) segment operator to temporarily move the map data into
% the value segment, run the bitwise operation on the data, then move it
% back to the map with the "Copy from Value" (kcpfromval) operator.
% When operating with two input objects, where at least one of the objects
% contains map data, data should be mapped prior to running the bitwise
% operator.  The "Map Data" operator (kmapdata) can be used to perform
% the mapping operation.
%
%  "Validity Mask - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "Validity Mask - Two Input Objects" 5
% Masking has not been implemented yet for the bitwise operators.  Therefore,
% only the mask from the first input object will be transferred to the output.
%
%  "Input Objects of Different Sizes" 5
% .cI $DATAMANIP/repos/shared/man/sections/resize_2input
% The value used to pad the data when the input files are not the same size is
% zero.
%
%  "Explicit Location and Time Data - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "Explicit Location and Time Data - Two Input Objects" 5
% The bitwise operations have not been extended to understand location and
% time data.  Therefore, only location and time data present in the first input
% object will be transferred to the output.
%
% Executing \fBOR Inverted" runs the program "kbitwise\fP with
% the -orinv flag.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::kbitwise
%
%  RESTRICTIONS
%
% THIS ROUTINE is still under development, and has not been modified to fully
% support the polymorphic data model.  See paragraphs above for discussions
% concerning the polymorphic data segments.
%
% All processing is currently performed as UNSIGNED long, so the operations
% may not work properly on negative values.
%
% The bitwise operators will not operate on float or complex data types.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkbitorinv(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkbitorinv(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'o', '__output';'i2', '__input';'real', 0};
maxval={0,0,1,0};
minval={0,0,1,0};
istoggle=[0,0,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','InputFile','Integer'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kbitwise"  -orinv'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB