text:
function plotDiff(SO3F1,SO3F2,varargin)
% difference plot between two SO3Funs or an SO3Fun and a pole figure
%
% Syntax
%   plotDiff(SO3F1,SO3F2,...,param,val,...)
%   plotDiff(SO3F1,pf,...,param,val,...)
%
% Input
%  SO3F1,SO3F2  - @SO3Fun
%  pf   - @PoleFigure
%
% Options
%  RP - calculate RP error (only for SO3Fun -- pole figure)
%  l1 - calculate $|pf1--pf2|$ error (only for SO3Fun -- pole figure)
%  l2 - calculate $|pf1--pf2|^2$ error (only for SO3Fun -- pole figure)
%
% See also
% S2Grid/plot PoleFigure/calcError SO3Fun/calcError savefigure
% Plotting Annotations_demo ColorCoding_demo PlotTypes_demo
% SphericalProjection_demo

if isa(SO3F2,'PoleFigure')
  plot(calcErrorPF(SO3F2,SO3F1,varargin{:}),'colorrange','equal',varargin{:})
else
  plotSection(SO3F1-SO3F2,varargin{:})
end

author: mtex-toolbox

repo: mtex

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


text:
function [all_theta] = oneVsAll(X, y, num_labels, lambda)
%ONEVSALL trains multiple logistic regression classifiers and returns all
%the classifiers in a matrix all_theta, where the i-th row of all_theta
%corresponds to the classifier for label i
%   [all_theta] = ONEVSALL(X, y, num_labels, lambda) trains num_labels
%   logisitc regression classifiers and returns each of these classifiers
%   in a matrix all_theta, where the i-th row of all_theta corresponds
%   to the classifier for label i

% Some useful variables
m = size(X, 1);
n = size(X, 2);

% You need to return the following variables correctly
all_theta = zeros(num_labels, n + 1);

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

% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the following code to train num_labels
%               logistic regression classifiers with regularization
%               parameter lambda.
%
% Hint: theta(:) will return a column vector.
%
% Hint: You can use y == c to obtain a vector of 1's and 0's that tell use
%       whether the ground truth is true/false for this class.
%
% Note: For this assignment, we recommend using fmincg to optimize the cost
%       function. It is okay to use a for-loop (for c = 1:num_labels) to
%       loop over the different classes.
%
%       fmincg works similarly to fminunc, but is more efficient when we
%       are dealing with large number of parameters.
%
% Example Code for fmincg:
%
%     % Set Initial theta
%     initial_theta = zeros(n + 1, 1);
%
%     % Set options for fminunc
%     options = optimset('GradObj', 'on', 'MaxIter', 50);
%
%     % Run fmincg to obtain the optimal theta
%     % This function will return theta and the cost
%     [theta] = ...
%         fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), ...
%                 initial_theta, options);
%

% Set Initial theta
initial_theta = zeros(n + 1, 1);

% Set options for fminunc
options = optimset('GradObj', 'on', 'MaxIter', 50);
for c = 1:num_labels
  [theta] = fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), initial_theta, options);
  all_theta(c, :) = theta';
end

% =========================================================================

end

author: benoitvallon

repo: coursera-machine-learning


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nickabattista@gmail.com
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs, non-invariant beams*)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "Hill+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nickabattista@gmail.com) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the spring attributes!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function springs_info = update_Springs(dt,current_time,xLag,yLag,springs_info)

%springs_info: col 1: starting spring pt (by lag. discretization)
%              col 2: ending spring pt. (by lag. discretization)
%              col 3: spring stiffness
%              col 4: spring resting lengths

%RL = springs_info(:,4); % resting-length vector

% Contraction Frequency
freq = 2.0;

% CHANGE RESTING LENGTH BTWN SIDES OF JELLYFISH BELL 1
springs_info(319:357,4) = abs( cos(freq*current_time*pi) );

% CHANGE RESTING LENGTH BTWN SIDES OF JELLYFISH BELL 2
%springs_info(676:714,4) = abs( cos(freq*current_time*pi) );

author: nickabattista


text:
classdef constrained_quadratic
% This file defines constrained quadratic problem class
%
% Inputs:
%           Q       a positive definite matrix of size nxn
%           p       a column vector of size n
%           A       a matrix of size nxm
%           b       a column vector of size m
%           G       a matrix of size nxp
%           h       a column vector of size p
%
% Output:
%           Problem     problem instance.
%
%
% The problem of interest is defined as
%
%           min f(x) = 1/2 * x^T * Q * x - p^T * x.
%           s.t.
%           A * x = b, G * x <= h
%
% "w" is the model parameter of size d vector.
%
%
% This file is part of GDLibrary.
%
% Created by H.Kasai on Oct. 22, 2020
% Modified by H.Kasai on Dec. 13, 2022

    properties
        d;
        n;
        name;
        dim;
        samples;
        Q;
        p;
        A;
        b;
        G;
        h;

        prox_flag;

    end

    methods
        function obj = constrained_quadratic(Q, p, A, b, G, h)
            obj.d = length(p);
            obj.n = obj.d;

            obj.name = 'constrained_quadratic';
            obj.dim = obj.d;
            obj.samples = obj.n;
            obj.Q = Q;
            obj.p = p;
            obj.A = A;
            obj.b = b;
            obj.G = G;
            obj.h = h;

            obj.prox_flag = true;
        end

        function v = prox(obj, w, t)
            % project y onto the feasible set Ax <= b, i.e. nearest feasible point x to y
            % minimize 0.5||x-y||^2
            % s.t.     Ax <= b

            % equivalent to
            % min J   = 0.5 x'x - x'y + 0.5 y'y
            % st. Ax <= b

            %n = length(y);
            H = eye(obj.d);
            w0 = w;
            qp_options = optimset('Display','off');
            v = quadprog(H,-w,obj.G,obj.h,[],[],[],[],w0, qp_options);
        end

        function f = cost(obj, x)

            f = 1/2 * x' * obj.Q * x - obj.p' * x;
        end

        function g = grad(obj, x)

            g = obj.Q * x - obj.p;

        end

        function g = full_grad(obj, x)

            g = obj.grad(x);

        end

        function h = hess(obj, x)

            h = obj.Q;

        end

        function h = full_hess(obj, x)

            h = obj.hess(x);

        end
    end

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function write_off(filename, pnt, plc)

% WRITE_OFF writes a set of geometrical planar forms (called piecewise linear complex, PLC)
% to an ascii *.off file, which is a file format created by Princeton Shape Benchmark
%
% Use as
%   write_stl(filename, pnt, tri)
%
% See also READ_OFF

% Copyright (C) 2010, Cristiano Micheli
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: write_off.m 2212 2010-11-27 11:55:07Z roboos $

nedges = 0;
fid  = fopen(filename, 'wb');
npnt = size(pnt,1);
nplc = size(plc,1);

fprintf(fid, 'OFF\n');
fprintf(fid, '%d %d %d\n',npnt,nplc,nedges);
for i=1:npnt
  fprintf(fid, '%f %f %f\n', pnt(i,1), pnt(i,2), pnt(i,3));
end
for i=1:nplc
  str = '%d ';
  nvert = size(plc(i,:),2);
  str2 = 'nvert,';
  for j=1:nvert
    str  = [str '%d '];
    str2 = [str2 'plc(i,' num2str(j) '),'];
  end
  str = [str '\n''']; str2 = str2(1:end-1);
  eval(['fprintf(fid,''' str ',' str2 ');']);
end

fclose(fid);

author: PatternRecognition

repo: OpenBMI


text:
function comp_infos = normal_mser2(src,info)
%this script use the normal version of mser extract regions
%default: one channel
img = imadjust(rgb2gray(src));
%mser
[bbox,pixelsList,bright_dark_flag] = mser(img,info);
%record
comp_infos={};
for kk=1:size(bbox,1)
    box=bbox(kk,:);
    comp_infos(kk).box=box;
    comp_infos(kk).center=floor([box(1)+box(3)/2,box(2)+box(4)/2]);
    comp_infos(kk).pixelList = pixelsList{kk,1};
    comp_infos(kk).bright_dark_flag = bright_dark_flag(kk);
end

%%%%%%%% function %%%%%%%%
function  [bbBox,pixelsList,bright_dark_flag] = mser(I,info)
 %
[h,w]=size(I);
%phase-1: 'BrightOnDark'
[r,f] = vl_mser(I,'MinDiversity',info.minDiversity,'MaxVariation',info.maxVariation,...
				'MaxArea',info.maxArea,'MinArea',info.minArea,'BrightOnDark',1,'DarkOnBright',0,...
				'Delta',info.delta);
M =zeros(size(I));
for x=r'
    s = vl_erfill(I,x);
    M(s) = M(s) + 1;%if s belong to one region, add 1.
end
mCount = max(max(M));
bbBox = [];
pixelsList = {};
bright_dark_flag=[];
nCount = 1;
for ii=mCount:-1:1
    MM = M;
    MM(find(MM<ii))=0;
    MM(find(MM>=ii))=1;
    mContours = bwlabel(MM,8);
    mNum = max(max(mContours));
    for jj=1:mNum
        [idx,idy,~] = find(mContours == jj);
        ww = max(idx) - min(idx) + 1;%note me: ww & hh are reverse
        hh = max(idy) - min(idy) + 1;
		%condition
		if (ww/(hh+eps)<0.3)||(hh/(ww+eps)<0.3)
			continue;
		end
		if (ww<4||hh<4) || (ww>400|| hh>400)
			continue;
        end
        if ww*hh< 20
            continue;
        end
		%update
		bbBox(nCount,:)=[min(idx),min(idy),ww,hh];
		pixelsList{nCount,1}=[idx,idy];
        bright_dark_flag(nCount)=1;
		nCount = nCount + 1;
    end
end

%phase-2: 'DarkOnBright'
[r,f] = vl_mser(I,'MinDiversity',info.minDiversity,'MaxVariation',info.maxVariation,...
				'MaxArea',info.maxArea,'MinArea',info.minArea,'BrightOnDark',0,'DarkOnBright',1,...
				'Delta',info.delta);
M =zeros(size(I));
for x=r'
    s = vl_erfill(I,x);
    M(s) = M(s) + 1;%if s belong to one region, add 1.
end
mCount = max(max(M));
for ii=mCount:-1:1
    MM = M;
    MM(find(MM<ii))=0;
    MM(find(MM>=ii))=1;
    mContours = bwlabel(MM,8);
    mNum = max(max(mContours));
    for jj=1:mNum
        [idx,idy,~] = find(mContours == jj);
        ww = max(idx) - min(idx) + 1;%note me: ww & hh are reverse
        hh = max(idy) - min(idy) + 1;
		%condition
		if (ww/(hh+eps)<0.3)||(hh/(ww+eps)<0.3)
			continue;
		end
		if (ww<4||hh<4) || (ww>400|| hh>400)
			continue;
        end
        if ww*hh< 20
            continue;
        end
		%update
		bbBox(nCount,:)=[min(idx),min(idy),ww,hh];
		pixelsList{nCount,1}=[idx,idy];
        bright_dark_flag(nCount)=0;
		nCount = nCount + 1;
    end
end
end

end

author: stupidZZ


text:
% DeInterleaver function
function [data_out]=deinterleav_d(data_in,rate_id)

switch (rate_id)
    case 0
       Ncbps=192;  % In BPSK no.of bits allocated to subcarrier in a OFDM symbol
       Ncpc=1;     % In cash of BPSK No. of coded bit per carrier.
    case {1,2}
      Ncbps=384;
      Ncpc=2;
    case {3,4}     %% note: no. of coded bit Ncpc=1,2,4,6 for bpsk,qpsk,16qam,64qam respectivly%%
       Ncbps=768;
       Ncpc=4;
    case {5,6}
       Ncbps=1152;
       Ncpc=6;
    otherwise
       display('error in interleaver give proper rate_id')
end

s=ceil(Ncpc/2);

data=data_in';
%%% Deinterleaving
  % j-->index of the bit encoded BEFORE the first permutation
  % mj-->index of this bit BEFORE the second permutation and AFTER the first one
  % kj-->index after the SECOND permutation, just before the mapping of the sign.

  j = 0:Ncbps-1;
  mj = s*floor(j/s) + mod((j + floor(12*j/Ncbps)),s);          % First permutation
  kj = 12*mj-(Ncbps-1)*floor(12*mj/Ncbps);                     % Second permutation

  % The indices are ordered to know what must be taken.
  [c d]= sort(kj);

% finally the bits are rearranged.

 i = 1:Ncbps;
 data_out = zeros(1,Ncbps);
 data_out(i) = data(d(i));
 data_out=data_out';

author: Sable

repo: mcbench-benchmarks


text:
function[imdb] = setupImdbClassification(trainName, testName, net)
% [imdb] = setupImdbClassification(trainName, testName, net)

global DATAopts;

% Setup the Imdb
% Get images and labels
set = trainName;
[trainIms, ~] = textread(sprintf(DATAopts.imgsetpath,set),'%s %d'); %#ok<DTXTRD>
trainLabs = zeros(length(trainIms), DATAopts.nclasses);

% Get all testlabels
if ~strcmp(set,'test') || DATAopts.year == 2007
    for classIdx = 1 : DATAopts.nclasses
        class = DATAopts.classes{classIdx};
        [~, trainLabs(:,classIdx)] = textread(sprintf(DATAopts.clsimgsetpath, class, set), '%s %d'); %#ok<DTXTRD>
    end

    trainLabs(trainLabs == 0) = -1;
end

set = testName;
[testIms, ~] = textread(sprintf(DATAopts.imgsetpath,set),'%s %d'); %#ok<DTXTRD>
testLabs = zeros(length(testIms), DATAopts.nclasses);

% Get all testlabels
if ~strcmp(set,'test') || DATAopts.year == 2007
    for classIdx = 1 : DATAopts.nclasses
        class = DATAopts.classes{classIdx};
        [~, testLabs(:, classIdx)] = textread(sprintf(DATAopts.clsimgsetpath, class, set), '%s %d'); %#ok<DTXTRD>
    end

    testLabs(testLabs == -1) = 0;
end

% Make train, val, and test set.
numValIms = 500;
allIms = cat(1, trainIms, testIms);
allLabs = cat(1, trainLabs, testLabs);
datasetIdx{1} = (1 : length(trainIms) - numValIms)';  % Last numValIms are val set
datasetIdx{2} = (length(trainIms) - numValIms + 1 : length(trainIms))';
datasetIdx{3} = (length(trainIms) + 1 : length(allIms))';

% Setup the classification imdb
imdb = ImdbClassification(DATAopts.imgpath(1:end-6), ...        % path
    DATAopts.imgpath(end-3:end), ...      % image extension
    allIms, ...                           % all images
    allLabs, ...                          % all labels
    datasetIdx, ...                       % division into train/val/test
    net.meta.normalization.averageImage, ...   % average image used to pretrain network
    20);                                  % num classes

% Store lists for use in eval
imdb.misc.trainName = trainName;
imdb.misc.testName  = testName;
imdb.misc.trainIms  = trainIms;
imdb.misc.trainLabs = trainLabs;
imdb.misc.testIms   = testIms;
imdb.misc.testLabs  = testLabs;

author: nightrome

repo: matconvnet-calvin


text:
function bytestream=savenifti(img, filename, varargin)
%
%    savenifti(img, filename)
%        or
%    savenifti(img, filename, rawhdr)
%    savenifti(img, filename, 'nifti2')
%    bytestream=savenifti(img)
%
%    Write an image to a NIfTI (*.nii) or compressed NIfTI file (.nii.gz)
%
%    author: Qianqian Fang (q.fang <at> neu.edu)
%
%    input:
%        img: this is a numerical array to be stored in the NIfTI file
%        filename: output file name, can have a suffix of '.nii' or '.nii.gz'
%                 if a .gz suffix is used, this function needs the JSONLab
%                 (http://gitlab.com/fangq/jsonlab) and ZMat (http://gitlab.com/fangq/zmat)
%                 to perform the compression.
%        rawhdr (optional): a struct, as a pre-created/loaded NIfTI header data structure
%                if rawhdr is 'nifti1' or 'nifti2', this function calls
%                nifticreate to create a default header.
%    output:
%        bytestream (optional): the output file byte stream. it only returns this output if
%                no filename is given.
%
%    example:
%        a=single(rand(10,20,30));
%        savenifti(a,'randnii.nii');
%        savenifti(a,'randnii2.nii.gz','nifti2'); % needs zmat
%
%
%    this file is part of JNIfTI specification: https://github.com/fangq/jnifti
%
%    License: Apache 2.0, see https://github.com/fangq/jnifti for details
%

if(~isempty(varargin))
    if(isstruct(varargin{1}))
        header=varargin{1};
    elseif(ischar(varargin{1}))
        header=nifticreate(img,varargin{1});
    end
else
    header=nifticreate(img);
end

names=fieldnames(header);
buf=[];
for i=1:length(names)
    buf=[buf,typecast(header.(names{i}),'uint8')];
end

if(length(buf)~=352 && length(buf)~=544)
    error('incorrect nifti-1/2 header %d',length(buf));
end

buf=[buf,typecast(img(:)','uint8')];

if(nargout>1 && nargin<2)
    bytestream=buf;
    return;
end

if(regexp(filename,'\.[Gg][Zz]$'))
    buf=gzipencode(buf);
end

fid=fopen(filename,'wb');
if(fid==0)
    error('can not write to the specified file');
end
fwrite(fid,buf);
fclose(fid);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [pos,picon,picoff,janelas,messages]=twave3D(w,timenew,i,freq,S,QRSoff,samp,rrmed,messages)
%
%[pos,picon,picoff,janelas]=twave3D(w,timenew,i,freq,S,QRSoff,samp,rrmed)
%
% multilead T wave delineation
%
%Input Parameters:
%   w: matrix with WT scales 1 to 5
%   timenew:  QRS times in inedexes refering the interval included in the current excerpt (borders excluded)
%   i: beat number
%   freq: sampling frequency (heasig.freq)
%   S: lattest S wave position found in any lead
%   QRSoff: lattest QRS end position found in any lead
%   samp: samples included in the current excerpt (borders excluded)
%   rrmed: Exponentially averaged RR
%
%Output Parameters:
%   pos: fiducial marks structure (position)
%   picon: position of the first relevant modulos maximum in the wavelet
%   picoff: position of the last relevant modulos maximum in the wavelet
%   janelas: T wave seach windows
%
% Rute Almeida
% Last update: Rute Almeida  05AGO2011
%
% Designed for MATLAB Version R12; tested with MATLAB Version R13
%

%%%%%% Constants and Thresholds  !!!!!!!!!!!!!!!!!!!!!!!!!
if ~isfield(messages.setup.wavedet,'umbraldetT')
    messages.setup.wavedet.umbraldetT = 0.25;  % We use umbraldet*sqrt(mean(w(time(i):time(i+1),4).^2))
end
if ~isfield(messages.setup.wavedet,'umbralsig')
    messages.setup.wavedet.umbralsig  =  1/8;
end
if ~isfield(messages.setup.wavedet,'Kton')
    messages.setup.wavedet.Kton = 4;      % 2 4!
end
if ~isfield(messages.setup.wavedet,'inivent_tol')
    messages.setup.wavedet.inivent_tol  =  0.1;
end
if ~isfield(messages.setup.wavedet,'inivent_tol_S')
    messages.setup.wavedet.inivent_tol_S=0.05; % sec
end
if ~isfield(messages.setup.wavedet,'finvent_tol')
    messages.setup.wavedet.finvent_tol=0.240;% sec
end
if ~isfield(messages.setup.wavedet,'finvent_max')
    messages.setup.wavedet.finvent_max=0.6;% sec
end
if ~isfield(messages.setup.wavedet,'min_vent')
    messages.setup.wavedet.min_vent=0.1;%sec
end
if ~isfield(messages.setup.wavedet,'Tmax_Tmin_time_min')
    messages.setup.wavedet.Tmax_Tmin_time_min=0.15;%sec
end
if ~isfield(messages.setup.wavedet,'Tmax_Tmin_bifasic')
    messages.setup.wavedet.Tmax_Tmin_bifasic=2.5;
end

Kton=messages.setup.wavedet.Kton;
Ktoff=messages.setup.wavedet.Kton;
messages.warnings=[messages.warnings {'Recall that in twave3D ktoff=kton.'}];
umbraldetT = messages.setup.wavedet.umbraldetT;
umbralsig = messages.setup.wavedet.umbralsig;
inivent_tol= messages.setup.wavedet.inivent_tol;
inivent_tol_S=messages.setup.wavedet.inivent_tol_S;
finvent_tol=messages.setup.wavedet.finvent_tol;
finvent_max= messages.setup.wavedet.finvent_max;
min_vent=messages.setup.wavedet.min_vent;
Tmax_Tmin_time_min= messages.setup.wavedet.Tmax_Tmin_time_min;%sec
Tmax_Tmin_bifasic=messages.setup.wavedet.Tmax_Tmin_bifasic;% extra criteria
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
janelas=[];%%31May05
% Initialization of auxiliary variables
T = []; Tprima=[]; picon=[]; picoff=[]; Ton=[]; Toff=[]; tipoT=[];
% minapos = []; minppos=[]; maxapos=[]; maxppos=[]; mina=[]; minp=[]; maxa=[];
% maxp=[];
picon_keep=[];%multileadchange
picoff_keep=[];%multileadchange
lead_keep=[];%multileadchange

inivent = round(inivent_tol*freq);   % Begining of window
if ~isempty(S),              % If there is an S wave
    inivent = max(inivent, S-timenew(i)+round(inivent_tol_S*freq));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% changed 13/06/02 Rute
if rrmed <= freq,
    %if rrmed >= freq,
    finvent = round(finvent_max*freq);     % End of window
else
    finvent = round(rrmed*finvent_max);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %% changed 13/06/02 Rute

% using timenew from diffferent leads ty ca dist less tahn 0.24 sec!!!
if i~=length(timenew),                 % For last beat in the segment
    finvent = min(finvent,timenew(i+1)-timenew(i)-round(finvent_tol*freq));
else
    finvent = min(finvent,timenew(i)-timenew(i-1)-round(finvent_tol*freq));
end

% We work at scale scale (in general).
if isempty(QRSoff),   % Should never happen, but...
    % begwin = inivent + timenew(i);
    begwin = inivent;
else
    %  begwin = max(inivent + timenew(i), QRSoff+1);
    begwin = max(inivent , QRSoff+1-timenew(i));
end
%endwin = min(finvent + timenew(i),length(w)); %% Rute 20/05/02
endwin = min(finvent ,length(w)); %% Rute 20/05/02
janelas=[janelas; i begwin endwin];%31May05
if ~isempty(begwin+1:endwin ) % 22.04.05

    % Positive maxima and negative minima in the window
    maxpos = begwin + modmax(w(begwin+1:endwin ),2,0,+1);
    minpos = begwin + modmax(w(begwin+1:endwin ),2,0,-1);
    [maxim ind] = max(w(maxpos ));   % The biggest of the positive
    maxpos = maxpos(ind);
    [minim ind] = min(w(minpos ));   % The biggest of the negative
    minpos = minpos(ind);

    if isempty(maxpos),               % If no local positive maximum
        % The maximum will be the first
        % or the last sample
        if (w(begwin )>=w(endwin )) && w(begwin )>0,
            maxpos = begwin; maxim = w(maxpos );
        elseif (w(endwin )>=w(begwin )) && w(endwin )>0,
            maxpos = endwin; maxim = w(maxpos );
        end
    end
    if isempty(minpos),               % if no local negative minimum
        % the minimum will be the first
        % or the last sample
        if (w(begwin )<=w(endwin )) && w(begwin )<0,
            minpos = begwin; minim = w(minpos );
        elseif (w(endwin )<=w(begwin )) && w(endwin )<0,
            minpos = endwin; minim = w(minpos );
        end
    end

    absmax = abs(maxim);
    absmin = abs(minim);

    if i<length(timenew),
        veficaz = sqrt(mean(w .^2)); % all w 03.Dec.04 interval was restricted before
    else                        % if last beat of the segment
        veficaz = sqrt(nanmean(w.^2)); % all w 03.Dec.04 interval was restricted before
    end

    hay_onda = ((absmax>umbraldetT*veficaz)|(absmin>umbraldetT*veficaz));

    % Rute 18/06/02
    if endwin-begwin<(min_vent*freq) % se a janela tem amplitude menor do que min_vent seg enato nao existe onda T
        hay_onda =0;
    end

    % Is there a wave?
    if hay_onda,
        if absmax >= absmin,    % the greatest modulus maximum is the maximum
            % Now we search the two minima nearest to maxpos, one before and one after
            minapos = max(modmax(w(begwin+1:maxpos-1 ),2,0,-1));
            minapos = begwin + minapos; % Position of the negative minimum before the maximum
            if isempty(minapos) && (maxpos ~= begwin) && (w(begwin )<0),
                minapos = begwin;% If no local minimum before the maximum, take the first sample
            end
            minppos = min(modmax(w(maxpos+1:endwin ),2,0,-1));
            minppos = maxpos + minppos; % Position of the positive maximum after the minimum
            if isempty(minppos) && (maxpos ~= endwin) && (w(endwin )<0),
                minppos = endwin;        % If no local minimum after the maximum, take the last sample
            end

            mina = abs(w(minapos ));     % Amplitude of minimum before maximum
            minp = abs(w(minppos ));     % Amplitude of minimum after maximum

            if (mina < umbralsig*absmax), % If mina is not big enough
                mina =[];                  % forget it
            elseif (maxpos-minapos>Tmax_Tmin_time_min*freq), %or if ther are more than 150 ms to maxpos
                mina = [];
            end
            if (minp < umbralsig*absmax), % If minp is not big enough
                minp =[];                  % forget it
            elseif (minppos-maxpos>Tmax_Tmin_time_min*freq), % and also if there are more than 150 ms to maxpos
                minp =[];
            end

            if ~isnan(mina)&~isnan(minp),      %#ok<AND2> %%% NUEVO JP
                if (mina >= minp)&&(minp < umbralsig*absmax*Tmax_Tmin_bifasic),
                    minp = [];
                elseif (minp> mina)&&(mina < umbralsig*absmax*Tmax_Tmin_bifasic),
                    mina = [];
                end
            end

            % Test which modulus maxima are significative and find zero crossings
            if isempty(mina),
                if isempty(minp),
                    tipoT = 2;      % only upwards T wave
                    if maxpos - minapos > 2,         % if not !!!!!!!!!!!
                        ind = zerocros(flipud(w(minapos:maxpos )));   %Scale 3 !!! % also in scale 4!!!!!!! 03.Dec.04
                        T = maxpos - ind +1;                           % Zero crossing = T wave position
                        picoff = maxpos;                               %wavelet  peak to detect offset
                    elseif isempty(minapos);  % If there were no minimum, there is no zero crossing
                        T = picant (w(begwin:maxpos ),maxpos);       % Take the minimum at scale 4
                        if ~isempty(T) %%%%%%%%%%%% 14/06/02 Rute
                            picoff = maxpos;
                        else
                            picoff = []; % if did not exist a peak in scale 4 there is no T
                        end %%%%%%%%%%%% 14/06/02 Rute
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                else      % minp exists (is significative) but mina not
                    tipoT = 0;  		%normal T wave
                    if minppos -maxpos >2,       % if not!!!!!!???
                        ind = zerocros(w(maxpos:minppos ));  %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind) %%%%%%%%%%%%%%Rute 03/09/02
                        %                             ind = zerocros(w(maxpos:minppos));
                        %                         end %%%%%%%%%%%%%%Rute 03/09/02
                        T = maxpos + ind -1;
                        picon = maxpos;		% For determining onset and offset
                        picoff = minppos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                end
            else
                if isempty(minp),   %  mina exists (is significative) but minp not
                    tipoT = 1;  	%inverted T wave
                    if maxpos -minapos >2,  % if not !!!!!!!!!
                        ind = zerocros(w(minapos:maxpos));  % wavelet zero crossing is T wave peak  %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind)%%%%%%%%%%%%%%Rute 03/09/02
                        %                             ind = zerocros(w(minapos:maxpos ));  % wavelet zero crossing is T wave peak  %% 03/09/02 Rute
                        %                         end %%%%%%%%%%%%%%Rute 03/09/02
                        T = minapos + ind -1;
                        picon = minapos;
                        picoff = maxpos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                else    % both mina and minp are significative.  Biphasic wave.
                    tipoT = 5;	% biphasic neg-pos T wave
                    if maxpos - minapos > 2,    %!!!!!!!!!!!!
                        ind = zerocros(flipud(w(minapos:maxpos ))); %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind) %%%%%%%%%%%%%%Rute 03/09/02
                        %                             ind = zerocros(flipud(w(minapos:maxpos ))); %% 03/09/02 Rute
                        %                         end%%%%%%%%%%%%%%Rute 03/09/02
                        T = maxpos - ind +1;
                        picon = minapos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                    if minppos - maxpos > 2,    %!!!!!!!!!!!!
                        ind = zerocros(flipud(w(maxpos:minppos ))); %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind)%%%%%%%%%%%%%%Rute 03/09/02
                        %                             ind = zerocros(flipud(w(maxpos:minppos ))); %% 03/09/02 Rute
                        %                         end %%%%%%%%%%%%%%Rute 03/09/02
                        %Tprima = maxpos + ind -1; %18.11.05
                        Tprima = minppos- ind +1;
                        picoff = minppos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                end
            end
        else        % If the greatest modulus maximum is the minimum
            % Search two maxima, one before and one after the minimum
            maxapos = max(modmax(w(begwin+1:minpos-1 ),2,0,1));
            maxapos = begwin + maxapos;
            if isempty(maxapos) && (minpos ~= begwin) && (w(begwin )>0),
                maxapos = begwin;
            end
            maxppos = min(modmax(w(minpos+1:endwin ),2,0,1));
            maxppos = minpos + maxppos;
            if isempty(maxppos) && (minpos ~= endwin) && (w(endwin )>0),
                maxppos = endwin;
            end
            maxa = abs(w(maxapos ));
            maxp = abs(w(maxppos )) ;                                   % See if they are significative
            if (maxa < umbralsig*absmin)
                maxa =[];
            elseif (minpos-maxapos>Tmax_Tmin_time_min*freq),
                maxa = [];
            end
            if (maxp < umbralsig*absmin),
                maxp =[];
            elseif (maxppos-minpos>Tmax_Tmin_time_min*freq),
                maxp = [];
            end
            if ~isnan(maxa)&~isnan(maxp),      %#ok<AND2> %%% NUEVO JP
                if (maxa >= maxp)&&(maxp < umbralsig*absmin*Tmax_Tmin_bifasic),
                    maxp = [];
                elseif (maxp> maxa)&&(maxa < umbralsig*absmin*Tmax_Tmin_bifasic),
                    maxa = [];
                end
            end
            % Test which modulus maxima are significative and find zero crossings
            if isempty(maxa),
                if isempty(maxp),
                    tipoT = 3;      % only downwards T wave
                    if minpos - maxapos > 2, %!!!!!!!!
                        ind = zerocros(flipud(w(maxapos:minpos )));   %Scale 3 % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind) %%%%%%%%%%%%%%Rute 03/09/02
                        %                             ind = zerocros(flipud(w(maxapos:minpos )));
                        %                         end %%%%%%%%%%%%%%Rute 03/09/02
                        T = minpos - ind +1;
                        picoff = minpos;
                    elseif isempty(maxapos);  % If there were no maximum, there is no zero crossing.
                        T = picant (w(begwin:minpos ),minpos);  % menimo en escala 4.
                        if ~isempty(T) %%%%%%%%%%%% 14/06/02 Rute
                            picoff = minpos;
                        else
                            picoff = []; % if did not exist a peak in scale 4 there is no T
                        end %%%%%%%%%%%% 14/06/02 Rute
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                else      % maxp is signficative, but not maxa
                    tipoT = 1;  %inverted T wave
                    if maxppos -minpos >2,  % !!!!!!
                        ind = zerocros(w(minpos:maxppos ));  %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind) %%%%%%%%%%%%%%Rute 03/09/02
                        %                             ind = zerocros(w(minpos:maxppos ));  %% 03/09/02 Rute
                        %                         end %%%%%%%%%%%%%%Rute 03/09/02
                        T = minpos + ind -1;
                        picon = minpos;		% For calculating onset and offset
                        picoff = maxppos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                end
            else
                if isempty(maxp),   % maxa is significative, but not maxp
                    tipoT = 0;       %normal T wave
                    if minpos -maxapos >2,
                        ind = zerocros(w(maxapos:minpos )); %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind) %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        %                             ind = zerocros(w(maxapos:minpos )); %% 03/09/02 Rute
                        %end %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        T = maxapos + ind -1;
                        picon = maxapos;
                        picoff = minpos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                else    % both maxa and maxp are significative.  Biphasic wave.
                    tipoT = 4;	% biphasic pos-neg T wave
                    if minpos - maxapos > 2,  %!!!!!!!!!!!
                        ind = zerocros(flipud(w(maxapos:minpos ))); %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind) %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        %                             ind = zerocros(flipud(w(maxapos:minpos ))); %% 03/09/02 Rute
                        %                         end %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        T = minpos - ind +1;
                        picon = maxapos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        %                         [i tipoT]
                        disp('caso nao previsto; onda T nao marcada!')
                    end
                    if maxppos - minpos > 2,    %!!!!!!!!!!!!
                        ind = zerocros(flipud(w(minpos:maxppos ))); %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind) %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        %                             ind = zerocros(flipud(w(minapos:maxpos ))); %% 03/09/02 Rute
                        %                         end %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        %Tprima = minpos + ind -1; %18.11.05 DUVIDA
                        Tprima = maxppos- ind +1;
                        picoff = maxppos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                end
            end
        end

        % T wave onset and offset detection
        %if isempty(T), picon=[]; picoff=[]; end

        picon_keep=[picon_keep picon]; %multileadchange
        picoff_keep=[picoff_keep picoff];%multileadchange
        lead_keep=[lead_keep 4];%multileadchange

        if ~isempty(picon),
            Ton = searchon (picon, w(max(begwin,picon-0.12*freq):picon ), Kton);
        end
        if ~isempty(picoff),
            Toff=searchoff(picoff, w(picoff:min([size(w,1) picoff+0.12*freq ])  ) , Ktoff);
            if (Toff > endwin),
                Toff = endwin;
            end
        end
    else
        tipoT=9;
        messages.warnings=[messages.warnings {'unknown case: unable to find T wave in multilead approach using scale 4.'}];
    end       %% utilizar a escala 5 07/06/02 Rute %%%%%%%%%%%%%%%%%%%%%%%%%%%

end
% Filling the structure with positions
if isempty(Ton), Ton = NaN; end;
if isempty(Toff), Toff = NaN;  end;
if isempty(T), T=NaN; end;
if isnan(T),
    picon_keep=[picon_keep NaN]; %multileadchange
    picoff_keep=[picoff_keep NaN];%multileadchange
    lead_keep=[lead_keep NaN];%#ok<NASGU> %multileadchange
end
if isempty(Tprima), Tprima=NaN; end;
if isempty(tipoT), tipoT=NaN; end;
pos.Ton= Ton+samp(1)-1+timenew(i)-1;
pos.Toff= Toff+samp(1)-1+timenew(i)-1;
pos.T= T+samp(1)-1+timenew(i)-1;
pos.Tprima= Tprima+samp(1)-1+timenew(i)-1;
pos.Ttipo= tipoT;
% T = []; Tprima=[]; picon=[]; picoff=[]; Ton=[]; Toff=[]; tipoT=[];
% minapos = []; minppos=[]; maxapos=[]; maxppos=[]; mina=[]; minp=[]; maxa=[];
% maxp=[];
%end
% position.Ton(intervalo(1):intervalo(2)) = pos.Ton+samp(1)-1;
% position.Toff(intervalo(1):intervalo(2))= pos.Toff+samp(1)-1;
% position.T(intervalo(1):intervalo(2)) = pos.T+samp(1)-1;
% position.Tprima(intervalo(1):intervalo(2)) = pos.Tprima+samp(1)-1;
% position.Ttipo(intervalo(1):intervalo(2)) = pos.Ttipo;
% pos.Ton
% Ton = pos.Ton(i)
% Toff= pos.Toff(i)
% T= pos.T(i)
% Tprima = pos.Tprima(i)
% Ttipo= pos.Ttipo(i);
%timenew(i)
%picon_keep
picoff=picoff_keep+timenew(i)-1;
picon=picon_keep+timenew(i)-1;

author: marianux

repo: ecg-kit

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


text:
classdef PTKAirwaysLabelledByLobe < PTKPlugin
    % PTKAirwaysLabelledByLobe. Plugin to visualise the airway tree labelled
    % according to lobe
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     PTKAirwaysLabelledByLobe calls the PTKAirways plugin to segment the
    %     airway tree and the PTKAirwayCentreline plguin to obtain the airway
    %     centreline. It then calls the library function PTKGetAirwaysLabelledByLobe
    %     to allocate each broncus to the unique lobe it serves. Bronchi serving
    %     more than one lobe are not displayed.
    %
    %     The resulting image is a labelled image with nonzero values
    %     representing bronchi allocated to the following lobes:
    %         1 - Upper right lobe
    %         2 - Mid right lobe
    %         4 - Lower right lobe
    %         5 - Upper left lobe
    %         6 - Lower left lobe
    %         3 - Lobe could not be determined with certainty
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Airways <br>by lobe'
        ToolTip = 'Shows airways coloured by lobe, derived from analysis of the airway centreline'
        Category = 'Airways'

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

        EnableModes = MimModes.EditMode
        SubMode = MimSubModes.ColourRemapEditing
        EditRequiresPluginResult = true
    end

    methods (Static)
        function results = RunPlugin(dataset, reporting)
            results_image = dataset.GetTemplateImage(PTKContext.LungROI);
            [airway_results, airway_image] = dataset.GetResult('PTKAirways');
            centreline_results = dataset.GetResult('PTKAirwayCentreline');

            % Generate a map of all the bronchi (one mapped colour for each bronchus)
            [airway_mapped_image, labelled_centreline] = PTKMapAirwayCentrelineToImage(centreline_results, airway_image);

            % Find the airways for each lobe
            start_branches = PTKGetAirwaysLabelledByLobe(results_image, labelled_centreline, reporting);

            airway_mapping = PTKAirwaysLabelledByLobe.GetAirwayMappingForLobes(start_branches);
            airway_mapped_image.ChangeColorLabelMap(airway_mapping);

            results = [];
            results.AirwaysByLobeImage = airway_mapped_image;
            results.StartBranches = start_branches;
        end

        function start_branches = GetLobarStartBranchesFromAirwayMap(current_start_branches, airway_mapped_image)

            tree_root = current_start_branches.Trachea;
            left_upper_lobe = [];
            left_lower_lobe = [];
            right_upper_lobe = [];
            right_mid_lobe = [];
            right_lower_lobe = [];
            left_uncertain = [];

            % Generate the mapping for the lobes
            bronchi_to_do = CoreStack(tree_root);
            while ~bronchi_to_do.IsEmpty
                next_bronchus = bronchi_to_do.Pop;
                bronchus_label = next_bronchus.BronchusIndex;
                bronchus_mapping = airway_mapped_image.ColorLabelMap(bronchus_label + 1);
                if bronchus_mapping == 7
                    bronchi_to_do.Push(next_bronchus.Children);
                else
                    switch bronchus_mapping
                        case PTKColormapLabels.LeftUpperLobe
                            left_upper_lobe = [left_upper_lobe, next_bronchus];

                        case PTKColormapLabels.LeftLowerLobe
                            left_lower_lobe = [left_lower_lobe, next_bronchus];

                        case PTKColormapLabels.RightUpperLobe
                            right_upper_lobe = [right_upper_lobe, next_bronchus];

                        case PTKColormapLabels.RightLowerLobe
                            right_lower_lobe = [right_lower_lobe, next_bronchus];

                        case PTKColormapLabels.RightMiddleLobe
                            right_mid_lobe = [right_mid_lobe, next_bronchus];

                        case 3
                            left_uncertain = [left_uncertain, next_bronchus];
                    end
                end
            end

            start_branches = [];
            start_branches.Trachea = current_start_branches.Trachea;
            start_branches.Left = current_start_branches.Left;
            start_branches.Right = current_start_branches.Right;
            start_branches.LeftLower = left_lower_lobe;
            start_branches.LeftUpper = left_upper_lobe;
            start_branches.RightUpper = right_upper_lobe;
            start_branches.RightLower = right_lower_lobe;
            start_branches.RightMid = right_mid_lobe;
            start_branches.LeftUncertain = left_uncertain;
        end

        function airway_mapping = GetAirwayMappingForLobes(start_branches)

            airway_mapping = [];

            % Map the whole tree initially to 0 to take into account any
            % parent parenches of the lobar bronchi
            airway_mapping = PTKAirwaysLabelledByLobe.MapTheseBranchesToLabel(airway_mapping, start_branches.Trachea, 0);

            % Map the labels to the lobar colours
            airway_mapping = PTKAirwaysLabelledByLobe.MapTheseBranchesToLabel(airway_mapping, start_branches.LeftUpper, PTKColormapLabels.LeftUpperLobe);
            airway_mapping = PTKAirwaysLabelledByLobe.MapTheseBranchesToLabel(airway_mapping, start_branches.LeftLower, PTKColormapLabels.LeftLowerLobe);
            airway_mapping = PTKAirwaysLabelledByLobe.MapTheseBranchesToLabel(airway_mapping, start_branches.RightUpper, PTKColormapLabels.RightUpperLobe);
            airway_mapping = PTKAirwaysLabelledByLobe.MapTheseBranchesToLabel(airway_mapping, start_branches.RightLower, PTKColormapLabels.RightLowerLobe);
            airway_mapping = PTKAirwaysLabelledByLobe.MapTheseBranchesToLabel(airway_mapping, start_branches.RightMid, PTKColormapLabels.RightMiddleLobe);
            airway_mapping = PTKAirwaysLabelledByLobe.MapTheseBranchesToLabel(airway_mapping, start_branches.LeftUncertain, 3);
            airway_mapping = PTKAirwaysLabelledByLobe.MapTheseBranchesToLabel(airway_mapping, start_branches.RightUncertain, 3);
            airway_mapping(airway_mapping == 0) = 7;
            airway_mapping(1) = 0;
        end

        function airway_mapping = MapTheseBranchesToLabel(airway_mapping, branch_list, label)
            branches_to_label = CoreStack(branch_list);
            while ~branches_to_label.IsEmpty
                branch = branches_to_label.Pop;
                airway_mapping(branch.BronchusIndex + 1) = label;
                branches_to_label.Push(branch.Children);
            end
        end

        function results = GenerateImageFromResults(airway_results, ~, ~)
            results = airway_results.AirwaysByLobeImage;
        end

        function results = GetEditedResult(saved_result, edited_result_image, ~)
            start_branches = PTKAirwaysLabelledByLobe.GetLobarStartBranchesFromAirwayMap(saved_result.StartBranches, edited_result_image);

            results = [];
            results.AirwaysByLobeImage = edited_result_image;
            results.StartBranches = start_branches;
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%% Initvvx
%%
%%

a = findobj('style','edit');
for i = 1:length(a),
   eval( [get(a(i),'tag') '=' get(a(i),'string') ';'] );
end

rho = 1000;
useoutput = get(findobj('tag','useoutput'),'value');
cocurrent = get(findobj('tag','cocurrent'),'value');


text:
function quiverSection(sVF1,sVF2,N,varargin)
% plot a vector field along another vectorfield
%
% Syntax
%
%   N = vector3d.Z;
%   quiverSection(sVF1,v,N)
%   quiverSection(sVF1,sVF2,N,pi/3)
%
% Input
%  sVF1, sVF2 - @S2VectorField
%  v - @vector3d
%  N - normal @vector3d of the section
%
% Options
%  normalized - draw unit length vectors

[mtexFig,isNew] = newMtexFigure(varargin{:});

% where to plot - compute circle positions
omega = linspace(0,2*pi,36);

if nargin > 2 && isnumeric(varargin{1})
  eta = varargin{1};
else
  eta = pi/2;
end

S1 = axis2quat(N,omega)*axis2quat(orth(N),eta)*N;
circ = reshape(sVF1.eval(S1),length(S1), []);

% what to plot
if isa(sVF2,'function_handle')
  v = sVF2(S1);
else
  v = sVF2.eval(S1);
end
v = v(:);
if check_option(varargin,'normalized'), v = v.normalize; end

% plot the vector field v at the positions circ
if v.antipodal
  opt = {'showArrowHead','off'};
  h = quiver3(circ.x,circ.y,circ.z,-v.x,-v.y,-v.z,'parent',mtexFig.gca,opt{:});
  set(get(get(h,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
else
  opt = {};
  h = [];
end

h = [h,quiver3(circ.x,circ.y,circ.z,v.x,v.y,v.z,'parent',mtexFig.gca,opt{:})];

% post process output
view(mtexFig.gca,squeeze(double(N)));
set(mtexFig.gca,'dataAspectRatio',[1 1 1]);
optiondraw(h,varargin{:});

if isNew, mtexFig.drawNow('figSize',getMTEXpref('figSize'),varargin{:}); end


text:
function varargout=lpcar2im(varargin)
%
% For calling details please see v_lpcar2im.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpcar2im.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpcar2im.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpcar2im(varargin{:});
else
    v_lpcar2im(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function coeffS = computeMcoeffs(coeffS,symmetric)
% function coeffS = computeMcoeffs(coeffS,symmetric)
%
% APA, 6/11/2018

if ~exist('symmetric','var')
    symmetric = true;
end

%% ComputeRemainingCoefficients(bool symmetric)

if symmetric

    coeffS.M1 = coeffS.N1 - coeffS.D1 * coeffS.N0;
    coeffS.M2 = coeffS.N2 - coeffS.D2 * coeffS.N0;
    coeffS.M3 = coeffS.N3 - coeffS.D3 * coeffS.N0;
    coeffS.M4 = -coeffS.D4 * coeffS.N0;

else

   coeffS.M1 = -( coeffS.N1 - coeffS.D1 * coeffS.N0 );
   coeffS.M2 = -( coeffS.N2 - coeffS.D2 * coeffS.N0 );
   coeffS.M3 = -( coeffS.N3 - coeffS.D3 * coeffS.N0 );
   coeffS.M4 = coeffS.D4 * coeffS.N0;
end

% Compute coefficients to be used at the boundaries ...
% in order to simulate edge extension boundary conditions.
coeffS.SN = coeffS.N0 + coeffS.N1 + coeffS.N2 + coeffS.N3;
coeffS.SM = coeffS.M1 + coeffS.M2 + coeffS.M3 + coeffS.M4;
coeffS.SD = 1.0 + coeffS.D1 + coeffS.D2 + coeffS.D3 + coeffS.D4;

coeffS.BN1 = coeffS.D1 * coeffS.SN / coeffS.SD;
coeffS.BN2 = coeffS.D2 * coeffS.SN / coeffS.SD;
coeffS.BN3 = coeffS.D3 * coeffS.SN / coeffS.SD;
coeffS.BN4 = coeffS.D4 * coeffS.SN / coeffS.SD;

coeffS.BM1 = coeffS.D1 * coeffS.SM / coeffS.SD;
coeffS.BM2 = coeffS.D2 * coeffS.SM / coeffS.SD;
coeffS.BM3 = coeffS.D3 * coeffS.SM / coeffS.SD;
coeffS.BM4 = coeffS.D4 * coeffS.SM / coeffS.SD;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function pth = ffpath(fname)
%   FFPATH    Find file path
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The function browses very fast current directory and directories known in
% 'matlabpath' and the system variable 'path'. It searches for the file,
% name of which is in the input argument 'fname'. If a directory is found,
% the output argument pth is filled by path to the file name from the input
% argument, otherwise pth is empty.
% File names should have their extensions, but MATLAB m-files.
%
% Arguments:
%   fname   file name
%   pth     path to the fname
%
% Examples:
%   pth = ffpath('gswin32c.exe')
%   pth =
%   c:\Program Files\gs\gs8.60\bin\
%
%   pth = ffpath('hgrc')
%   pth =
%   C:\PROGRA~1\MATLAB\R2006b\toolbox\local

% Miroslav Balda
% miroslav AT balda DOT cz
% 2008-12-15    v 0.1   only for system variable 'path'
% 2008-12-20    v 1.0   for both 'path' and 'matlabpath'

if nargin<1
    error('The function requires one input argument (file name)')
end
pth = filesepStandard(pwd);
if exist([pth '/' fname],'file')
    return
end % fname found in current dir

tp = matlabpath;
t  = 0;
if isunix() | ismac()
    I = [t, findstr(tp,':'), length(tp)+1];
elseif ispc()
    I = [t, findstr(tp,';'), length(tp)+1];
end
for k = 1:length(I)-1               %   search in path's directories
    pth = tp(I(k)+1:I(k+1)-1);
    % fprintf('%s\n', [pth '/' fname]);
    if exist([pth '/' fname],'file')
        return;
    end
end
[status,tp] = system('path');
t = 5;
pth = '';

author: BUNPC


text:
function out=fix(x)

precision=x(1).precision;
out_rval=cell(size(x));
out_ival=cell(size(x));

for ii=1:numel(x)
 imag=false;
 [xrval,xival]=getVals(x,ii);
 if hasimag(xival), imag=true; end
 if imag
  out_rval{ii}=mpfr_fix(precision,xrval);
  out_ival{ii}=mpfr_fix(precision,xival);
 else
  out_rval{ii}=mpfr_fix(precision,xrval);
 end
end % for ii=1:max(ex,
out=class(struct('rval',out_rval,...
                  'ival',out_ival,...
                  'precision',precision),'mp');

author: opencobra

repo: cobratoolbox


text:
function voronoi_neighbor_test01 ( )

%*****************************************************************************80
%
%% VORONOI_NEIGHBOR_TEST01 demonstrates VORONOI_NEIGHBORS
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'VORONOI_NEIGHBORS_TEST01:\n' );
  fprintf ( 1, '  Select a random set of points in the unit square.\n' );
  fprintf ( 1, '  Compute the Voronoi diagram.\n' );
  fprintf ( 1, '  Have VORONOI_NEIGHBORS determine the neighbors.\n' );

  n = 8;
  x = rand ( n, 2 );
%
%  Compute and display the Voronoi diagram.
%
  figure ( );

  hold on
  voronoi ( x(:,1), x(:,2) );
  for i = 1 : n
    txt = sprintf ( '%d', i );
    text ( x(i,1), x(i,2), txt, 'FontSize', 12 );
  end
  hold off

  filename = 'voronoi_neighbors_test01.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Voronoi diagram saved as "%s".\n', filename );
%
%  Compute and list the finite Voronoi edges.
%
  [ V, C ] = voronoin ( x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Voronoi edges:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : length ( C )
    disp ( C{i} )
  end
%
%  Compute the Voronoi neighbors.
%
  vn = voronoi_neighbors ( x );
%
%  To print the matrix, make a full version.
%
  vn = full ( vn );
  i4mat_print ( n, n, vn, '  Voronoi adjacency:' )

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ s, t ] = segment_point_coords_2d ( p1, p2, p )

%*****************************************************************************80
%
%% SEGMENT_POINT_COORDS_3D: coordinates of a point on a line segment in 3D.
%
%  Discussion:
%
%    A line segment is the finite portion of a line that lies between
%    two points P1 and P2.
%
%    By the coordinates of a point P with respect to a line segment [P1,P2]
%    we mean numbers S and T such that S gives us the distance from the
%    point P to the nearest point PN on the line (not the line segment!),
%    and T gives us the position of PN relative to P1 and P2.
%
%    If S is zero, then P lies on the line.
%
%    If 0 <= T <= 1, then PN lies on the line segment.
%
%    If both conditions hold, then P lies on the line segment.
%
%    If E is the length of the line segment, then the distance of the
%    point to the line segment is:
%
%      sqrt ( S**2 +  T**2    * E**2 )     if      T <= 0;
%             S                            if 0 <= T <= 1
%      sqrt ( S**2 + (T-1)**2 * E**2 )     if 1 <= T
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 August 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(3), P2(3), the endpoints of the line segment.
%
%    Input, real P(3), the point to be considered.
%
%    Output, real S, the distance of P to the nearest point PN
%    on the line through P1 and P2.  (S will always be nonnegative.)
%
%    Output, real T, the relative position of the point PN
%    to the points P1 and P2.
%
  dim_num = 3;
%
%  If the line segment is actually a point, then the answer is easy.
%
  if ( p1(1:dim_num) == p2(1:dim_num) )

    t = 0.0;

  else

    bot = sum ( ( p2(1:dim_num) - p1(1:dim_num) ).^2 );

    t = ( p(1:dim_num)  - p1(1:dim_num) ) ...
      * ( p2(1:dim_num) - p1(1:dim_num) )' / bot;

  end

  pn(1:dim_num) = p1(1:dim_num) + t * ( p2(1:dim_num) - p1(1:dim_num) );

  s = sqrt ( sum ( ( p(1:dim_num) - pn(1:dim_num) ).^2 ) );

  return
end


text:
function ir = sofa_get_data_fir(sofa,idx)
%SOFA_GET_DATA_FIR returns a FIR data matrix from a SOFA file or struct
%
%   Usage: ir = sofa_get_data_fir(sofa,[idx])
%
%   Input parameters:
%       sofa    - impulse response data set (SOFA struct/file)
%       idx     - index of the single impulse responses that should be returned.
%                 idx could be a single value, then only one impulse response
%                 will be returned, or it can be a vector then all impulse
%                 responses for the corresponding index positions will be
%                 returned (default: return all impulse responses)
%
%   Output parameters:
%       ir      - impulse response (M,2,N), where
%                   M ... number of impulse responses
%                   N ... samples
%
%   SOFA_GET_DATA_FIR(sofa,idx) returns impulse response of the given
%   SOFA file or struct, specified by idx. If no idx is specified all data
%   contained in sofa is returned.
%   For the struct the SOFA file has to loaded before with SOFAload().
%   For a description of the SOFA file format see: https://sofaconventions.org
%
%   See also: sofa_get_data_fire, sofa_get_header, get_ir, SOFAload

%*****************************************************************************
% 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 of input  parameters ==================================
nargmin = 1;
nargmax = 2;
narginchk(nargmin,nargmax)
if nargin==nargmax-1
    idx = [];
else
    isargvector(idx);
end

%% ===== Computation ====================================================
if length(idx)==0
    if sofa_is_file(sofa)
        sofa = SOFAload(sofa);
    end
    ir = sofa.Data.IR;
else
    header = sofa_get_header(sofa);
    if sofa_is_file(sofa)
        ir = zeros(length(idx),2,header.API.N);
        for ii=1:length(idx)
            tmp = SOFAload(sofa,[idx(ii) 1]);
            ir(ii,:,:) = tmp.Data.IR;
        end
    else
        ir = sofa.Data.IR(idx,:,:);
    end
end

author: sfstoolbox

repo: sfs-matlab


text:
function y = minus(X,Y)
%MINUS (overloaded)

global FACTORTRACKING

% Cannot use isa here since blkvar is marked as sdpvar
X_class = class(X);
Y_class = class(Y);
X_is_spdvar = strcmp(X_class,'sdpvar');
Y_is_spdvar = strcmp(Y_class,'sdpvar');

% Convert block objects
if ~X_is_spdvar && ~strcmp(X_class,'double')
    if isa(X,'blkvar')
        X = sdpvar(X);
        X_is_spdvar = isa(X,'sdpvar');
    elseif isa(X,'intval')
        X_is_spdvar = 0;
        Y.basis = intval(Y.basis);
    elseif isa(X,'uint8') || isa(X,'uint16') || isa(X,'uint32') || isa(X,'uint64')
        X = double(X);
    end
end

if ~Y_is_spdvar && ~strcmp(Y_class,'double')
    if isa(Y,'blkvar')
        Y = sdpvar(Y);
        Y_is_spdvar = isa(Y,'sdpvar');
    elseif isa(Y,'intval')
        Y_is_spdvar = 0;
        X.basis = intval(X.basis);
    elseif isa(Y,'uint8') || isa(Y,'uint16') || isa(Y,'uint32') || isa(Y,'uint64')
        Y = double(Y);
    end
end

if X_is_spdvar
    if X.typeflag == 40
        y = X + uminus(Y);
        return
    end
else
    if any(isnan(X))
        error('Adding NaN to an SDPVAR makes no sense.');
    end
end
if Y_is_spdvar
    if Y.typeflag == 40
        y =X + uminus(Y);
        return
    end
else
    if any(isnan(Y))
        error('Adding NaN to an SDPVAR makes no sense.');
    end
end

switch 2*X_is_spdvar+Y_is_spdvar
    case 1
        if isempty(X)
            try
                y = full(X - reshape(Y.basis(:,1),Y.dim(1),Y.dim(2)));
            catch
                error(lasterr);
            end
            return
        end

        y = Y;
        n_Y = Y.dim(1);
        m_Y = Y.dim(2);
        [n_X,m_X] = size(X);
        x_isscalar = (n_X*m_X==1);
        y_isscalar = (n_Y*m_Y==1);
        any_scalar = x_isscalar | y_isscalar;

        % Speeeeeeed
        if x_isscalar && y_isscalar
            y.basis = -y.basis;
            tmp = y.basis(1)+X;
            if (isequal(class(tmp),'gem') || isequal(class(tmp),'sgem')) && ~isequal(class(y.basis), class(tmp))
                y.basis = gemify(y.basis);
            end
            y.basis(1) = tmp;
            % Reset info about conic terms
            y.conicinfo = [0 0];
            y.extra.opname='';
            if FACTORTRACKING, y = addfactors(y,X,-Y);end
            return
        end

        if any_scalar || all([n_Y m_Y]==[n_X m_X])
            if y_isscalar
                y.basis = repmat(y.basis,n_X*m_X,1);
                y.dim(1) = n_X;
                y.dim(2) = m_X;
            end
            y.basis = -y.basis;
            if nnz(X)~=0
                tmp = y.basis(:,1)+X(:);
                if (isequal(class(tmp),'gem') || isequal(class(tmp),'sgem')) && ~isequal(class(y.basis), class(tmp))
                    y.basis = gemify(y.basis);
                end
                y.basis(:,1) = tmp;
            end
        else
            error('Matrix dimensions must agree.');
        end
        % Reset info about conic terms
        y.conicinfo = [0 0];
        y.extra.opname='';
        y.extra.createTime = definecreationtime;
        if FACTORTRACKING, y = addfactors(y,X,-Y);end
    case 2

        if isempty(Y)
            try
                y = full(reshape(X.basis(:,1),X.dim(1),X.dim(2))-Y);
            catch
                error(lasterr);
            end
            return
        end
        y = X;
        n_X = X.dim(1);
        m_X = X.dim(2);
        [n_Y,m_Y] = size(Y);
        x_isscalar = (n_X*m_X==1);
        y_isscalar = (n_Y*m_Y==1);
        any_scalar = x_isscalar | y_isscalar;

        % Silly hack
        % Taking X-scalar(0) takes unnecessary time
        % and is used in most definitions of LMIs
        if (y_isscalar && (Y==0))
            return
        end

        % Speeeeeeed
        if x_isscalar && y_isscalar
            tmp = y.basis(1)-Y;
            if (isequal(class(tmp),'gem') || isequal(class(tmp),'sgem')) && ~isequal(class(y.basis), class(tmp))
                y.basis = gemify(y.basis);
            end
            y.basis(1) = tmp;
            % Reset info about conic terms
            y.conicinfo = [0 0];
            y.extra.opname='';
            if FACTORTRACKING, y = addfactors(y,X,-Y);end
            return
        end

        if any_scalar || all([n_Y m_Y]==[n_X m_X])
            if x_isscalar
                y.basis = repmat(y.basis,n_Y*m_Y,1);
                y.dim(1) = n_Y;
                y.dim(2) = m_Y;
            end
            tmp = y.basis(:,1)-Y(:);
            if (isequal(class(tmp),'gem') || isequal(class(tmp),'sgem')) && ~isequal(class(y.basis), class(tmp))
                y.basis = gemify(y.basis);
            end
            y.basis(:,1) = tmp;
        else
            error('Matrix dimensions must agree.');
        end

        % Update information about conic terms
        % This information is used in DUALIZE to
        % speed up some checks, and to facilitate some
        % advanced dualization features. It also
        % speeds up checking for symmetry in some other code
        % Ugly, but the best way at the moment
        % For a description of this field, check SDPVAR code
        % if (y.conicinfo(1)~=0) && isequal(Y,Y') && (y.conicinfo(2) ~= 2)
        %     y.conicinfo(2) = max(1,y.conicinfo(2));
        % else
        y.conicinfo = [0 0];
        y.extra.opname='';
        if FACTORTRACKING, y = addfactors(y,X,-Y);end
        % end

    case 3

        %	if (X.typeflag~=0) || (Y.typeflag~=0)
        %		error('Relational objects cannot be manipulated')
        %	end

        n_X = X.dim(1);
        m_X = X.dim(2);
        n_Y = Y.dim(1);
        m_Y = Y.dim(2);
        x_isscalar = (n_X*m_X==1);
        y_isscalar = (n_Y*m_Y==1);
        any_scalar = x_isscalar | y_isscalar;

        if ~any_scalar
            if (~((n_X==n_Y) && (m_X==m_Y)))
                error('Matrix dimensions must agree.')
            end
        end

        Xlmi_variables = X.lmi_variables;
        Ylmi_variables = Y.lmi_variables;
        yFirst = 0;
        xFirst = 0;
        if Ylmi_variables(end) < Xlmi_variables(1)
            all_lmi_variables = [Ylmi_variables Xlmi_variables];
            yFirst = 1;
        elseif Xlmi_variables(end) < Ylmi_variables(1)
            all_lmi_variables = [Xlmi_variables Ylmi_variables];
            xFirst = 1;
        else
            all_lmi_variables = uniquestripped([Xlmi_variables Ylmi_variables]);
        end

        y = X;
        %X.basis = []; % Returns memory?
        y.lmi_variables = all_lmi_variables;

        if isequal(all_lmi_variables,Xlmi_variables)
            in_X_logical = ones(1,length(Xlmi_variables));
            in_X = 1:length(Xlmi_variables);
        else
            in_X_logical = ismembcYALMIP(all_lmi_variables,Xlmi_variables);
            in_X = find(in_X_logical);
        end

        if isequal(all_lmi_variables,Ylmi_variables)
            in_Y_logical = ones(1,length(Ylmi_variables));
            in_Y = 1:length(Ylmi_variables);
        else
            in_Y_logical = ismembcYALMIP(all_lmi_variables,Ylmi_variables);
            in_Y = find(in_Y_logical);
        end

        if isequal(Xlmi_variables,Ylmi_variables) && n_Y==n_X && m_Y==m_X
            y.basis = y.basis - Y.basis;
            % Super special case f(scalar)-f(scalar)
            if length(X.lmi_variables)==1
                if all(y.basis(:,2)==0)
                    y = full(reshape(y.basis(:,1),n_Y,m_Y));
                else
                    y.conicinfo = [0 0];
                    y.extra.opname='';
                    if FACTORTRACKING, y = addfactors(y,X,-Y);end
                end
                return
            end
        elseif max(Xlmi_variables) < min(Ylmi_variables) &&  n_Y==n_X && m_Y==m_X
            % Disjoint variables in X - Y
            y.basis = [y.basis(:,1) - Y.basis(:,1) y.basis(:,2:end) -Y.basis(:,2:end)];
        elseif max(Ylmi_variables) < min(Xlmi_variables) &&  n_Y==n_X && m_Y==m_X
            % Disjoint variables in X - Y
            y.basis = [y.basis(:,1) - Y.basis(:,1) -Y.basis(:,2:end) y.basis(:,2:end)];
        else
           % [ix,jx,sx] = find(y.basis);y.basis = [];
           % [iy,jy,sy] = find(Y.basis);%Y.basis = [];
            mapX = [1 1+in_X];
            mapY = [1 1+in_Y];
           % basis_X = sparse(ix,mapX(jx),sx,n_X*m_X,1+length(all_lmi_variables));ix=[];jx=[];sx=[];
           % basis_Y = sparse(iy,mapY(jy),sy,n_Y*m_Y,1+length(all_lmi_variables));iy=[];jy=[];sy=[];
            basis_X = X.basis*(sparse(1:length(mapX),mapX,1,size(X.basis,2),length(all_lmi_variables)+1));
            basis_Y = Y.basis*(sparse(1:length(mapY),mapY,1,size(Y.basis,2),length(all_lmi_variables)+1));

            % Fix addition of matrix+scalar
            if n_X*m_X<n_Y*m_Y
                y.dim(1) = n_Y;
                y.dim(2) = m_Y;
                try
                    y.basis = bsxfun(@plus,basis_X,-basis_Y);
                catch
                    basis_X = repmat(basis_X,n_Y*m_Y,1);
                    y.basis = basis_X - basis_Y;
                end

            elseif n_Y*m_Y<n_X*m_X
                y.dim(1) = n_X;
                y.dim(2) = m_X;
                try
                     y.basis = bsxfun(@minus,basis_X,basis_Y);
                catch
                    basis_Y = repmat(basis_Y,n_X*m_X,1);
                    y.basis = basis_X - basis_Y;
                end
            else
                      y.basis = basis_X - basis_Y;
            end
        end

        y.conicinfo = [0 0];
        y.extra.opname='';
        y.extra.createTime = definecreationtime;
        if FACTORTRACKING, y = addfactors(y,X,-Y);end
        if nnz(in_Y_logical & in_X_logical)>0
            y = clean(y);
        end

    otherwise
end

% Update info on KYP objects
if X_is_spdvar && Y_is_spdvar
  if  X.typeflag==9  && Y.typeflag==9
    error('Substraction of KYP objects currently not supported')
  end
end
if Y_is_spdvar
  if  Y.typeflag==9
    y.extra.M = -Y.extra.M+X;
    y.extra.negated = ~Y.extra.negated;
    return
  end
end
if X_is_spdvar
 if X.typeflag==9
    y.extra.M = y.extra.M-Y;
    return
 end
end

author: shamilmamedov


text:
% Mathematics Q3079400
% https://math.stackexchange.com/questions/3079400
% Numerical Implementation: Solution for the Euler Lagrange Equation Of the Rudin Osher Fatemi (ROF) Total Variation Denoising Model
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     29/03/2019
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0; %<! Continue from Question 1
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Parameters

numRows = 5;
numCols = 4;

%% Load / Generate Data

mX = reshape(1:(numRows * numCols), numRows, numCols);

%% Gradient by MATLAB Notation

vDv = mX(1:end - 1, :) - mX(2:end, :);
vDv = vDv(:);
vDh = mX(:, 1:end - 1) - mX(:, 2:end);
vDh = vDh(:);

vDRef = [vDv; vDh];

%% Gradient by MATLAB Matrix Operation

mD = CreateGradientOperator(numRows, numCols);
vD = mD * mX(:);

%% Analysis

vE = abs(vD - vD);
max(vE)

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function cdf = hypergeometric_cdf ( x, n, m, l )

%*****************************************************************************80
%
%% HYPERGEOMETRIC_CDF evaluates the Hypergeometric CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X, the argument of the CDF.
%
%    Input, integer N, the number of balls selected.
%    0 <= N <= L.
%
%    Input, integer M, the number of white balls in the population.
%    0 <= M <= L.
%
%    Input, integer L, the number of balls to select from.
%    0 <= L.
%
%    Output, real CDF, the value of the CDF.
%
  c1_log = binomial_coef_log ( l - m, n );
  c2_log = binomial_coef_log ( l, n );

  pdf = exp ( c1_log - c2_log );
  cdf = pdf;

  for x2 = 0 : x - 1

    pdf = pdf * ( m - x2 ) * ( n - x2 ) ...
      / ( ( x2 + 1 ) * ( l - m - n + x2 + 1 ) );

    cdf = cdf + pdf;

  end

  return
end


text:
function res = imProjectedArea(img, shifts, varargin)
% Total projected area of a 3D region in a given direction.
%
%   AREA = imProjectedArea(IMG, SHIFT)
%   IMG is a 3D binary image, SHIFT is a 1-by-3 row vector indicating the
%   shift between two voxels to test.
%
%   Example
%     % generate binary image of a 3D ellipsoid
%     elli = [50.12 50.23 50.34 50 35 20 30 40 50];
%     img = discreteEllipsoid(1:100, 1:100, 1:100, elli);
%     % compute projected area in main directions
%     D1 = imProjectedArea(img, [1 0 0]);
%     D2 = imProjectedArea(img, [0 1 0]);
%     D3 = imProjectedArea(img, [0 0 1]);
%     % compute projected area in a less common direction
%     D4 = imProjectedArea(img, [2 1 1]);
%
%   See also
%     imSurface, imProjectedDiameter

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2015-05-27,    using Matlab 8.4.0.150421 (R2014b)
% Copyright 2015 INRAE - Cepia Software Platform.

dim = size(img);

dx = shifts(1);
dy = shifts(2);
dz = shifts(3);

dl = hypot(hypot(dx, dy), dz);
% vol = dx * dy * dz;

% iterate over pixels in image to count number of transitions
% count = 0;
% for z = 3:dim(3)-2
%     for y = 3:dim(1)-2
%         for x = 3:dim(2)-2
%             if img(y, x, z) ~= img(y+dy, x+dx, z+dz)
%                 count = count + 1;
%             end
%         end
%     end
% end

ix = 3:dim(2)-2;
iy = 3:dim(1)-2;
iz = 3:dim(3)-2;
count = sum(sum(sum( img(iy, ix, iz) ~= img(iy+dy, ix+dx, iz+dz) )));

% number of connected components
count = count / 2;

% normalize with line density
res = count * 1 / dl;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [J, grad] = costFunction(theta, X, y)
%COSTFUNCTION Compute cost and gradient for logistic regression
%   J = COSTFUNCTION(theta, X, y) computes the cost of using theta as the
%   parameter for logistic regression and the gradient of the cost
%   w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Note: grad should have the same dimensions as theta
%
z=X*theta;
h=sigmoid(z);
%logisf=-y*log(h)-(1-y)*log(1-h);
%logisf=log(h)*(-y)-log(1-h)*(1-y);
logisf=(-y)'*log(h)-(1-y)'*log(1-h);
J=(1/m)*sum(logisf);		% sum will only work for a column, will not sum over a row

grad=1/m*((X'*h-X'*y)');

% =============================================================

end

author: yhyap

repo: machine-learning-coursera


text:
function [hyp, flogtheta, i,simy,simse2] = trainGPoe(hyp, inf, mean, cov, lik, input, target, lag, simf, Nsamples)
% Function for the optimisation (training) of GP model hyperparameters
%
%% Syntax
% [hyp, flogtheta, i] = trainGPoe(hyp, inf, mean, cov, lik, input, target,
%                                 simf, lag, Nsamples);
%
%% Description
% Function for the optimisation (training) of GP model hyperparameters
% based on the training data via maximum marginal likelihood.
% Uses routines gp and minimize.
% Based on the work of C.E.Rasmussen.
%
% Input:
% * hyp      ... the structure of initial hyperparameters
% * inf      ... the function specifying the inference method
% * cov      ... the prior covariance function (see below)
% * mean     ... the prior mean function
% * lik      ... the likelihood function
% * input    ... the input part of the training data,  NxD matrix
% * target   ... the output part of the training data (ie. target), Nx1 vector
% * simf     ... the function handle of simulation function (e.g. @simulGPmc)
% * lag      ... the order of the model (number of used lagged outputs)
% * Nsamples ... the number of samples for MCMC simulation (optional)
% * minf     ... the function handle of the minimization method to be used
%                (optional, default=@minimize_and_return). The minimization
%				  method should return the following values:
%				  [X, fX, i,fout] = minimize(X, f, length, P1, P2, P3, ... ),
%				  where fout is the additional argument (along fitness and gradient)
%				  of the evaulation function f
%
% Output:
% * hyp       ... optimized hyperparameters
% * flogtheta ... the minus log likelihood for the different runs (init. to 0)
% * i         ... the number of iterations needed for the last optimization
%
% Examples:
% demo_example_gp_training.m
%
% See Also:
% gp, minimize, covFunctions, trainlgmp
%

if(nargin < 11)
  minf = @minimize;
end
if (nargin < 10)
  Nsamples = 100; % default value of Nsamples for MC simulation
end
if (nargin < 9)
  simf = @simulGPnaive; %choose the default simulation routine
end

if (nargin < 8)
  warning(['The number of dynamical system model order was not provided. ' ...
  'Guessing from the input data...']);
  lag=size(input,2)/2;
  assert(mod(lag,1)==0);
  fprintf('Guess: lag=%d\n',lag);
end

if (nargin < 7)
  error('Too few parameters are given.');
end

MIN_DIFF = 0.002;
flogtheta = [0,-1/eps];

% we reconstruct the input signal from the regressor matrix as
%  the two following commented lines show:
%%N = length(input,1) +    lag-1
%%[u(1) ... u(N-lag)]';	[u(N-lag+1) ...  u(N)]

u=[input(:,lag+1)      ;	input(end,lag+2:end);NaN];

simulated_y=[input(1:lag,1);target];

while (abs(flogtheta(end) - flogtheta(end-1))>MIN_DIFF)
	disp(' ');
	disp(strcat(['delta flogtheta: ', num2str(abs(flogtheta(end) - flogtheta(end-1)))]));
	disp(' ')
	[hyp, flogthetatmp,i] = feval(minf, hyp, @simLL, -100,...
				inf, mean, cov, lik, u, simulated_y, target, lag, simf);

	%in the last iteration of minimisation algorithm the simulation is evaluated again for obtaining
	% the simulated response 'simulated_y'
	disp('saving simulated response.')
    [~,~,simulated_y] = simLL(hyp,inf, mean, cov, lik, u, simulated_y, target, lag, simf);

    if isempty(flogthetatmp) % no improvement: at minimum
        disp('oops')
        break
    end
	flogtheta = [flogtheta flogthetatmp(end)];
end

function [nlZ,dnlZ,simulated_y]=simLL(hyp, inf, mean, cov, lik, u, simulated_y, target,lag,simf)
	% simLL simulates the response and forms the regressor vectors from the response. Despite the response signal in regressor vectors is the simulated one, the regressands are always the original measurements.
	%
	% log-likelihood is calculated from the gp model with input data from simulation (in regressor matrix) and measurements (as regressands).
	%
	% simulated_y is a NxS or 2NxS matrix where N is the signal length and S is the number of past iterations. The values from N+1 to 2N represent the predicted variances.
	%
	N=size(target,1)+lag;
	if size(simulated_y,1)==N
		simulated_y=[simulated_y;simulated_y*NaN];
	end

	simulated_input0 = construct(lag, u, simulated_y(1:N,end));

	[y_mu, se2] = feval(simf,hyp, inf, mean, cov, lik, simulated_input0, ...
	target, simulated_input0, lag);

	%y_mu_full combines measured initial values and simulated response:
	simulated_y(1:N,end+1)=[simulated_y(1:lag,end);y_mu(:)];
	simulated_y(N+1:2*N,end)=[zeros(lag,1);se2(:)];

	simulated_input1 = construct(lag, u, simulated_y(1:N,end));
	[nlZ,dnlZ]=gp(hyp,inf,mean,cov,lik,simulated_input1,target);
end

simy=simulated_y(1:end/2,:);
simse2=simulated_y(end/2:end,:);
end

author: Dynamic-Systems-and-GP

repo: GPdyn


text:
%% im2MRcart
% Below is a demonstration of the features of the |im2MRcart| function

%%
clear; close all; clc;

%% Syntax
% |[Xc,Yc,Zc]=im2MRcart(I,J,K,v,OR,r,c);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function fileList = image_search(dirname, search_filter_handle)
% IMAGE_SEARCH Searches the specified directory for images satisfying the
% given filter.
%
% This function reads the metadata from each image in the specified
% directory and calls the provided function handle on each one.  If the
% (boolean) result of the function is true the file name and metadata are
% returned in the resulting list.  If it evaluates to false, the image is
% ignored.
%
% Inputs:
%       dirname:        The string name of the directory in which the
%                       images live.  All images will be opened, checked,
%                       and closed.
%
%       search_filter_handle:  The search filter to be applied to each
%                       image found.  This is a function pointer (handle)
%                       to a function that takes on SICD structure and
%                       returns a boolean value.  True means "keep the
%                       image", false means "ignore the image".
%
% Return:
%       IMAGE_SEARCH returns an array, each element being a structure
%       with elements two elements:
%               filename: The full path to the image
%               meta:     The SICD metadata structure for the image
%
% Written by: Tom Krauss, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

fullFileNameList = {};
metadataList     = {};
wb_hand=waitbar(0,'Enumerating files to search');
allFilesToReview = recurdir(dirname);
for i = 1:length(allFilesToReview)
    waitbar((i-1)/length(allFilesToReview),wb_hand,'Searching for recognized formats');
    try
        % disp(allFilesToReview{i}); % For debugging
        filename = fullfile(dirname,allFilesToReview{i});
        if ~isempty(guess_ph_format(filename))
            newreader = open_ph_reader(filename);
            % Maybe should convert to SICD metadata here (meta2sicd_cphdx)
            % for consistent filtering?  For now, leave as CPHD.
        else
            newreader = open_reader(filename);
        end
        if iscell(newreader)
            % Multiple images in file.  Assume each is an independant image and
            % tack on its metadata
            for iii=1:length(newreader)
                metadataList{end+1}     = newreader{iii}.get_meta();
                if isfield(metadataList{end},'native')
                    metadataList{end}=rmfield(metadataList{end},'native'); % Helps clear memory
                end
                fullFileNameList{end+1} = allFilesToReview{i};
                newreader{iii}.close();
            end
        else
            metadataList{end+1}     = newreader.get_meta();
            if isfield(metadataList{end},'native')
                metadataList{end}=rmfield(metadataList{end},'native'); % Helps clear memory
            end
            fullFileNameList{end+1} = allFilesToReview{i};
            newreader.close();
        end
    catch % Nothing to do.  Not a recognized format.  Just go to next file
    end
    clear newreader; % Helps clear memory
end

% Build the output array...
% Check each image to see if it passes the user-supplied filter.  If it
% does we'll append the file name to the end of the returned list,
% otherwise we'll just ignore/drop it.
fileList = [];
for i = 1:length(metadataList)
    filename = fullFileNameList{i};
    sicdmeta = metadataList{i};
    try
        if search_filter_handle(sicdmeta)
            result.filename = filename;
            result.meta     = sicdmeta;
            fileList = [fileList result];
        end
    end
    waitbar(1,wb_hand,'Filtering results');
end
close(wb_hand);

end

%RECURDIR List all files found recursively under given directory
function [ filelist ] = recurdir( dir2search )

if(nargin<1||isempty(dir2search))
    d = dir;
    dir2search=[];
else
    d = dir(dir2search);
end

filelist={};
for j=1:length(d)
    if ~d(j).isdir
        filelist{end+1,1}=d(j).name; % Record this file
    elseif ~strcmp(d(j).name(1),{'.','..','.svn'})
        sublist=recurdir(fullfile(dir2search, d(j).name)); % Find all of the files under this subdirectory
        sublist=cellfun(@(x) [d(j).name filesep x],sublist,'UniformOutput',false); % Append current directory
        filelist=cat(1,filelist,sublist); % Append to all previous files
    end
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function spm_voice_P300
% Illustrate voice recognition with lexical priors
% FORMAT spm_voice_P300
%
% loads the global variable VOX.mat
%
% VOX.LEX(w,k) -  structure array for k variants of w words
% VOX.PRO(p)   -  structure array for p aspects of prosody
% VOX.WHO(w)   -  structure array for w aspects of idenity
%
% This routine demonstrates the basic functionality of voice recognition or
% active listening with a special focus on segmentation and the simulated
% neurophysiological correlates of belief updating. It starts by
% demonstrating segmentation; either in response to some spoken sentences
% (read from prompts in the script or by loading exemplar sentences). It
% then moves on to demonstrating the effect of changing the precision of
% prior beliefs about lexical content and how this is expressed in terms of
% simulated belief updating via the minimisation of variational free
% energy.
%
% This routine assumes the necessary files are located in a particular
% (Sound files) directory; that can be specified by editing the script
% below.
%__________________________________________________________________________
% Copyright (C) 2019 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_voice_P300.m 7750 2019-12-05 17:54:29Z spm $

%% demo mode loads sentence (.mat) files
%==========================================================================
DEMO = 1;

% get lexical and prosody arrays in sound file directory
%--------------------------------------------------------------------------
clear global VOX
global VOX

NAME  = which('VOX.mat');
PATH  = fileparts(NAME);
load(NAME)

VOX.analysis = 0;
VOX.graphics = 0;
VOX.interval = 0;
VOX.mute     = 0;
VOX.onsets   = 0;
VOX.depth    = 0;

%% set up priors a succession of 'triangle' or 'square'
%==========================================================================
nw    = numel(VOX.LEX);                                   % number of words
[i,P] = spm_voice_i(repmat({{'triangle','square'}},1,6)); % prior words

% record corresponding sequence and save - or load a preprepared sentence
%--------------------------------------------------------------------------
if DEMO

    load(fullfile(PATH,'triangle_square.mat'))

else

    % SAY: {'triangle','square','triangle','square','triangle','square'};
    %----------------------------------------------------------------------
    spm_voice_read
    Y  = getaudiodata(VOX.audio);
    save(fullfile(PATH,'triangle_square.mat'),'Y');
    [FS,read] = spm_voice_FS(wfile)

end

% illustrate candidate intervals (and F0) for the first word
%--------------------------------------------------------------------------
VOX.FS       = spm_voice_FS;
VOX.IT       = 1;
VOX.onsets   = 1;
VOX.interval = 1;
spm_voice_get_word(Y);
VOX.interval = 0;
spm_voice_fundamental(Y,VOX.FS);
VOX.onsets   = 0;

% segment without priors
%--------------------------------------------------------------------------
VOX.depth    = 0;
spm_figure('GetWin','Segmentation: no priors');  clf;
SEG0 = spm_voice_read(Y);
EEG0 = spm_voice_segmentation(Y,SEG0);

% segment with priors
%--------------------------------------------------------------------------
VOX.depth    = 1;
spm_figure('GetWin','Segmentation: with priors'); clf;
SEG1 = spm_voice_read(Y,P);
EEG1 = spm_voice_segmentation(Y,SEG1);
VOX.depth    = 0;

% illustrate the role of prior precision
%==========================================================================

%% lexical priors; i.e., plausible sequence of words
%--------------------------------------------------------------------------
str{1} = {'is'};
str{2} = {'there'};
str{3} = {'a'};
str{4} = {'triangle','square'};
str{5} = {'below','above','there'};
[i,P]  = spm_voice_i(str);                   % get priors
P      = spm_softmax(log(P));                % and double their precision

% record corresponding sequence and save - or load a preprepared sentence
%--------------------------------------------------------------------------
if DEMO

    load(fullfile(PATH,'sentence.mat'))

else

    % SAY: {'is','there','a','square','above'};
    %----------------------------------------------------------------------
    spm_voice_read
    Y  = getaudiodata(VOX.audio);
    save(fullfile(PATH,'sentence.mat'),'Y');

end

%% change VOX.noise to simulate speech in noise (e.g., VOX.noise = 2)
%--------------------------------------------------------------------------
VOX.noise = 1;                                  % noise level (default: 1)
s         = 3;                                  % index of surprising word

%% segment with priors
%--------------------------------------------------------------------------
spm_figure('GetWin','Segmentation (P300): with priors'); clf;
SEG1  = spm_voice_read(Y,P);
EEG1  = spm_voice_segmentation(Y,SEG1);

% segment without priors
%--------------------------------------------------------------------------
spm_figure('GetWin','Segmentation (P300): no priors'); clf;
SEG0  = spm_voice_read(Y);
EEG0  = spm_voice_segmentation(Y,SEG0);

% add evoked responses with priors, to highlight more exuberant ERPs
%--------------------------------------------------------------------------
subplot(4,1,4), hold on
plot(VOX.PST,EEG1,'-.');

% remove priors from a single word (s)
%--------------------------------------------------------------------------
spm_figure('GetWin','Segmentation: Bayesian surprise'); clf;
Q     = P; Q(:,s) = ones(nw,1)/nw;
SEG2  = spm_voice_read(Y,Q);
EEG2  = spm_voice_segmentation(Y,SEG2);

% illustrate surprise responses (e.g.,P300)
%==========================================================================

% simulations of P300: responses to a single word (s)
%--------------------------------------------------------------------------
spm_figure('GetWin','P300'); clf;
i   = spm_voice_i(SEG1(s).str);
E0  = EEG0(:,i);
E1  = EEG1(:,i);
E2  = EEG2(:,i);

% peristimulus time for plotting (250 ms before and after offset)
%--------------------------------------------------------------------------
t   = SEG2(s).IT/VOX.FS - 1/4;
x   = [t,t + 1/2,t + 1/2,t];
y   = [-1,-1,1,1];
c   = spm_softmax(rand(3,1))';

% responses with and without priors - all words
%--------------------------------------------------------------------------
subplot(3,2,1), plot(VOX.PST,E0,'-.',VOX.PST,E1)
xlabel('time (sec)'), ylabel('amplitude'), title('Waveforms','FontSize',16)
axis square, spm_axis tight, set(gca,'YLim',[-1 1]), box off

% response differentials - all words
%--------------------------------------------------------------------------
subplot(3,2,3), plot(VOX.PST,E0 - E1)
xlabel('time (sec)'), ylabel('amplitude'), title('Differences','FontSize',16)
axis square, spm_axis tight, set(gca,'YLim',[-1 1]/3), box off

% responses with and without priors - one word
%--------------------------------------------------------------------------
subplot(3,2,2), plot(VOX.PST,E2,'-.',VOX.PST,E1)
xlabel('time (sec)'), ylabel('amplitude'), title('Waveforms','FontSize',16)
axis square, spm_axis tight, set(gca,'YLim',[-1 1]), box off

% show peristimulus time
%--------------------------------------------------------------------------
hold on, h = fill(x,y,c); hold off
set(h,'Facealpha',1/8,'EdgeAlpha',1/8);

% response differentials - one word
%--------------------------------------------------------------------------
subplot(3,2,4), plot(VOX.PST,E2 - E1)
xlabel('time (sec)'), ylabel('amplitude'), title('Differences','FontSize',16)
axis square, spm_axis tight, set(gca,'YLim',[-1 1]/3), box off

% show peristimulus time
%--------------------------------------------------------------------------
hold on, h = fill(x,y/3,c); hold off
set(h,'Facealpha',1/8,'EdgeAlpha',1/8);

%% illustrate the relationship between belief updating and RMS responses
%==========================================================================
for i = 1:numel(SEG2)
    q     = spm_softmax(SEG2(i).L{1});
    p     = SEG2(i).p;
    KL(i) = q'*(log(q + exp(-8)) - log(p + exp(-8)));
end

% show RMS responses as a function of time
%--------------------------------------------------------------------------
spm_figure('GetWin','Segmentation: Bayesian surprise');
RMS  = var(EEG2,[],2);
RMS  = 4*RMS/max(RMS);
subplot(4,2,3), hold off, plot(VOX.PST,RMS), hold on
xlabel('time (sec)'), ylabel('Power/KL (nats)')
title('Evoked power','FontSize',16)
spm_axis tight

% and overlay associated belief updating in terms of KL divergence
%--------------------------------------------------------------------------
for i = 1:numel(SEG2)
    t = SEG2(i).IT/VOX.FS + 1/8;
    plot([t,t],[0 KL(i)],'r','LineWidth',4)
end

i     = find(diff(RMS(1:end - 1)) > 0 & diff(RMS(2:end)) < 0);
[R,i] = sort(RMS(i),'descend');
R     = R(1:numel(KL));
KL    = sort(KL,'descend');
B     = [ones(numel(KL),1) KL(:)]\R;

subplot(4,2,4), plot(KL,R,'.r','MarkerSize',16), hold on
kl    = -1:4;
plot(kl,B(1) + kl*B(2),':r','MarkerSize',16)
xlabel('Bayesian suprise (KL - nats)'), ylabel('Evoked power')
title('Belief updating','FontSize',16)
spm_axis square

author: spm


text:
function [A,B,D,E,I]=mdds(X)
% gets A, B, D, E & I arrays for magnet driver
N=2;U=3;M=1;Y=1;
A1=zeros(U+N);B2=zeros(U+N,N+M);
%
% X = [R1 R2 R3 C1 C2]
%
R1=X(1);R2=X(2);R3=X(3);C1=X(4);C2=X(5);
%
% V1 V2 V3 iC1 iC2  column labels for A1
% 1  2  3  4	5
%
A1(1,1)=-1/R1;A1(1,4)=-1;A1(1,5)=-1; % From 1st equation
A1(2,1)=1/R2;A1(2,2)=-1/R2;A1(2,4)=-1; % 2nd equation
A1(3,3)=1;A1(3,4)=-R3;A1(3,5)=-R3; % V3 - iC1*R3 - iC2*R3 = 0
A1(4,2)=1;A1(4,3)=-1; % V2 - V3 = E1
A1(5,1)=1;A1(5,3)=-1; % V1 - V3 = E2
%
E1=1;E2=1;Ein=1;
% E1  E2   Ein  Column labels for B2
%  1   2    3
B2(4,1)=E1;
B2(5,2)=E2;
B2(1,3)=-Ein/R1;
%
P=diag([C1 C2]);
%
% Template statements (same for every circuit):
V=A1\B2;H=V(U+1:U+N,1:N+M);I=eye(N);
AB=P\H;A=AB(1:N,1:N);B=AB(1:N,N+1:N+M);
D=V(Y:Y,1:N);E=V(Y:Y,N+1:N+M);

author: Sable


text:
%% Pendulum example
%  This example applies energy-preserving piDMD to reconstruct the trajectory of
%  a double pendulum from noisy measurements.
addpath('../src')
rng(1); % Set random seed

% Define parameters of problem
l1=1; l2=1.5; % Lengths of rods
m1=1 ; m2=1.5; g=9.81; % Masses and gravity
params = [l1, l2, m1, m2, g]; % Concatenate parameters

% Construct linearised energy inner product
W(1,1) = (m1/2+m2/2)*g*l1;
W(2,2) = m2/2*g*l2;
W(3,3) = (m1/2+m2/2)*l1^2;
W(4,3) = m2/2*l1*l2;
W(3,4) = W(4,3);
W(4,4) = m2/2*l2^2;
C = chol(W); % Calculate inner product

% Set number of samples and span
tend = 30; nt = 1000;
tspan= linspace(0,tend,nt);

% Set initial conditions
theta1= 0.4; theta1_prime=0;
theta2= 0.7; theta2_prime=0;
y0=[theta1 theta1_prime theta2 theta2_prime];

% Solve ODE
[t,y]=ode45(@(t,y) pendulum(t,y,params), tspan,y0);

% Extract data
th1 = y(:,1); th2 = y(:,3); th1dt = y(:,2); th2dt = y(:,4);
x = [th1'; th2'; th1dt'; th2dt'];
xn = x + 1e-1*std(x,[],2).*randn(size(x)); % Add noise
data = C*xn; % Rescale measurements into energy norm
nTrain = nt-1;
X = data(:,1:nTrain); Y = data(:,2:nTrain+1);

% Train models
[piA,piVals] = piDMD(X,Y,'orthogonal');
[exA,exVals] = piDMD(X,Y,'exact');

% Perform reconstructions
piRec = zeros(4,nt); piRec(:,1) = data(:,1);
exRec = zeros(4,nt); exRec(:,1) = data(:,1);
for j = 2:nt
piRec(:,j) = piA(piRec(:,j-1));
exRec(:,j) = exA(exRec(:,j-1));
end

% Rescale reconstructions back into physical norm
fpiRec = C\piRec; fexRec = C\exRec;

%% Plot results
figure(1); LW = 'LineWidth'; IN = 'Interpreter'; LT = 'Latex'; FS = 'FontSize';
subplot(3,1,1)
c1 = .8*[1 1 1]; c2 = .8*[1 1 1];
plot(tspan,xn(1,:),LW,2,'Color',c1)
hold on
plot(tspan,xn(2,:),LW,2,'Color',c2)
hold off; trajPlot('measurements'); xticklabels([])
subplot(3,1,2)
plot(tspan,x(1,:),LW,3,'Color', c1)
hold on
plot(tspan,fpiRec(1,:),'b--',LW,2)
plot(tspan,fexRec(1,:),'r--',LW,2)
ylabel('$\theta_1$',IN,LT)
hold off; trajPlot('$\theta_1$'); xticklabels([])
subplot(3,1,3)
l1=plot(t,x(2,:),LW,3,'Color', c2);
hold on
l2=plot(t,fpiRec(2,:),'b--',LW,2);
l3=plot(t,fexRec(2,:),'r--',LW,2);
hold off; xlabel('time',FS,20,IN,LT)
trajPlot('$\theta_2$')
legend([l1,l2,l3],{'truth','piDMD','exact DMD'},IN,LT)
function yp = pendulum(~, y, params)

l1=params(1);  l2=params(2);
m1=params(3);  m2=params(4);
g=params(5);

a = (m1+m2)*l1 ;
b = m2*l2*cos(y(1)-y(3)) ;
c = m2*l1*cos(y(1)-y(3)) ;
d = m2*l2 ;
e = -m2*l2*y(4)* y(4)*sin(y(1)-y(3))-g*(m1+m2)*sin(y(1)) ;
f = m2*l1*y(2)*y(2)*sin(y(1)-y(3))-m2*g*sin(y(3)) ;
yp=zeros(4,1);
yp(1) = y(2);
yp(3)= y(4) ;
yp(2)= (e*d-b*f)/(a*d-c*b) ;
yp(4)= (a*f-c*e)/(a*d-c*b) ;
end

function trajPlot(j) % Nice plot of trajectories
yticks([-pi/4,0,pi/4]); yticklabels([{'$-\pi/4$'},{'0'},{'$\pi/4$'}])
set(gca,'TickLabelInterpreter','Latex','FontSize',20);grid on
ylim([-1,1])
ylabel(j,'Interpreter','latex','FontSize',20)
end

author: baddoo

repo: piDMD

save_path: github-repos/MATLAB/baddoo-piDMD