text:
function getPart()
partNum=5;
partList={'bird','cat','cow','dog','horse','sheep'};
theConf=configurations;

for i=1:numel(partList)
    Name_batch=partList{i};
    objset=readAnnotation(Name_batch,theConf);

    tarSize=[224,224];
    getIMDBMask(Name_batch,objset,partNum,tarSize,theConf);
    tarSize=[227,227];
    getIMDBMask(Name_batch,objset,partNum,tarSize,theConf);
end
end

function getIMDBMask(Name_batch,objset,partNum,tarSize,theConf)
objNum=numel(objset);
mask=zeros(tarSize(1),tarSize(2),objNum,'uint8');
partPos=cell(objNum,1);
for i=1:objNum
    tmp=getPart_cate(objset(i),theConf,tarSize);
    mask(:,:,i)=uint8(tmp);
    pos=ones(2,partNum).*(-1);
    for part=1:partNum
        [hL,wL]=find(tmp==part);
        if(isempty(hL))
            continue;
        end
        pos(1,part)=mean(hL);
        pos(2,part)=mean(wL);
    end
    partPos{i}=pos;
    fprintf('%s   scale %d   %d\n',Name_batch,tarSize(1),i);
end
filename=sprintf('../results/imdb/%d/%s/partMask.mat',tarSize(1),Name_batch);
save(filename,'mask','partPos');
end

function mask=getPart_cate(obj,theConf,tarSize)
Batch_Name = obj.name;
switch(Batch_Name)
    case 'bird'
        partMerge={[1,2,3,4],[5,6,7,8],[9,10,11,12],13};
    case 'cat'
        partMerge={[1,2,3,4,5,6],[7,8],[9,10,11,12],[13,14,15,16],17};
    case 'cow'
        partMerge={[1,2,3,4,5,6,7,8],[9,10],[11,12,13,14],[15,16,17,18]};
    case 'dog'
        partMerge={[1,2,3,4,5,6,20],[7,8],[9,10,11,12],[13,14,15,16],17};
    case 'horse'
        partMerge={[1,2,3,4,5,6],[9,10],[11,12,13,14,30,31],[15,16,17,18,32,33]};
    case 'sheep'
        partMerge={[1,2,3,4,5,6,7,8],[9,10],[11,12,13,14],[15,16,17,18]};
    otherwise
        error('Errors here.');
end
[mask,bndbox]=getMask(obj,theConf,partMerge);
mask=mask(bndbox.ymin:bndbox.ymax,bndbox.xmin:bndbox.xmax);
yL=round(linspace(1,size(mask,1),tarSize(1)));
xL=round(linspace(1,size(mask,2),tarSize(2)));
mask=mask(yL,xL);
end

function [mask,bndbox]=getMask(obj,theConf,partMerge)
for i=1:numel(partMerge)
    pL=partMerge{i};
    for j=1:numel(pL)
        [tmp,bndbox]=getPart_one(obj,pL(j),theConf);
        if((i==1)&&(j==1))
            mask=zeros(size(tmp));
        end
        mask(tmp>0)=i;
    end
end
end

function [part_mask,bndbox]=getPart_one(obj, partID, theConf)

images_path = [theConf.data.imgdir, obj.folder];
anno_file = [theConf.data.annodir,'/%s.mat'];
imname = obj.filename;
batch_name = obj.name;
bndbox = obj.bndbox;
bndbox.ymax=str2double(bndbox.ymax);
bndbox.ymin=str2double(bndbox.ymin);
bndbox.xmax=str2double(bndbox.xmax);
bndbox.xmin=str2double(bndbox.xmin);
batch2classID=containers.Map('KeyType','char','ValueType','int32');
batch2classID('bird')=3;
batch2classID('cat')=8;
batch2classID('cow')=10;
batch2classID('dog')=12;
batch2classID('horse')=13;
batch2classID('sheep')=17;
classID=batch2classID(batch_name);

pimap = part2ind();
img = imread([images_path, '/', imname]);
% load annotation -- anno
load(sprintf(anno_file, imname(1:end-4)));

[cls_mask, inst_mask, part_mask] = mat2map(anno, img, pimap);

tmp=(cls_mask==classID).*(part_mask==partID);
tmp([1:bndbox.ymin,bndbox.ymax:size(tmp,1)],:)=0;
tmp(:,[1:bndbox.xmin,bndbox.xmax:size(tmp,2)])=0;
inst_mask=inst_mask.*uint8(tmp);
part_mask=part_mask.*uint8(tmp);

for h=1:size(part_mask,1)
    for w=1:size(part_mask,2)
        if cls_mask(h,w)~=classID
            part_mask(h,w)=0;
            inst_mask(h,w)=0;
        elseif ((part_mask(h,w)~=partID) || (h>bndbox.ymax) || (h<bndbox.ymin) || (w<bndbox.xmin) || (w>bndbox.xmax))
            % skip if no part annotation or pixel outside of bndbox
            part_mask(h,w)=0;
            inst_mask(h,w)=0;
        end
    end
end
end

function [cls_mask, inst_mask, part_mask] = mat2map(anno, img, pimap)
% Read the annotation and present it in terms of 3 segmentation mask maps (
% i.e., the class maks, instance maks and part mask). pimap defines a
% mapping between part name and index (See part2ind.m).

cls_mask = zeros(size(img,1), size(img,2), 'uint8');
inst_mask = zeros(size(img,1), size(img,2), 'uint8');
part_mask = zeros(size(img,1), size(img,2), 'uint8');

for oo = 1:numel(anno.objects)
    % The objects are ordered such that later objects can occlude previous
    % objects
    obj = anno.objects(oo);
    class_ind = obj.class_ind;
    silh = obj.mask;            % the silhouette mask of the object
    assert(size(silh,1) == size(img,1) && size(silh,2) == size(img,2));
    inst_mask(silh) = oo;
    cls_mask(silh) = class_ind;

    for pp = 1:numel(obj.parts)
        part_name = obj.parts(pp).part_name;
        assert(isKey(pimap{class_ind}, part_name));     % must define part index for every part
        assert(all(silh(obj.parts(pp).mask)));          % all part region is a subregion of the object
        pid = pimap{class_ind}(part_name);
        part_mask(obj.parts(pp).mask) = pid;
    end
end
end

function pimap = part2ind()
% Define the part index of each objects.
% One can merge different parts by using the same index for the
% parts that are desired to be merged.
% For example, one can merge
% the left lower leg (llleg) and the left upper leg (luleg) of person by setting:
% pimap{15}('llleg')      = 19;               % left lower leg
% pimap{15}('luleg')      = 19;               % left upper leg

pimap = cell(20, 1);
% Will define part index map for the 20 PASCAL VOC object classes in ascending
% alphabetical order (the standard PASCAL VOC order).
for ii = 1:20
    pimap{ii} = containers.Map('KeyType','char','ValueType','int32');
end

% [aeroplane]
pimap{1}('body')        = 1;
pimap{1}('stern')       = 2;
pimap{1}('lwing')       = 3;                % left wing
pimap{1}('rwing')       = 4;                % right wing
pimap{1}('tail')        = 5;
for ii = 1:10
    pimap{1}(sprintf('engine_%d', ii)) = 10+ii; % multiple engines
end
for ii = 1:10
    pimap{1}(sprintf('wheel_%d', ii)) = 20+ii;  % multiple wheels
end

% [bicycle]
pimap{2}('fwheel')      = 1;                % front wheel
pimap{2}('bwheel')      = 2;                % back wheel
pimap{2}('saddle')      = 3;
pimap{2}('handlebar')   = 4;                % handle bar
pimap{2}('chainwheel')  = 5;                % chain wheel
for ii = 1:10
    pimap{2}(sprintf('headlight_%d', ii)) = 10 + ii;
end

% [bird]
pimap{3}('head')        = 1;
pimap{3}('leye')        = 2;                % left eye
pimap{3}('reye')        = 3;                % right eye
pimap{3}('beak')        = 4;
pimap{3}('torso')       = 5;
pimap{3}('neck')        = 6;
pimap{3}('lwing')       = 7;                % left wing
pimap{3}('rwing')       = 8;                % right wing
pimap{3}('lleg')        = 9;                % left leg
pimap{3}('lfoot')       = 10;               % left foot
pimap{3}('rleg')        = 11;               % right leg
pimap{3}('rfoot')       = 12;               % right foot
pimap{3}('tail')        = 13;

% [boat]
% only has silhouette mask

% [bottle]
pimap{5}('cap')         = 1;
pimap{5}('body')        = 2;

% [bus]
pimap{6}('frontside')   = 1;
pimap{6}('leftside')    = 2;
pimap{6}('rightside')   = 3;
pimap{6}('backside')    = 4;
pimap{6}('roofside')    = 5;
pimap{6}('leftmirror')  = 6;
pimap{6}('rightmirror') = 7;
pimap{6}('fliplate')    = 8;                % front license plate
pimap{6}('bliplate')    = 9;                % back license plate
for ii = 1:10
    pimap{6}(sprintf('door_%d',ii)) = 10 + ii;
end
for ii = 1:10
    pimap{6}(sprintf('wheel_%d',ii)) = 20 + ii;
end
for ii = 1:10
    pimap{6}(sprintf('headlight_%d',ii)) = 30 + ii;
end
for ii = 1:20
    pimap{6}(sprintf('window_%d',ii)) = 40 + ii;
end

% [car]
keySet = keys(pimap{6});
valueSet = values(pimap{6});
pimap{7} = containers.Map(keySet, valueSet);         % car has the same set of parts with bus

% [cat]
pimap{8}('head')        = 1;
pimap{8}('leye')        = 2;                % left eye
pimap{8}('reye')        = 3;                % right eye
pimap{8}('lear')        = 4;                % left ear
pimap{8}('rear')        = 5;                % right ear
pimap{8}('nose')        = 6;
pimap{8}('torso')       = 7;
pimap{8}('neck')        = 8;
pimap{8}('lfleg')       = 9;                % left front leg
pimap{8}('lfpa')        = 10;               % left front paw
pimap{8}('rfleg')       = 11;               % right front leg
pimap{8}('rfpa')        = 12;               % right front paw
pimap{8}('lbleg')       = 13;               % left back leg
pimap{8}('lbpa')        = 14;               % left back paw
pimap{8}('rbleg')       = 15;               % right back leg
pimap{8}('rbpa')        = 16;               % right back paw
pimap{8}('tail')        = 17;

% [chair]
% only has sihouette mask

% [cow]
pimap{10}('head')       = 1;
pimap{10}('leye')       = 2;                % left eye
pimap{10}('reye')       = 3;                % right eye
pimap{10}('lear')       = 4;                % left ear
pimap{10}('rear')       = 5;                % right ear
pimap{10}('muzzle')     = 6;
pimap{10}('lhorn')      = 7;                % left horn
pimap{10}('rhorn')      = 8;                % right horn
pimap{10}('torso')      = 9;
pimap{10}('neck')       = 10;
pimap{10}('lfuleg')     = 11;               % left front upper leg
pimap{10}('lflleg')     = 12;               % left front lower leg
pimap{10}('rfuleg')     = 13;               % right front upper leg
pimap{10}('rflleg')     = 14;               % right front lower leg
pimap{10}('lbuleg')     = 15;               % left back upper leg
pimap{10}('lblleg')     = 16;               % left back lower leg
pimap{10}('rbuleg')     = 17;               % right back upper leg
pimap{10}('rblleg')     = 18;               % right back lower leg
pimap{10}('tail')       = 19;

% [diningtable]
% only has silhouette mask

% [dog]
keySet = keys(pimap{8});
valueSet = values(pimap{8});
pimap{12} = containers.Map(keySet, valueSet);         	% dog has the same set of parts with cat,
                                            		% except for the additional
                                            		% muzzle
pimap{12}('muzzle')     = 20;

% [horse]
keySet = keys(pimap{10});
valueSet = values(pimap{10});
pimap{13} = containers.Map(keySet, valueSet);        	% horse has the same set of parts with cow,
                                                        % except it has hoof instead of horn
remove(pimap{13}, {'lhorn', 'rhorn'});
pimap{13}('lfho') = 30;
pimap{13}('rfho') = 31;
pimap{13}('lbho') = 32;
pimap{13}('rbho') = 33;

% [motorbike]
pimap{14}('fwheel')     = 1;
pimap{14}('bwheel')     = 2;
pimap{14}('handlebar')  = 3;
pimap{14}('saddle')     = 4;
for ii = 1:10
    pimap{14}(sprintf('headlight_%d', ii)) = 10 + ii;
end

% [person]
pimap{15}('head')       = 1;
pimap{15}('leye')       = 2;                    % left eye
pimap{15}('reye')       = 3;                    % right eye
pimap{15}('lear')       = 4;                    % left ear
pimap{15}('rear')       = 5;                    % right ear
pimap{15}('lebrow')     = 6;                    % left eyebrow
pimap{15}('rebrow')     = 7;                    % right eyebrow
pimap{15}('nose')       = 8;
pimap{15}('mouth')      = 9;
pimap{15}('hair')       = 10;

pimap{15}('torso')      = 11;
pimap{15}('neck')       = 12;
pimap{15}('llarm')      = 13;                   % left lower arm
pimap{15}('luarm')      = 14;                   % left upper arm
pimap{15}('lhand')      = 15;                   % left hand
pimap{15}('rlarm')      = 16;                   % right lower arm
pimap{15}('ruarm')      = 17;                   % right upper arm
pimap{15}('rhand')      = 18;                   % right hand

pimap{15}('llleg')      = 19;               	% left lower leg
pimap{15}('luleg')      = 20;               	% left upper leg
pimap{15}('lfoot')      = 21;               	% left foot
pimap{15}('rlleg')      = 22;               	% right lower leg
pimap{15}('ruleg')      = 23;               	% right upper leg
pimap{15}('rfoot')      = 24;               	% right foot

% [pottedplant]
pimap{16}('pot')        = 1;
pimap{16}('plant')      = 2;

% [sheep]
keySet = keys(pimap{10});
valueSet = values(pimap{10});
pimap{17} = containers.Map(keySet, valueSet);        % sheep has the same set of parts with cow

% [sofa]
% only has sihouette mask

% [train]
pimap{19}('head')       = 1;
pimap{19}('hfrontside') = 2;                	% head front side
pimap{19}('hleftside')  = 3;                	% head left side
pimap{19}('hrightside') = 4;                	% head right side
pimap{19}('hbackside')  = 5;                 	% head back side
pimap{19}('hroofside')  = 6;                	% head roof side

for ii = 1:10
    pimap{19}(sprintf('headlight_%d',ii)) = 10 + ii;
end

for ii = 1:10
    pimap{19}(sprintf('coach_%d',ii)) = 20 + ii;
end

for ii = 1:10
    pimap{19}(sprintf('cfrontside_%d', ii)) = 30 + ii;   % coach front side
end

for ii = 1:10
    pimap{19}(sprintf('cleftside_%d', ii)) = 40 + ii;   % coach left side
end

for ii = 1:10
    pimap{19}(sprintf('crightside_%d', ii)) = 50 + ii;  % coach right side
end

for ii = 1:10
    pimap{19}(sprintf('cbackside_%d', ii)) = 60 + ii;   % coach back side
end

for ii = 1:10
    pimap{19}(sprintf('croofside_%d', ii)) = 70 + ii;   % coach roof side
end

% [tvmonitor]
pimap{20}('screen')     = 1;
end

repo: interpretableCNN


text:
%................................................................

function outputDisplacementsReactionsPretty...
    (displacements,stiffness,GDof,prescribedDof,force)

% output of displacements and reactions in
% tabular form

% GDof: total number of degrees of freedom of
% the problem

disp('Displacements')
fprintf('node\t\tdisplacements\n')
% displacements
for jj=1:GDof
    fprintf('%2.0f:        %10.4e\n', jj, displacements(jj))
end

% reactions
F=stiffness*displacements;
reactions=F(prescribedDof)-force(prescribedDof);
disp('Reactions')
fprintf('node\t\treactions\n')
for jj=1:size(prescribedDof,1)
    fprintf('%2.0f:        %10.4e\n', prescribedDof(jj), reactions(jj))
end

repo: Aero-matlab


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = psabr_4_2(a, b, r, n, f, k, t,m, mu, nu, l, u)
% sabr prices using an admissible region kl <= k <= ku where sabr is used
% for 0 <= k <= kl    we use a put pricing function
%           f(x) = x^mu exp(a + bx + cx^2)
% for ku < k < +infty we use a call pricing function
%           f(x) = x^(-nu) exp(a + bx^(-1) + cx(-2)
% a big range (small kl and big ku) guarantee that the prices are in line
% with the observed calls and puts
%
% this is the preferred version of the density for computations
%

eps = 1e-004;                              % .1 bp

s = @(x) psabr(a, b, r, n, f, x, t);
index = find(s(k)>0,1,'first');

if (isempty(index) || index ==1)
   kl = l *f;
else
    kl = max(l * f,k(index)+(f-k(index))/m);% lower strike level
end
s1 = s(kl-eps);                          % for calc of derivatives
s2 = s(kl);                              % for calc of derivatives
s3 = s(kl+eps);                          % for calc of derivatives

V1 = log(s2);                            % log density

U2 = (s3-s1)/(2*eps);                    % derivative of density
V2 = U2/s2;                              % derivative of log density

U3 = (s3-2*s2+s1)/eps^2;                 % second derivative of density
V3 = U3/s2 - V2^2;                       % second derivative of log density

% fix mu and solve
%       V1 = mu log kl + a + b kl + c kl^2
%       V2 = mu / kl + b + 2c kl
%       V3 = - mu / kl^2 + 2c

%mu = 1.5;                               % controls the left tails
cl = .5*(V3+mu/kl^2);
bl = V2-mu/kl - 2*cl*kl;
al = V1 - mu * log(kl) - bl*kl - cl*kl^2;

% upper strike level (the bigger the strike level the better fit the call
% prices to Hagan formula!)
ku = u * f;
s = @(x) psabr(a, b, r, n, f, x, t);        % sabr price (1) for calls
s1 = s(ku-eps);                             % for calc of derivatives
s2 = s(ku);                                 % for calc of derivatives
s3 = s(ku+eps);                             % for calc of derivatives

V1 = log(s2);                               % log call price

U2 = (s3-s1)/(2*eps);                       % derivative of  density
V2 = U2/s2;                                 % derivative of log density

U3 = (s3-2*s2+s1)/eps^2;                    % second derivative of density
V3 = U3/s2 - (U2/s2)^2;                     % second derivative of density

% fix nu and solve
%       V1 = -nu log ku + a + b/ku + c/ku^2
%       V2 = -nu / ku - b / ku^2 - 2c/ku^3
%       V3 = nu / ku^2 + 2 b / ku^3 - 6 c / ku^4

cu = (-1.5*nu / ku + .5*V3 * ku - V2)*ku^3/5;
bu = -ku^2*(V2 + nu/ku +2*cu/ku^3);
au = V1 + nu * log(ku) - bu / ku - cu / ku^2;

    yl = real(k(k<kl).^mu .* exp(al + bl .* k(k<kl) + cl * k(k<kl).^2));
    ym = real(psabr(a,b,r,n,f,k((kl<=k)&(k<=ku)),t));
    yu = real(k(k>ku).^(-nu) .* exp(au+bu./k(k>ku)+cu./k(k>ku).^2));

y = [yl ym yu];                            % output
end

author: Sable

repo: mcbench-benchmarks


text:
hydro = struct();

hydro = readWAMIT(hydro,'rm3.out',[]);
hydro = radiationIRF(hydro,20,[],[],[],[]);
hydro = radiationIRFSS(hydro,[],[]);
hydro = excitationIRF(hydro,20,[],[],[],[]);
writeBEMIOH5(hydro)
plotBEMIO(hydro)

author: WEC-Sim

repo: WEC-Sim

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


text:
function S=triplyPeriodicMinimal(varargin)

% function S=triplyPeriodicMinimal(X,Y,Z,typeStr)
% ------------------------------------------------------------------------
% This function creates the image S which can be used to define triply
% periodic minimal surfaces. The input consists of a grid of coordinates
% (X,Y,Z) and the surface type (typeStr).
%
% 2009 Created
% 2016 Updated input handling
% 2018/12/13 Added comments at top of function
% ------------------------------------------------------------------------

%%
%Parse input
switch nargin
    case 2
        P=varargin{1};
        X=P(:,1); Y=P(:,2); Z=P(:,3); %Get coordinates
        typeStr=varargin{2}; %Get the surface type string
    case 4
        X=varargin{1};
        Y=varargin{2};
        Z=varargin{3};
        typeStr=varargin{4};
    otherwise
        error('False number of input arguments.');
end

%Evaluate metric on coordinates
switch typeStr
    case 'p' %Schwarz P
        S=(cos(X)+cos(Y)+cos(Z));
    case 'd' %Schwarz D
        S=(sin(X).*sin(Y).*sin(Z))...
            +(sin(X).*cos(Y).*cos(Z))...
            +(cos(X).*sin(Y).*cos(Z))...
            +(cos(X).*cos(Y).*sin(Z));
    case 'g' %Schoen Gyroid
        S=(sin(X).*cos(Y))+(sin(Y).*cos(Z))+(cos(X).*sin(Z));
    case 'n' %Neovius
        S=3*(cos(X)+ cos(Y)+ cos(Z))+ (4*cos(X).*cos(Y).*cos(Z));
    case 'w'
        S=2*(cos(X).*cos(Y)+cos(Z).*cos(X)+cos(Y).*cos(Z))-(cos(2*X)+cos(2*Y)+cos(2*Z));
    case 'pw'
        S=(4.*(cos(X).*cos(Y)+cos(Y).*cos(Z)...
            +cos(Z).*cos(X))-3.*cos(X).*cos(Y).*cos(Z))+2.4;
    otherwise
        error('unknown surface type requested')
end

%%
% _*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 spok_install
%SPOK_INSTALL compiles and installs the SPOK mexFunction
% Your current working directory must be the "spok" directory for this function
% to work.
%
% Example:
%   spok_install
%
% See also sparse, spok, spok_test

% Copyright 2008, Tim Davis, University of Florida

help spok_install
mex spok.c spok_mex.c
addpath (pwd)

fprintf ('Added the following directory to your path:\n') ;
disp (pwd) ;
fprintf ('Use pathtool and click "save" to save this for future sessions.\n') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Convert a string graph representation to an adjacency matrix
%                                          (see also adj2str.m)
%
% INPUTs: string graph representation: .i1.j1.k1,.i2.j2.k2,....
% OUTPUTs: adjacency matrix, nxn, n - number of nodes
%
% Note 1: Valid for a general graph.
% Note 2: This is the reverse routine for adj2str.m.
% Note 3: The string nomenclature is arbitrarily chosen.
%
% GB: last updated, Sep 25, 2012

function adj = str2adj(str)

commas=find(str==',');
n=length(commas); % number of nodes
adj=zeros(n); % initialize adjacency matrix

if commas(1)>1
  % Extract the neighbors of the first node only
  neigh=str(1:commas(1)-1);
  dots=find(neigh=='.');
  for d=1:length(dots)-1; adj(1,str2num(neigh(dots(d)+1:dots(d+1)-1)))=1; end
  adj(1,str2num(neigh(dots(length(dots))+1:length(neigh))))=1;
end

% Extract the neighbors of the remaining 2:n nodes
for i=2:n
    neigh=str(commas(i-1)+1:commas(i)-1);
    if isempty(neigh); continue; end

    dots=find(neigh=='.');
    for d=1:length(dots)-1; adj(i,str2num(neigh(dots(d)+1:dots(d+1)-1)))=1; end

    adj(i,str2num(neigh(dots(length(dots))+1:length(neigh))))=1;

end

author: aeolianine

repo: octave-networks-toolbox


text:
suite = matlab.unittest.TestSuite.fromFolder(pwd);
suite.run;

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function blas1_s_test09 ( )

%*****************************************************************************80
%
%% TEST09 tests SROT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2006
%
%  Author:
%
%    John Burkardt
%
  n = 6;

  for i = 1 : n
    x(i) = i;
  end

  for i = 1 : n
    y(i) = i * i - 12;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST09\n' );
  fprintf ( 1, '  SROT carries out a Givens rotation.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X and Y\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %12f  %12f\n', i, x(i), y(i) );
  end

  c = 0.5;
  s = sqrt ( 1.0 - c * c );
  [ x, y ] = srot ( n, x, 1, y, 1, c, s );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  SROT ( N, X, 1, Y, 1, %f, %f )\n', c, s );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %12f  %12f\n', i, x(i), y(i) );
  end

  for i = 1 : n
    x(i) = i;
  end

  for i = 1 : n
    y(i) = i * i - 12;
  end

  c = x(1) / sqrt ( x(1) * x(1) + y(1) * y(1) );
  s = y(1) / sqrt ( x(1) * x(1) + y(1) * y(1) );
  [ x, y ] = srot ( n, x, 1, y, 1, c, s );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  SROT ( N, X, 1, Y, 1, %f, %f )\n', c, s );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %12f  %12f\n', i, x(i), y(i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% CONVERT ROTATION ANGLES INTO A QUATERNION
function [q] = GetQuaternionFromEulers(eulers)
% THis function converts euler angle rotations into a
% quaternion via its components. Associated block:
% "Rotation Angles to Quaternions".
eulers = 0.5*eulers;
% Build vector of trigonometric arguements
trigArgs = [sin(eulers);cos(eulers)];
% Assemble Quaternion components
q = zeros(4,1);
q(1) = trigArgs(4)*trigArgs(5)*trigArgs(6) + trigArgs(1)*trigArgs(2)*trigArgs(3);
q(2) = trigArgs(4)*trigArgs(5)*trigArgs(3) - trigArgs(1)*trigArgs(2)*trigArgs(6);
q(3) = trigArgs(4)*trigArgs(2)*trigArgs(6) + trigArgs(1)*trigArgs(5)*trigArgs(3);
q(4) = trigArgs(1)*trigArgs(5)*trigArgs(6) - trigArgs(4)*trigArgs(2)*trigArgs(3);
% Renormalise
q = unit(q);
end

repo: OpenMAS


text:
function H = ldivide( F, G )
%.\   Pointwise DISKFUNV left divide.
%   F.\G if G is a DISKFUNV and F is a double this returns (1/F)*G
%
%   F.\G if G is a double and F is a DISKFUNV this returns G\F, but this
%   does not work if F becomes numerically close to zero.
%
%   F.\G is not allowed if both F and G are DISKFUNV objects.
%
%   F.\G is the same as the command LDIVIDE(F, G)
%
% See also RDIVIDE.

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

% Empty check:
if ( ( isempty(F) ) || ( isempty(G) ) )
    H = diskfunv();
    return
end

if ( isa( F, 'double' ) )                % DOUBLE .\ DISKFUNV
    nG = G.nComponents;
    Gc = G.components;
    H = G;
    for j = 1:nG
        H.components{j} = ldivide(F, Gc{j});
    end
elseif ( isa(G, 'double') )              % DISKFUNV .\ DOUBLE
    nF = F.nComponents;
    Fc = F.nComponents;
    H = F;
    for j = 1:nF
        H.components{j} = ldivide(Fc{j}, G);
    end
else                                     % DISKFUNV.\DISKFUNV
    nF = F.nComponents;
    nG = G.nComponents;
    if ( nF ~= nG )
        error('DISKFUN:DISKFUNV:ldivide:dimensionMismatch', ...
            'DISKFUNV do not have the same number of components.')
    end
    H = F;
    for j = 1 : nF
        H.components{j} = ldivide( F.components{j}, G.components{j} );
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function scale = FindChromaticyScale(M, I)
%
%       scale = FindChromaticyScale(M, I)
%
%
%        Input:
%
%
%        Output:
%
%
%     Copyright (C) 2016 Francesco Banterle
%
%     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/>.
%

l_m = length(M);
l_I = length(I);

if((l_m ~= l_I) || isempty(M) || isempty(I))
    error('FindChromaticyScale: input colors have different color channels.');
end

    function err = residualFunction(p)

        I_c = I .* p;

        I_c_n = I_c / norm(I_c);
        M_n = M / norm(M);

        err = sum((I_c_n - M_n).^2);
    end

    opts = optimset('Display', 'none', 'TolFun', 1e-8, 'TolX', 1e-8);
    scale = fminsearch(@residualFunction, ones(1, l_m), opts);

end

author: banterle


text:
function mColormap = gui_Colormap_BlueYellowRed(nSize)

% If size is not specified, set it to 256
if nargin < 1
  nSize = 256;
end

mColormap = [];

nSteps = floor(nSize * 0.25);
mColormap = [mColormap; gui_Interpolate([0.000000, 0.000000, 1.000000], [0.500000, 0.500000, 0.500000], nSteps)];

nSteps = floor(nSize * 0.25);
mColormap = [mColormap; gui_Interpolate([0.500000, 0.500000, 0.500000], [1.000000, 1.000000, 0.000000], nSteps)];

nSteps = floor(nSize * 0.25);
mColormap = [mColormap; gui_Interpolate([1.000000, 1.000000, 0.000000], [1.000000, 0.500000, 0.000000], nSteps)];

nSteps = floor(nSize * 0.25);
mColormap = [mColormap; gui_Interpolate([1.000000, 0.500000, 0.000000], [1.000000, 0.000000, 0.000000], nSteps)];

author: CelsoReyes


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
    = up(fid, data, varargins, opts)
% This file is an entry for the up strategy.
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
%           = up_start(fid, data, varargins, opts)
% cum_ret: cumulative wealth achived at the end of a period.
% cumprod_ret: cumulative wealth achieved till the end each period.
% daily_ret: daily return achieved by a strategy.
% daily_portfolio: daily portfolios
%
% data: market sequence vectors
% fid: handle for write log file
% varargins: variable parameters
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
%          = up_start(fid, data, {0}, opts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract the parameters
tc = varargins{1};      % transaction cost fee rate

% Run the greedy algorithm
[cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
    = up_run(fid, data, tc, opts);

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
classdef (Abstract) nme_buslink < mp.nm_element %& mp.form_ac

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

%     properties
%     end

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

        function np = np(obj)
            np = 4;     %% this is a 4 port element
        end

        function nz = nz(obj)
            nz = 3;     %% 3 complex non-voltage states per element
        end

        function obj = add_zvars(obj, nm, dm, idx)
            p = idx{1};
            ng = obj.nk;

            if p == 1
                nm.init_indexed_name('zr', 'Plink', {obj.nz});
                nm.init_indexed_name('zi', 'Qlink', {obj.nz});
            end
            nm.add_var('zr', 'Plink', {p}, obj.nk, 0, -Inf, Inf);
            nm.add_var('zi', 'Qlink', {p}, obj.nk, 0, -Inf, Inf);
        end

        function obj = build_params(obj, nm, dm)
            build_params@mp.nm_element(obj, nm, dm);    %% call parent
            I = (dm.base_kva / dm.base_mva / 1000) * speye(obj.nk);
            obj.N = [ repmat(I, 1, obj.nz);
                     -speye(obj.nk * obj.nz) ];
        end

        function [A, b_va, b_vm, Istack_] = voltage_constraints(obj)
            %% form constraint matrices for matching voltages
            nk = obj.nk;

            %% basic constraint matrix for voltage equality (all nodes)
            Istack = repmat(speye(nk), obj.nz, 1);  %% stacked identities
            A = [ Istack -speye(nk * obj.nz) ] * obj.C';
            ang120 = 2*pi/3*ones(nk, 1);

            %% RHS
            b_va = [zeros(nk, 1); ang120; -ang120]; %% angle constraints
            b_vm = zeros(nk*obj.nz, 1);             %% magnitude constraints

            if nargout > 3
                Istack_ = Istack;
            end
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
 function out = make_block(in, nblock)
%function out = make_block(in, nblock)
% construct blocks of sparse matrix for Gtomo2_sparse object
% called by Gblock.m
% Copyright 2002-2-19, Jeff Fessler, The University of Michigan

if nargin < 2, ir_usage, end

out = in; % copy entire object

nb = in.nb;
na = in.na;
G = in.G;

for iblock=1:nblock
	ia = iblock:nblock:na;
	ii = outer_sum(1:nb, (ia-1)*nb);
	out.blocks{iblock} = G(ii(:),:);
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
classdef DualUpdaterFactory < handle

    methods (Access = public, Static)

        function d = create(cParams)
            switch cParams.type
                case 'AlternatingPrimalDual'
                    d = DualUpdater_AugmentedLagrangian(cParams);
                case 'LagrangeMultiplierEstimation'

                case 'NullSpace'
                    d = DualUpdater_NullSpace(cParams);
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
clear all; close all;
[x, y]=meshgrid(-128:2:127, -128:2:127);
D=sqrt((x).^2+(y).^2);
D0=50;
W=30;
n=5;
H=1./(1+(D*W./(D.^2-D0^2)).^(2*n));
H2=1-H;
figure;
mesh(double(H));
% axis tight;

author: HuangCongQing


text:
function isargequalsize(x1,varargin)
%ISARGEQUALSIZE throws an error if not the same size
%
%   Usage: isargequalsize(x1,x2,...)
%
%   Input parameters:
%       x1,x2,...   - matrices to test
%
%   See also: isargequallength

%*****************************************************************************
% 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 ===================================
if nargin<2
    error('%s: you need at least two input arguments.',upper(mfilename));
end

%% ===== Checking for equal size =========================================
for ii = 1:nargin-1
    if ~isequal(size(varargin{ii}),size(x1))
        error('%s and %s have not the same size',inputname(1),inputname(ii+1));
    end
end

author: sfstoolbox

repo: sfs-matlab


text:
function [H, L] = storeR( F, levelF, o, H, L)
%------------------------------------------------------------------------------
%
% This function stores a two-dimensional gridfunction F into H for
% future extraction. Properties and location are stored in the book-
% keeping matrix L.
% This function is a two-dimensional lifting scheme utility.
%
% F = 2D gridfunction of coefficients
%
% levelF = integer designated as the level of F
%
% o = character, should be either 'a' or 'd',
%     describing the type of F:
%     'a' relates to approximation (coefficients) and
%     'd' relates to detail (coefficients)
%     other characters 'n' may be added to this list infuture versions
%
% L = 2D integer array of bookkeeping, its contents are updated by calls
%     to functions as storeQ, storeQ.... and does not require the
%     attention of the user.
%     L is organised as follows:
%       (for the j-th gridfunction F appended to H)
%       L(j, 1) = level of F as it has been assigned by the user.
%       L(j, 2) = 0 means that F corresponds to a rectangular domain,
%                   storeR always assigns this value. Functions like
%                   storeQ, storeQ.... assign values other than 0.
%       L(j, 3) = 0 stands for approximation type coefficients, 'a',
%                 1 stands for detail type coefficients, 'd'.
%       L(j, 4) = (integer) 1st dimension of gridfunction F.
%       L(j, 5) = (integer) 2nd dimension of gridfunction F.
%       L(j, 6) = integer that points to the end of the j-th
%                 addition to the heap.
%     With each call to storeR, storeQ, storeQ.... L obtains additional
%     row(s).
%
% H = 1D array that functions as storage (heap). The coefficients of F
%     are appended to H as a result of calling storeR. Bookkeeping is
%     maintained in L. No attention of the user is required.
%
% See also: retrieveR, storeQ
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: November 12, 2001.
% (c) 1999-2002 Stichting CWI, Amsterdam.
%------------------------------------------------------------------------------
addend = [0 0 0 0 0 0];
if isempty(L)
  if ~isempty(H)
    error(' storeR - books empty but heap is not ')
  else
    nL = 0;
%   mL = 6;
    heaptr = 1;
  end
else
  if isempty(H)
    error(' storeR - books not empty but heap is ')
  else
    [nL, mL] = size(L);
    if mL ~= 6
      error(' storeR - books do not fit ')
    else
      heaptr = L(nL, 6);
      if heaptr < 2
        error(' storeR - bookkeeping error ')
      end
    end
  end
end
%
L = [L; addend];
nL = nL + 1;
%
L(nL, 1) = levelF;
L(nL, 2) = 0;
switch o
    case 'a' , L(nL, 3) = 0;
    case 'd' , L(nL, 3) = 1;
%   case 'n' , L(nL, 3) = 2;
    otherwise
        error(' storeR - unknown type of coefficients ')
end
if isempty(F)
  error(' storeR - addition is empty ')
else
  [nF, mF] = size(F);
end
L(nL, 4) = nF;
L(nL, 5) = mF;
L(nL, 6) = heaptr + nF * mF;
%
H = addtoheap(F, H);
%------------------------------------------------------------------------------


text:
% Created for the blood glucose step response
% by John D. Hedengren
% 02 Sept 2003

global u A

% Steady State Initial Conditions for the States
% Basal values of glucose and insulin conc.
G_ss = 4.5; % mmol/L
X_ss = 15; % mU/L
I_ss = 15; % mU/L

x_ss = [G_ss;I_ss;X_ss];

% Steady State Initial Condition for the Control
u_ss = 16.667; % mU/min

% Steady State for the Disturbance
d_ss = 0; % mmol/L-min

% Final Time (min)
tf = 400;

[t,x] = ode15s('blood_glucose',[0 tf],x_ss);

% Separate out the state values
G = x(:,1);
X = x(:,2);
I = x(:,3);

% Plot the results
figure(1);
plot(t,G);
legend('Glucose');
xlabel('Time (min)');
ylabel('mmol/L');

figure(2);
plot(t,I,t,X);
legend('Plasma Insulin','Remote Compartment Insulin');
xlabel('Time (min)');
ylabel('mU/L');


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtHmxHelmholtz2dDt.m                         |
%|    #    |   VERSION    : 0.50                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal & Martin Averseng             |
%|  ( # )  |   CREATION   : 25.11.2018                                    |
%|  / 0 \  |   LAST MODIF :                                               |
%| ( === ) |   SYNOPSIS   : Solve neumann scatering problem with double   |
%|  `---'  |                layer transpose potential                     |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Parameters
N   = 1e3
tol = 1e-3
typ = 'P1'
gss = 3
X0  = [0 -1 0]
k   = 5

% Boundary mesh
mesh = mshCircle(N,1);

% Radiative mesh
radiat = mshSquare(10*N,[5 5]);

% Mesh representation
figure
plot(mesh)
hold on
plot(radiat)
plotNrm(mesh)
axis equal
axis(2.5*[-1 1 -1 1 -1 1])

% Frequency adjusted to maximum esge size
stp  = mesh.stp;
kmax = 1/stp(2)
% k    = kmax
if (k > kmax)
    warning('Wave number is too high for mesh resolution')
end
f = (k*340)/(2*pi);

% Incident wave
PW = @(X) exp(1i*k*X*X0');
gradxPW{1} = @(X) 1i*k*X0(1) .* PW(X);
gradxPW{2} = @(X) 1i*k*X0(2) .* PW(X);
gradxPW{3} = @(X) 1i*k*X0(3) .* PW(X);

% Incident wave representation
plot(radiat,real(PW(radiat.vtx)))
title('Incident wave')
xlabel('X');   ylabel('Y');   zlabel('Z');
hold off
alpha(0.80)
view(0,90)

%%% PREPARE OPERATOR
disp('~~~~~~~~~~~~~ PREPARE OPERATOR ~~~~~~~~~~~~~')

% Green kernels
dyGxy{1} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]1',k);
dyGxy{2} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]2',k);
dyGxy{3} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]3',k);

% Domain
sigma = dom(mesh,gss);

% Finite elements
u = fem(mesh,typ);
v = fem(mesh,typ);

% Mass matrix
Id = integral(sigma,u,v);

% Finite element boundary operator --> \int_Sx \int_Sy psi(x)' dnyG(x,y) psi(y) dx dy
tic
D = (1i/4) .* integral(sigma,sigma,u,dyGxy,ntimes(v),tol);
toc

% Regularization
tic
Dr = -1/(2*pi) .* regularize(sigma,sigma,u,'grady[log(r)]',ntimes(v));
toc

% Operator [-Id/2 + Dt]
LHS = - 0.5*Id + (D + Dr).';

% Finite element incident wave trace --> \int_Sx psi(x) dnx(pw(x)) dx
RHS = - integral(sigma,ntimes(u),gradxPW);

% Structure
figure
subplot(2,2,1:2)
spy(LHS)

%%% SOLVE LINEAR PROBLEM
disp('~~~~~~~~~~~~~ SOLVE LINEAR PROBLEM ~~~~~~~~~~~~~')

% LU factorization
tic
[Lh,Uh] = lu(LHS);
toc
subplot(2,2,3)
spy(Lh)
subplot(2,2,4)
spy(Uh)

% Solve linear system [-Id/2 + Dt] * lambda = dnP0
tic
lambda = Uh \ (Lh \ RHS); % LHS \ RHS;
toc

%%% INFINITE SOLUTION
disp('~~~~~~~~~~~~~ INFINITE RADIATION ~~~~~~~~~~~~~')

% Plane waves direction
theta = 2*pi/1e3 .* (1:1e3)';
nu    = [sin(theta),cos(theta),zeros(size(theta))];

% Green kernel function
xdoty = @(X,Y) X(:,1).*Y(:,1) + X(:,2).*Y(:,2);
Ginf  = @(X,Y) exp(-1i*k*xdoty(X,Y));

% Finite element infinite operator --> \int_Sy exp(ik*nu.y) * psi(y) dx
Sinf = 1i/4 * integral(nu,sigma,Ginf,v);

% Finite element radiation
sol = Sinf * lambda;

% Analytical solution
ref = diskHelmholtz('inf','neu',1,k,nu);
norm(ref-sol,2)/norm(ref,2)
norm(ref-sol,'inf')/norm(ref,'inf')

% Graphical representation
figure
plot(theta,log(abs(sol)),'b',theta,log(abs(ref)),'--r')

%%% DOMAIN SOLUTION
disp('~~~~~~~~~~~~~ RADIATION ~~~~~~~~~~~~~')

% Green kernels
Gxy = @(X,Y) femGreenKernel(X,Y,'[H0(kr)]',k);

% Finite element radiative operator --> \int_Sy G(x,y) psi(y) dy
tic
Sdom = 1i/4 .* integral(radiat.vtx,sigma,Gxy,v,tol);
toc

% Regularization
tic
Sreg = -1/(2*pi) .* regularize(radiat.vtx,sigma,'[log(r)]',v);
Sdom = Sdom + Sreg;
toc

% Domain solution
Psca = Sdom * lambda;
Pinc = PW(radiat.vtx);
Pdom = Psca + Pinc;

% Annulation mesh interieure
r             = sqrt(sum(radiat.vtx.^2,2));
Pdom(r<=1.01) = Pinc(r<=1.01);

% Graphical representation
figure
plot(radiat,abs(Pdom))
axis equal
title('Total field solution')
colorbar

%%% ANAYTICAL SOLUTIONS FOR COMPARISONS
% Analytical solution
Pdom = diskHelmholtz('dom','neu',1,k,radiat.vtx) + PW(radiat.vtx);

% Solution representation
figure
plot(radiat,abs(Pdom))
axis equal;
title('Analytical solution')
colorbar
view(0,90)

disp('~~> Michto gypsilab !')


text:
function cc_levels_constrained_display ( )

%*****************************************************************************80
%
%% CC_LEVELS_CONSTRAINED_DISPLAY displays grids generated by CC_LEVELS_CONSTRAINED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 November 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CC_LEVELS_CONSTRAINED_DISPLAY:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Display the 2D Clenshaw-Curtis grids\n' );
  fprintf ( 1, '  generated by CC_LEVELS_CONSTRAINED.\n' );

  dim_num = 2;

  while ( 1 )
%
%  Get user input.
%
    q_max = input ( 'Enter Q_MAX or RETURN to exit;' );

    if ( isempty ( q_max ) )
      break
    end

    if ( q_max < dim_num )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  We require DIM_NUM <= Q_MAX!\n' );
      continue
    end

    alpha = input ( 'Enter [ ALPHA1, ALPHA2 ] or RETURN to exit;' );

    if ( isempty ( alpha ) )
      break
    end

    level_min = input ( 'Enter [ LEVEL_MIN1, LEVEL_MIN2 ] or RETURN to exit;' );

    if ( isempty ( level_min ) )
      break
    end

    level_max = input ( 'Enter [ LEVEL_MAX1, LEVEL_MAX2 ] or RETURN to exit;' );

    if ( isempty ( level_max ) )
      break
    end
%
%  Compute data.
%
    [ grid_num, point_num ] = cc_levels_constrained_size ( dim_num, ...
      q_max, alpha, level_min, level_max );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Number of grids is %d\n', grid_num );
    fprintf ( 1, '  Number of points is %d\n', point_num );

    [ grid_level, grid_point ] = cc_levels_constrained ( dim_num, ...
      q_max, alpha, level_min, level_max, grid_num, point_num );

    clf
%
%  We have to name the axes in order to control the grid.
%
    axes_handle = axes;
%
%  Plot the points.
%
    handle = scatter ( grid_point(1,:), grid_point(2,:), 'filled' );
%
%  Force the plotting region to be square, not rectangular.
%
    axis square
%
%  Request grid lines.
%
    grid on
%
%  Specify the location of the grid lines, and suppress labeling.
%
    set ( axes_handle, 'xtick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
    set ( axes_handle, 'xticklabel', [] );
    set ( axes_handle, 'ytick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
    set ( axes_handle, 'yticklabel', [] );
%
%  Make the plotting region slightly bigger than the data.
%
    axis ( [ -1.1, 1.1, -1.1, 1.1 ] )
%
%  Title
%
    s = sprintf ( 'ALPHA = [ %f, %f]', alpha(1), alpha(2) );
    title ( s );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CC_LEVELS_CONSTRAINED_DISPLAY:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

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

  return
ed


text:
function A = double(X)
%DOUBLE Convert tensor to double array.
%
%   A = double(X) converts X to a standard multidimensional array.
%
%   See also TENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

A = double(X.data);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%  Figure 10.72      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% Fig. 10.72
% Data for RTP Demo 3-3-99
% Data provided by Dr. Gwen van der Linden
% Data is from System Identification Studies
InputFlux=[3.460064464376177e-1 1.177299050104922e-1 2.838023866104041e-2;
   3.880303397347619e-11 8.024902450324316e-2 1.807231516460469e-2;
   8.004191616976514e-9 2.721604310757543e-3 3.171348842079633e-2];
M_inv=diag([1.000040130716728 5.557442686788876 13.63821806414694]);
Radiation=[5.47621193859299e-2 -8.570695054070524e-3 -8.296135532988507e-4...
      -4.536181077856052e-2;
   -8.570695054070524e-3 8.570946319867835e-3 -1.621311365067015e-7...
      -8.913466080455817e-8;
   -8.296135532988507e-4 -1.621311365067015e-7 8.299854517643017e-4...
      -2.097673289443245e-7];
Conduction=[3.559939609150268e-7 -1.113667477845243e-7 -1.976161155515125e-7...
      -4.701109757899004e-8;
   -1.113667477845243e-7 1.160207476868843e-2 -2.502736022145532e-3...
      -9.099227379795117e-3;
   -1.976161155515125e-7 -2.502736022145532e-3 6.37364815665867e-3...
      -3.870714518397587e-3];
ScaleTemp=diag([0.01 0.01 0.01 0.01]);

clf;
%
sim('fig10_71')
%plot(tout,r,'-');
%hold on;
%plot(tout,y,'--');
%xlabel('Time (sec)');
%ylabel('Temperature (K)');
%hold off;
%pause;
ii=240:876;
plot(tout(ii),r(ii),'-',tout(ii),y(ii,2),'--');
legend('r','y');
grid on;
xlabel('Time (sec)');
ylabel('Temperature (K)');
title('Fig. 10.72(a) Temperature tracking response');
pause;
hold off;
ii=240:876;
plot(tout(ii),u(ii),'-');
xlabel('Time (sec)');
ylabel('Lamp voltage (V)');
grid on;
legend('u');
title('Fig. 10.72(b) Control effort');

author: Sable


text:
%%%%%script to generate figure 7 from the paper 'Invariant Scattering Convolution Networks'
%%%%%feb 2012%%%%%%%
%%%%%%Joan Bruna and Stephane Mallat$$$$$$$$$

%%%%%This figure displays the scattering coefficients of
%%%%%% a digit taken from the MNIST dataset.

Nim=32;
copts.renorm_process=0;
copts.l2_renorm=1;
foptions.J=3;
foptions.L=8;
soptions.M=2;
soptions.oversampling = 0;
[Wop,filters]=wavelet_factory_2d([Nim Nim], foptions, soptions);

dirac=zeros(Nim);
dirac(1)=1;
dirac=fftshift(dirac);
[scdirac]=scat(dirac,Wop);

tmp=load('mnist_sample.mat');
im=tmp.im;

%compute scattering coefficients
[sc_digit]=scat(im,Wop);

sc=format_scat(sc_digit);
S=size(sc);

%construct the grid with the scattering displays

drawing1=[];
drawing2=[];
for s1=1:S(2)
	row1=[];
	row2=[];
	for s2=1:S(3)
		tempo=scat_display(sc_digit,scdirac,copts,sc(:,s1,s2));
		row1=[row1 tempo{1}];
		row2=[row2 tempo{2}];
	end
	drawing1=[drawing1 ; row1];
	drawing2=[drawing2 ; row2];
end

[N,M]=size(drawing1);

digit_bg=imresize(im,[N M]);
digit_bg=min(255,max(0,digit_bg));

atten=0.0;

fact1=atten* max(drawing1(:)) / max(digit_bg(:));
fact2=atten* max(drawing2(:)) / max(digit_bg(:));

gg=colormap(gray);
invgg=1-gg;
figure
imagesc(drawing1)
set(gca,'XTick',[]);
set(gca,'YTick',[]);
axis square

figure
imagesc(drawing2)
set(gca,'XTick',[]);
set(gca,'YTick',[]);
axis square

figure
imagesc(digit_bg);
colormap(invgg);
set(gca,'XTick',[]);
set(gca,'YTick',[]);
axis square


text:
function lengths = trk_length(tracks)
%TRK_LENGTH - Calculate the lengths of tracks
%
% Syntax: lengths = trk_length(tracks)
%
% Inputs:
%    tracks - TrackVis track group. For tracks all of the same length (i.e.
%    TRK_INTERP has been run), either structure or matrix form is fine.
%
% Outputs:
%    lengths [1 x nTracks]
%
% Example:
%    [header tracks] = read_trk(trkPath);
%    lengths         = trk_length(tracks);
%    mean(lengths), std(lengths)
%
% Other m-files required: trk_restruc
% Subfunctions: none
% MAT-files required: none
%
% See also: TRK_READ

% Author: John Colby (johncolby@ucla.edu)
% UCLA Developmental Cognitive Neuroimaging Group (Sowell Lab)
% Apr 2010

% Put in matrix form if possible
if isstruct(tracks) && length(unique(cat(tracks.nPoints)))==1
    tracks = trk_restruc(tracks);
end

% Fast matrix operation if all tracks are the same length
if isnumeric(tracks)
    lengths = sum(sqrt(squeeze(sum((tracks(2:end,1:3,:) - tracks(1:(end-1),1:3,:)).^2, 2))), 1);

% Slow forloop if tracks are not the same length
else
    lengths = zeros(1,length(tracks));
    for i=1:length(tracks)
        lengths(i) = sum(sqrt(squeeze(sum((tracks(i).matrix(2:end,1:3) - tracks(i).matrix(1:(end-1),1:3)).^2, 2))), 1);
    end
end

author: yetianmed

repo: subcortex

save_path: github-repos/MATLAB/yetianmed-subcortex


text:
function [pb,theta] = pb2MM2(im,sigma)
% function [pb,theta] = pb2MM2(im,sigma)
%
% Compute probability of boundary using the spatially averaged
% second moment matrix at multiple scales.
%
% See also det2MM.
%
% David R. Martin <dmartin@eecs.berkeley.edu>
% March 2003

if nargin<2, sigma=1; end

switch sigma, % from logistic fits (train2MM2.m)
 case 1, beta = [ -3.2369080e+00 ...
                  -7.9668057e+00 -7.1098318e-01 ...
                  -4.6178498e-01  3.1575066e+01 ];
 case 2, beta = [ -3.5512404e+00 ...
                  -8.0101784e+00  1.3527093e+01 ...
                   1.0531737e+01  1.9292447e+01 ];
 otherwise,
  error(sprintf('no parameters for sigma=%g\n',sigma));
end

h = size(im,1);
w = size(im,2);
[a1,b1,t1] = det2MM(im,sigma);
[a2,b2,t2] = det2MM(im,sigma*2);
a1 = sqrt(a1(:));
b1 = sqrt(b1(:));
a2 = sqrt(a2(:));
b2 = sqrt(b2(:));
x = [ ones(size(a1)) a1 b1 a2 b2 ];
pb = 1 ./ (1 + exp(-x*beta'));
pb = reshape(pb,[h w]);

% average orientations for nonmax suppression
dt = mod(t2-t1,2*pi); % [0,2pi)
dt = (dt<pi).*dt + (dt>=pi).*(dt-2*pi); % [-pi,pi)
pb = nonmax(pb,t1+dt/2);

author: Cloud-CV

repo: object-proposals


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

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [codebook,scheme] = code_OneVsAll(m)
% Generate the codebook for multiclass classification with One-Versus-All encoding.
%
% codebook = code_OneVsAll(m)
%
%  see also:
%    code, codedist_hamming

% (c) SCD-KULeuven, rights & help @ http://www.esat.kuleuven.be/sista/lssvmlab

codebook = eye(m).*2-1;
scheme = []; for i=1:m, scheme = [scheme 'b']; end

author: peterhcharlton

repo: RRest


text:
function g = gama_dct(I)

import bliinds2.*;

temp1=dct2(I);
temp2=temp1(:);
temp3=temp2(2:end);

%g=kurtosis(temp3);
g=gama_gen_gauss(temp3);

author: dsoellinger


text:
function [pyFiltS,cerrFiltS] = comparePyradFilters
% Function to compare filtered images computed using CERR vs. Pyradiomics.
% Supported filer types: LoG, Wavelets.
%------------------------------------------------------------------------
% AI 07/02/2020

%% Load sample data
fpath = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing/data_for_cerr_tests/IBSI1_CT_phantom/IBSILungCancerCTImage.mat.bz2');
planC = loadPlanC(fpath,tempdir);
planC = updatePlanFields(planC);
planC = quality_assure_planC(fpath,planC);
indexS = planC{end};
scanNum=1;

%% Apply filters using pyradiomics

% 1. LoG
filtName = 'LoG';
testSigmaV = [3.0,4.0];
pyFiltParam1S.sigma = testSigmaV;
pyFiltS = processImageUsingPyradiomics(planC,'',filtName,pyFiltParam1S);
pyFieldLoGC = fieldnames(pyFiltS);

%2. Wavelets
filtName = 'wavelet';
pyFiltParam2S.wavetype = 'coif1';
pyFilt2S = processImageUsingPyradiomics(planC,'',filtName,pyFiltParam2S);
pyFieldWavC = fieldnames(pyFilt2S);
for n = 1:length(pyFieldWavC)
    pyFiltS.(pyFieldWavC{n}) = pyFilt2S.(pyFieldWavC{n});
end

%% Apply filters using CERR

%Get scan array
scan3M = getScanArray(scanNum,planC);
CToffset = planC{indexS.scan}(scanNum).scanInfo(1).CTOffset;
scan3M = double(scan3M)-CToffset;
wholeScanMask3M = ones(size(scan3M));

% 1. LoG
filtName = 'LoG';
cerrFiltS = struct();
scanS = planC{indexS.scan}(scanNum);
[xV,yV,zV] = getScanXYZVals(scanS);
dx = median(abs(diff(xV)));
dy = median(abs(diff(yV)));
dz = median(diff(zV));
voxelSizeV = [dx,dy,dz];
cerrFiltParam1S.VoxelSize_mm.val = voxelSizeV*10;%convert to mm
for n = 1:length(testSigmaV)
        cerrFiltParam1S.Sigma_mm.val = testSigmaV(n);
        outFieldname = pyFieldLoGC{n};
        outS = processImage(filtName,scan3M,wholeScanMask3M,cerrFiltParam1S);
        fieldsC = fieldnames(outS);
        cerrFiltS.(outFieldname) = outS.(fieldsC{1});
end

% 2. Wavelets
cerrFiltParam2S = struct();
cerrFiltParam2S.Wavelets.val='coif';
cerrFiltParam2S.Index.val='1';
cerrFiltParam2S.NormFlag.val=0;
dirC = {'LLH','LHL','LHH','HLL','HLH','HHL','HHH','LLL'};
filtName = 'Wavelets';
for n = 1:length(dirC)
    %Convert python dictionary to matlab struct
    cerrFiltParam2S.Direction.val = dirC{n};
    outS = processImage(filtName,scan3M,wholeScanMask3M,cerrFiltParam2S);
    outFieldname = pyFieldWavC{n};
    fieldsC = fieldnames(outS);
    cerrFiltS.(outFieldname) = outS.(fieldsC{1});
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function isargmatrix(varargin)
%ISARGMATRIX throws an error if not a matrix
%
%   Usage: isargmatrix(arg1,arg2,...)
%
%   Input parameters:
%       args        - list of args
%
%   See also: isargvector, isargscalar

%*****************************************************************************
% 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 for matrix =============================================
for ii = 1:nargin
    if ~isnumeric(varargin{ii}) || ndims(varargin{ii})~=2
        error('%s need to be a matrix.',inputname(ii));
    end
end


text:
function order = level_to_order_open ( dim_num, level )

%*****************************************************************************80
%
%% LEVEL_TO_ORDER converts a level to an order for open rules.
%
%  Discussion:
%
%    Sparse grids can naturally be nested.  A natural scheme is to use
%    a series of one-dimensional rules arranged in a series of "levels"
%    whose order roughly doubles with each step.
%
%    The arrangement described here works naturally for the Fejer Type 1,
%    Fejer Type 2, Newton Cotes Open, Newton Cotes Half Open,
%    and Gauss-Patterson rules.  It also can be used, partially, to describe
%    the growth of Gauss-Legendre rules.
%
%    The idea is that we start with LEVEL = 0, ORDER = 1 indicating the single
%    point at the center, and for all values afterwards, we use the relationship
%
%      ORDER = 2**(LEVEL+1) - 1.
%
%    The following table shows how the growth will occur:
%
%    Level    Order
%
%    0          1
%    1          3 =  4 - 1
%    2          7 =  8 - 1
%    3         15 = 16 - 1
%    4         31 = 32 - 1
%    5         63 = 64 - 1
%
%    For the Fejer Type 1, Fejer Type 2, Newton Cotes Open,
%    Newton Cotes Open Half, and Gauss-Patterson rules, the point growth is
%    nested.  If we have ORDER points on a particular LEVEL, the next level
%    includes all these old points, plus ORDER+1 new points, formed in the
%    gaps between successive pairs of old points plus an extra point at each
%    end.
%
%    Level    Order = New + Old
%
%    0          1   =  1  +  0
%    1          3   =  2  +  1
%    2          7   =  4  +  3
%    3         15   =  8  +  7
%    4         31   = 16  + 15
%    5         63   = 32  + 31
%
%    If we use a series of Gauss-Legendre rules, then there is almost no
%    nesting, except that the central point is shared.  If we insist on
%    producing a comparable series of such points, then the "nesting" behavior
%    is as follows:
%
%    Level    Order = New + Old
%
%    0          1   =  1  +  0
%    1          3   =  2  +  1
%    2          7   =  6  +  1
%    3         15   = 14  +  1
%    4         31   = 30  +  1
%    5         63   = 62  +  1
%
%    Moreover, if we consider ALL the points used in such a set of "nested"
%    Gauss-Legendre rules, then we must sum the "NEW" column, and we see that
%    we get roughly twice as many points as for the truly nested rules.
%
%    In this routine, we assume that a vector of levels is given,
%    and the corresponding orders are desired.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL(DIM_NUM), the nesting level.
%
%    Output, integer ORDER(DIM_NUM,1), the order (number of points) of the rule.
%
  order = zeros ( dim_num, 1 );

  for dim = 1 : dim_num

    if ( level(dim) < 0 )
      order(dim,1) = -1;
    elseif ( level(dim) == 0 )
      order(dim,1) = 1;
    else
      order(dim,1) = 2^( level(dim) + 1 ) - 1;
    end

  end

  return
end

author: johannesgerer


text:
function w = dwt2D(x, J, af)

% discrete 2-D wavelet transform
%
% USAGE:
%   w = dwt2D(x, stages, af)
% INPUT:
%   x - N by M matrix
%       1) M, N are both even
%       2) min(M,N) >= 2^(J-1)*length(af)
%   J - number of stages
%   af - analysis filters
% OUPUT:
%   w - cell array of wavelet coefficients
% EXAMPLE:
%   [af, sf] = farras;
%   x = rand(128,64);
%   J = 3;
%   w = dwt2D(x,J,af);
%   y = idwt2D(w,J,sf);
%   err = x - y;
%   max(max(abs(err)))
%
% WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
% http://taco.poly.edu/WaveletSoftware/

for k = 1:J
    [x w{k}] = afb2D(x, af, af);
end
w{J+1} = x;

author: thomaskuestner


text:
function [X_poly] = polyFeatures(X, p)
%POLYFEATURES Maps X (1D vector) into the p-th power
%   [X_poly] = POLYFEATURES(X, p) takes a data matrix X (size m x 1) and
%   maps each example into its polynomial features where
%   X_poly(i, :) = [X(i) X(i).^2 X(i).^3 ...  X(i).^p];
%

% You need to return the following variables correctly.
X_poly = zeros(numel(X), p);

% ====================== YOUR CODE HERE ======================
% Instructions: Given a vector X, return a matrix X_poly where the p-th
%               column of X contains the values of X to the p-th power.
%
%

for i=1:p
	X_poly(:, i) = X(:,1).^i;
end

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

end

author: vugsus

repo: coursera-machine-learning


classdef empty2dLayer < nnet.layer.Layer
    % 主要想用于导入yolov3中route层的单个跳层连接
    % 2019.9.8
    % 2020.4.29加入属性connectID，用于连接前面的层
    properties
        connectID %connectID 是以cfg文件中第一个非[net]开始的module为0开始的计数
    end

    methods
        function layer = empty2dLayer(name,con)
            layer.Name = name;
            text = ['[', num2str(1), ' ', num2str(1), '] emptyLayer '];
            layer.Description = text;
            layer.Type = 'empty2dLayer';
            layer.connectID= con;
        end

        function Z = predict(layer, X)
               Z = X;
        end

        function [dX] = backward( layer, X, ~, dZ, ~ )
            dX = dZ;
        end
    end
end
%%
% 参考：https://ww2.mathworks.cn/matlabcentral/fileexchange/71277-deep-learning-darknet-importer


text:
%% Analyzing Neural Time Series Data
% Matlab code for Chapter 16
% Mike X Cohen
%
% This code accompanies the book, titled "Analyzing Neural Time Series Data"
% (MIT Press). Using the code without following the book may lead to confusion,
% incorrect data analyses, and misinterpretations of results.
% Mike X Cohen assumes no responsibility for inappropriate or incorrect use of this code.

%% Figure 16.1

load sampleEEGdata

channel2plot = 'o1';
timewin      = 400; % in ms

timewinidx = round(timewin/(1000/EEG.srate));
tapers     = dpss(timewinidx,5); % this line will crash without matlab signal processing toolbox

% extract a bit of EEG data
d = detrend(squeeze(EEG.data(strcmpi(channel2plot,{EEG.chanlocs.labels}),200:200+timewinidx-1,10)));

% plot EEG data snippet
figure
subplot(5,2,1)
plot(d)
axis tight,axis off

% plot tapers
for i=1:5
    subplot(5,2,(2*(i-1))+2)
    plot(tapers(:,i))
    axis tight,axis off
end

% plot taper.*data
figure
for i=1:5
    subplot(5,2,(2*(i-1))+1)
    plot(tapers(:,i).*d')
    axis tight,axis off
end

% plot fft of taper.*data
f=zeros(5,timewinidx);
for i=1:5
    subplot(5,2,(2*(i-1))+2)
    f(i,:)=fft(tapers(:,i).*d');
    plot(abs(f(i,1:timewinidx/2)).^2)
    axis tight,axis off
end

figure
subplot(5,2,2)
plot(mean(abs(f(:,1:timewinidx/2)).^2,1))
axis tight, axis off

subplot(5,2,3)
hann = .5*(1-cos(2*pi*(1:timewinidx)/(timewinidx-1)));
plot(hann)
axis tight, axis off

subplot(525)
plot(hann.*d)
axis tight, axis off

subplot(526)
ff=fft(hann.*d);
plot(mean(abs(ff(1:timewinidx/2)).^2,1))
axis tight, axis off

%% Figure 16.2

channel2plot    = 'p7';
frequency2plot  = 15;  % in Hz
timepoint2plot  = 200; % ms

nw_product      = 3;  % determines the frequency smoothing, given a specified time window
times2save      = -300:50:1000;
baseline_range  = [-200 -00];
timewin         = 400; % in ms

% convert time points to indices
times2saveidx = dsearchn(EEG.times',times2save');
timewinidx    = round(timewin/(1000/EEG.srate));

% find baselinetimepoints
baseidx = zeros(size(baseline_range));
[~,baseidx(1)] = min(abs(times2save-baseline_range(1)));
[~,baseidx(2)] = min(abs(times2save-baseline_range(2)));
% note that the following line is equivalent to the previous three
%baseidx = dsearchn(times2save',baseline_range');

% define tapers
tapers = dpss(timewinidx,nw_product); % note that in practice, you'll want to set the temporal resolution to be a function of frequency
% define frequencies for FFT
f = linspace(0,EEG.srate/2,floor(timewinidx/2)+1);

% find logical channel index
chanidx = strcmpi(channel2plot,{EEG.chanlocs.labels});

% initialize output matrix
multitaper_tf = zeros(floor(timewinidx/2)+1,length(times2save));

% loop through time bins
for ti=1:length(times2saveidx)

    % initialize power vector (over tapers)
    taperpow = zeros(floor(timewinidx/2)+1,1);

    % loop through tapers
    for tapi = 1:size(tapers,2)-1

        % window and taper data, and get power spectrum
        data      = bsxfun(@times,squeeze(EEG.data(chanidx,times2saveidx(ti)-floor(timewinidx/2)+1:times2saveidx(ti)+ceil(timewinidx/2),:)),tapers(:,tapi));
        pow       = fft(data,timewinidx)/timewinidx;
        pow       = pow(1:floor(timewinidx/2)+1,:);
        taperpow  = taperpow + mean(pow.*conj(pow),2);
    end

    % finally, get power from closest frequency
    multitaper_tf(:,ti) = taperpow/tapi;
end

% db-correct
db_multitaper_tf = 10*log10( multitaper_tf ./ repmat(mean(multitaper_tf(:,baseidx(1):baseidx(2)),2),1,length(times2save)) );

% plot time courses at one frequency band
figure
subplot(121)
[junk,freq2plotidx]=min(abs(f-frequency2plot)); % can replace "junk" with "~"
plot(times2save,mean(log10(multitaper_tf(freq2plotidx-2:freq2plotidx+2,:)),1))
title([ 'Sensor ' channel2plot ', ' num2str(frequency2plot) ' Hz' ])
axis square
set(gca,'xlim',[times2save(1) times2save(end)])

subplot(122)
[junk,time2plotidx]=min(abs(times2save-timepoint2plot));
plot(f,log10(multitaper_tf(:,time2plotidx)))
title([ 'Sensor ' channel2plot ', ' num2str(timepoint2plot) ' ms' ])
axis square
set(gca,'xlim',[f(1) 40])

% plot full TF map
figure
contourf(times2save,f,db_multitaper_tf,40,'linecolor','none')
set(gca,'clim',[-2 2])
xlabel('Time (ms)'), ylabel('Frequency (Hz)')
title([ 'Power via multitaper from channel ' channel2plot ])

%% end.

author: mikexcohen

repo: AnalyzingNeuralTimeSeries


text:
function feat = spp_poolX_to_fcX(feat, layer, spp_model, use_gpu)
% feat = spp_poolX_to_fcX(feat, layer, spp_model, use_gpu)
%   On-the-fly conversion of last pool features to fcX
%   using the weights and biases stored in spp_model.cnn.layers.
%   feat is transformed in columns for continuous memory access and fast
%   speed
%
% Adapted from spp code written by Ross Girshick
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Shaoqing Ren
%
% This file is part of the SPP code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

fc_layer_ind = layer - (spp_model.cnn.first_fc_idx - 1);

if fc_layer_ind <= 0
    return;
end

if use_gpu
    maxNumel = 1024*1024*256; % around 4G memory
    if numel(feat) < maxNumel
        feat_gpu = gpuArray(feat);
        for i = 1:fc_layer_ind
            % weights{1} = matrix of CNN weights [input_dim x output_dim]
            % weights{2} = column vector of biases
            feat_gpu = max(0, bsxfun(@plus, spp_model.cnn.layers(i).weights_gpu{1} * feat_gpu, ...
                              spp_model.cnn.layers(i).weights_gpu{2}));
        end
        feat = gather(feat_gpu);
    else
        nSampleEach = floor(maxNumel / size(feat, 1));
        nSplits = ceil(size(feat, 2) / nSampleEach);
        splits = ones(nSplits, 1) * nSampleEach;
        splits(end) = size(feat, 2) - sum(splits(1:end-1));
        assert(sum(splits) == size(feat, 2));
        feats = mat2cell(feat, size(feat, 1), splits);
        for is = 1:length(feats)
            feat_gpu = gpuArray(feats{is});
            for i = 1:fc_layer_ind
                % weights{1} = matrix of CNN weights [input_dim x output_dim]
                % weights{2} = column vector of biases
                feat_gpu = max(0, bsxfun(@plus, spp_model.cnn.layers(i).weights_gpu{1} * feat_gpu, ...
                                  spp_model.cnn.layers(i).weights_gpu{2}));
            end
            feats{is} = gather(feat_gpu);
        end
        feat = cell2mat(feats);
    end

else
    for i = 1:fc_layer_ind
        % weights{1} = matrix of CNN weights [input_dim x output_dim]
        % weights{2} = column vector of biases
        feat = max(0, bsxfun(@plus, spp_model.cnn.layers(i).weights{1} * feat, ...
                              spp_model.cnn.layers(i).weights{2}));
    end
end

author: ShaoqingRen


text:
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc % clean the command prompt
L = 0.1; % 10 cm length
W = 0.01; % 1 cm thickness
Nx = 50; % number of cell in x direction
Ny = 20; % number of cells in the y direction
m = createMesh2D(Nx, Ny, L, W); % creates a 2D Cartesian grid
% m = createMesh3D(Nx, Ny,Ny, L, W,W); % creates a 3D Cartesian grid, activate this line if you are curious
T_inf = 25+273.15; % [K] ambient temperature
T_base = 100+273.15; % [K] temperature at the base of the fin
k_val = 237; % W/(m.K) thermal conductivity
h_val = 10; % W/(m^2.K) heat transfer coefficient
k = createCellVariable(m, k_val); % assign thermal cond. value to all cells
k_face = geometricMean(k); % geometric average of the thermal conductivity values on the cell faces
BC = createBC(m); % creates a BC structure for the domain m; all Neumann boundaries
BC.left.a(:)=0; BC.left.b(:)=1; BC.left.c(:)=T_base; % convert the left boundary to constant temperature
BC.right.a(:)=k_val/h_val; BC.right.b(:)=1; BC.right.c(:)=T_inf; % right boundary to Robin
BC.top.a(:)=k_val/h_val; BC.top.b(:)=1; BC.top.c(:)=T_inf; % top boundary to Robin
BC.bottom.a(:)=-k_val/h_val; BC.bottom.b(:)=1; BC.bottom.c(:)=T_inf; % bottom boundary to Robin (don't forget the normal vector sign)
M_cond = diffusionTerm(k_face); % matrix of coefficients for the heat diffusion
[M_bc, RHS_bc] = boundaryCondition(BC); % matrix of coefficients and RHS vector for the boundary conditions
T = solvePDE(m, M_cond+M_bc, RHS_bc); % solve the linear system of discretized PDE
visualizeCells(T); % visualize the results

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function d = distmat3(x, y)
%DISTMAT3 Distance matrix (three nested for-loops).
%
%   D = DISTMAT3(X, Y) returns the distance matrix with all distances
%   between the points represented by the rows of X and Y.
%
%   DISTMAT3(X) is equivalent to DISTMAT3(X, X), but the former computes the
%   distance matrix faster.
%
%   Distance is Euclidean.
%
%   The calculation is done with three for-loops.
%
%   See also DISTMAT0, DISTMAT1, DISTMAT2.

%   Author:      Peter J. Acklam
%   Time-stamp:  2002-03-03 13:51:22 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   error(nargchk(1, 2, nargin));

   if nargin == 1               % DISTMAT3(X)

      if ndims(x) ~= 2
         error('Input must be a matrix.');
      end

      [m, p] = size(x);
      d = zeros(m, m);          % initialise output matrix

      for i = 1:m-1
         for j = i+1:m
            ssq = 0;            % sum of squares
            for k = 1:p
               ssq = ssq + abs(x(i,k) - x(j,k))^2;
            end
            d(i,j) = sqrt(ssq);
            d(j,i) = d(i,j);
         end
      end

   else                         % DISTMAT3(X)

      if ndims(x) ~= 2 | ndims(y) ~= 2
         error('Input must be two matrices.');
      end

      [mx, nx] = size(x);
      [my, ny] = size(y);

      if nx ~= ny
         error('Both matrices must have the same number of columns.');
      end

      m = mx;                   % number of rows in distance matrix
      n = my;                   % number of columns in distance matrix
      p = nx;                   % dimension of each point
      d = zeros(m, n);          % initialise output matrix

      for i = 1:m
         for j = 1:n
            ssq = 0;         % Sum of squares.
            for k = 1:p
               ssq = ssq + abs(x(i,k) - y(j,k))^2;
            end
            d(i,j) = sqrt(ssq);
         end
      end

   end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
% update the elements of a sparse matrix .
% sparse_update(S, v, length(v));
% Note: the following must be hold or the MATLAB crashes
%     nnz(S) == length(v).
% TODO: write dimensionality check in c.

repo: lrslibrary


text:
function Fderandomized = derandomize(F)

chanceDeclarations = find(is(F,'chance'));
randomDeclarations = find(is(F,'random'));
if isempty(randomDeclarations)
    error('Cannot derandomize, cannot find declaration of random variables');
end

keep = ones(length(F),1);
keep(randomDeclarations)=0;
keep(chanceDeclarations)=0;
randomVariables = extractRandomDefinitions(F(randomDeclarations));
groupedChanceConstraints = groupchanceconstraints(F);

Fderandomized = deriveChanceModel(groupedChanceConstraints,randomVariables)

Fderandomized = Fderandomized + F(find(keep));

function Fderandomized = deriveChanceModel(groupedChanceConstraints,randomVariables);

Fderandomized = [];
for ic = 1:length(groupedChanceConstraints)
    if length(groupedChanceConstraints{ic})>1
        error('Joint chance still not supported');
    end
    if ~is(groupedChanceConstraints{ic},'elementwise')
        error('Only elementwise chance constraints supported')
    end
    X = sdpvar(groupedChanceConstraints{ic});
    if length(X)>1
        error('Only single elementwise chance constraints supported')
    end

    % OK, simple linear inequality
    allVars = depends(X);
    wVars = [];for j = 1:length(randomVariables);wVars = [wVars getvariables(randomVariables{j}.variables)];end
    xVars = setdiff(allVars,wVars);
    x = recover(xVars);
    w = recover(wVars);

    b = [];
    A = [];
    % Some pre-calc
    xw = [x;w];
    xind = find(ismembc(getvariables(xw),getvariables(x)));
    wind = find(ismembc(getvariables(xw),getvariables(w)));
    [Qs,cs,fs,dummy,nonquadratic] = vecquaddecomp(X,xw);
    c_wTbase = [];
    AAA = [];
    ccc = [];
    for i = 1:length(X)
        Q = Qs{i};
        c = cs{i};
        f = fs{i};
        if nonquadratic
            error('Constraints can be at most quadratic, with the linear term uncertain');
        end
        Q_ww = Q(wind,wind);
        Q_xw = Q(xind,wind);
        Q_xx = Q(xind,xind);
        c_x = c(xind);
        c_w = c(wind);

        %b = [b;f + c_w'*w];
        %A = [A;-c_x'-w'*2*Q_xw'];
        % A = [A -c_x-2*Q_xw*w];
        AAA = [AAA;sparse(-2*Q_xw)];
        ccc = [ccc;-sparse(c_x)];
        b = [b;f];
        c_wTbase = [c_wTbase;c_w'];
    end
   % b = b + c_wTbase*w;
   % A = reshape(ccc + AAA*w,size(c_x,1),[]);

   j = 1;
   confidencelevel = struct(groupedChanceConstraints{ic}).clauses{1}.confidencelevel;
   theMean    = randomVariables{j}.distribution.parameters{1};
   covariance = randomVariables{j}.distribution.parameters{2};
   gamma = icdf('normal',confidencelevel,0,1);

   switch randomVariables{j}.distribution.name
       case 'normal'
           theMean    = randomVariables{j}.distribution.parameters{1};
           covariance = randomVariables{j}.distribution.parameters{2};
           gamma = icdf('normal',confidencelevel,0,1);
           if isa(covariance,'sdpvar')
               error('Covariance cannot be an SDPVAR in normal distribution. Maybe you meant to use factorized covariance in ''normalf''');
           end
           Fderandomized = [Fderandomized, b + c_wTbase*theMean - (ccc + AAA*theMean)'*x >= gamma*norm(chol(covariance)*(AAA'*x+c_wTbase'))];
       case 'normalf'
           theMean    = randomVariables{j}.distribution.parameters{1};
           covariance = randomVariables{j}.distribution.parameters{2};
           gamma = icdf('normal',confidencelevel,0,1);
           Fderandomized = [Fderandomized, b + c_wTbase*theMean - (ccc + AAA*theMean)'*x >= gamma*norm(covariance*(AAA'*x+c_wTbase'))];
       otherwise
           error('Distribution not supported');
   end

end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function up = setup_universal_params(period_orig)
% SETUP_UNIVERSAL_PARAMS sets the universal parameters for use throughout
% "RRest".
%
%               setup_universal_params('mimic')
%
%	Inputs:
%		period_orig - the period of Vortal recording on which to do the
%		analysis, e.g. 'mimic'
%
%	Outputs:
%       up      - a struct of universal parameters
%
%   Context:    called by "RRest.m"
%
%   Further Information:
%       This version is specifically designed to facilitate reproduction of
%       the analysis performed in:
%           Charlton P.H. et al., "Waveform Analysis to Estimate
%           Respiratory Rate" [In Press]
%       Further information on this study can be obtained at:
%           http://peterhcharlton.github.io/RRest/waveform_analysis.html
%       In addition, further information on RRest, including future
%       versions, can be obtained at:
%           http://peterhcharlton.github.io/RRest/index.html
%
%   Comments, Questions, Criticisms, Feedback, Contributions:
%       See: http://peterhcharlton.github.io/RRest/contributions.html
%
%   Version:
%       v.1 - published on 23rd Feb 2016 by Peter Charlton
%
%   Licence:
%       please see the accompanying file named "LICENSE"
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% FILE PATHS TO BE SPECIFIED %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% (requires editing) %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Please note that your system may require the slashes in file paths to be
% of the opposite direction. In which case, change the following:
up.paths.slash_direction = '\';     % usually a backslash for Windows, forward slash for Linux

% Specify path of data root folder. For instance, if you specify
% "up.paths.root_folder = 'C:\Documents\Data\';", then data will be saved
% in the directory located at 'C:\Documents\Data\DATASETNAME\', where
% "DATASETNAME" is the name of the dataset being analysed, e.g. "MIMICII".
up.paths.root_folder = 'C:\Documents\Data\';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% LOAD PREVIOUS UNIVERSAL PARAMETERS %%%%%%%%
%%%%%%%%%%%%%%%% (no editing required) %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Provide licence details
provide_licence_details
% set current directory to that of this file
cd(fileparts(mfilename('fullpath')))
% Add all functions corresponding to options for each component to the path
addpath(genpath(fileparts(mfilename('fullpath'))));
% Specify name of universal params file:
up_foldername = 'universal_params';
up_filename = 'up';
up_path = [fileparts(mfilename('fullpath')), up.paths.slash_direction, up_foldername, up.paths.slash_direction, up_filename];
% Specify whether or not to re-make the universal_params (needs to be
% remade if any settings have changed)
redo_up = true;
% Don't remake the universal_params file unless needed:
if exist([up_filename, '.mat'], 'file') && ~redo_up
    fprintf('\n--- Loading Universal Parameters ');
    load(up_filename);
    return
end
% Otherwise make the universal_params file:
fprintf('\n--- Creating Universal Parameters ');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% ALGORITHMS %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% (edit if you wish) %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Specify the stages of the algorithms (best left alone):
up.al.key_components = {'extract_resp_sig', 'estimate_rr', 'fuse_rr'};      %   % To run the analysis in full this should be: {'extract_resp_sig', 'estimate_rr', 'fuse_rr'}
% Specify methods for extraction of respiratory signals (feature / filter, ecg / ppg):
up.al.options.extract_resp_sig = {'ppg_feat', 'ekg_feat'};
% Specify the components for feature-based extraction of respiratory signals:
up.al.sub_components.ppg_feat = {'EHF', 'PDt', 'FPt', 'FMe', 'RS', 'ELF'};  % Should read: {'EHF', 'PDt', 'FPt', 'FMe', 'RS', 'ELF'}
up.al.sub_components.ekg_feat = {'EHF', 'RDt', 'FPt', 'FMe', 'RS', 'ELF'};  % Should read: {'EHF', 'RDt', 'FPt', 'FMe', 'RS', 'ELF'}
% Specify the interchangeable technique(s) to be used for each component of feature-based extraction of respiratory signals:
up.al.options.PDt = {'IMS'};                                                % Possible methods: 'DCl', 'COr', 'IMS'
up.al.options.RDt = {'GC'};                                                 % Possible methods: 'GC', 'ME'
up.al.options.FMe = {'am', 'fm', 'bw'};                                     % Possible methods: 'pulW', 'am', 'fm', 'bw', 'pk', 'on', 'bwm', 'qrsW', 'qrsA', 'pca'
up.al.options.RS = {'linB'};                                                % Possible methods: 'cub', 'cubB', 'brg', 'lin', 'brgB', 'linB'
% Specify the interchangeable technique(s) to be used for filter-based respiratory signal extraction:
up.al.options.ekg_filt = {'Wfm', 'Wam', 'CCF', 'BFi'};                      % Possible methods: 'Wfm', 'Wam', 'CCF', 'BFi'
up.al.options.ppg_filt = {'Wfm', 'Wam', 'CCF', 'BFi'};                      % Possible methods: 'Wfm', 'Wam', 'CCF', 'BFi'
% Specify the interchangeable technique(s) for RR Estimation
up.al.options.estimate_rr = {'FTS', 'CtO'};
% Different methods for fusion of RR estimates:
up.al.options.fuse_rr = {'fus_mod'};                                        % Possible methods: 'fus_mod', 'fus_temp'
% Components for each method of extraction of RR fusion:
up.al.sub_components.fus_mod = {'SFu'};                                     % Possible methods: 'SFu', 'PMC', 'PRC', 'SPA'
up.al.sub_components.fus_temp = {'TFu'};                                    % Possible methods: 'TFu'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% DIRECTORIES %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% (no editing required) %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Specify the folders in which to save results files
up.paths.root_data_folder = [up.paths.root_folder, period_orig, up.paths.slash_direction];
up.paths.data_load_filename = [period_orig, '_data'];

% Specify file paths
up.paths.data_load_folder = [up.paths.root_data_folder, 'Analysis_files', up.paths.slash_direction, 'Data_for_Analysis', up.paths.slash_direction];
up.paths.data_save_folder = [up.paths.root_data_folder, 'Analysis_files', up.paths.slash_direction, 'Component_Data', up.paths.slash_direction];
up.paths.results_folder = [up.paths.root_data_folder, 'Analysis_files' , up.paths.slash_direction, 'Results', up.paths.slash_direction];
up.paths.plots_save_folder = [up.paths.results_folder, 'Figures', up.paths.slash_direction];
up.paths.tables_save_folder = [up.paths.results_folder, 'Tables', up.paths.slash_direction];
up.paths.filenames.respSigs = '_respSigs.mat';
up.paths.filenames.int_respSigs = '_int_respSigs.mat';
up.paths.filenames.pulse_peaks = '_PDt';
up.paths.filenames.qrss = '_RDt';
up.paths.filenames.fid_pts = '_FPt';
up.paths.filenames.elim_vhf = '_EHF';
up.paths.filenames.feat_meas = '_FMe';
up.paths.filenames.resampler = '_RS';
up.paths.filenames.elim_vlf2 = '_ELF';
up.paths.filenames.filt = '_flt';
up.paths.filenames.rrEsts = '_rrEsts';
up.paths.filenames.fuse = 'fuse';
up.paths.filenames.ref_rrs = '_rrRef';
up.paths.filenames.sqi = '_sqi';
up.paths.filenames.stats = '_stats';
up.paths.filenames.win_timings = '_wins';
up.paths.filenames.global_stats = 'stats_global';
up.paths.filenames.ehv_stats = 'stats_ehv';
up.paths.filenames.yhv_stats = 'stats_yhv';
up.paths.filenames.group_stats = 'stats_';
up.paths.filenames.global_data = 'data_global';
up.paths.filenames.yhv_data = 'data_yhv';
up.paths.filenames.ehv_data = 'data_ehv';
up.paths.filenames.group_data = 'data_';
up.paths.filenames.global_BA = 'BA_global';
up.paths.filenames.imp_BA = 'BA_imp';
up.paths.filenames.yhv_BA = 'BA_yhv';
up.paths.filenames.ehv_BA = 'BA_ehv';
up.paths.filenames.group_BA = 'BA_';
up.paths.filenames.study_stats = 'study_stats';
up.paths.filenames.alg_names = 'alg_names';
up.paths.filenames.imp_alg_names = 'imp_alg_names';
up.paths.filenames.results_table = 'results_table';
up.paths.filenames.BA_results_table = 'BA_results_table';
up.paths.filenames.stats_results_table = 'stats_results_table';
up.paths.filenames.subj_results_table = 'subj_results_table';
up.paths.filenames.win_results_table = 'win_results_table';
up.paths.filenames.win_results_headers = 'win_results_headers';
up.paths.filenames.win_data = 'win_data';
up.paths.filenames.win_data_imp = 'win_data_imp';
up.paths.filenames.win_data_raw_clin = 'win_data_raw_clin';
up.paths.filenames.win_data_raw_clin_imp = 'win_data_raw_clin_imp';
up.paths.filenames.vortal_results = 'vortal_results';
up.paths.filenames.synth_results = 'synth_results';
up.paths.filenames.hr_rr_scatter = 'hr_rr_scatter';
up.paths.filenames.stacked_bar = 'stacked_bar';
up.paths.filenames.joint_vortal_BA = 'BA_rest_and_rec';
up.paths.filenames.rest_and_rec_win_data = 'rest_and_rec_win_data';
up.paths.filenames.temp_prec_plot ='temp_prec_plot';
up.paths.filenames.feat_filt_plot = 'feat_filt_plot';

% Make directories
directories_to_make = {up.paths.data_load_folder, up.paths.data_save_folder, up.paths.results_folder, up.paths.plots_save_folder, up.paths.tables_save_folder};
for s = 1 : length(directories_to_make)
    rel_path = directories_to_make{s};
    if ~exist(rel_path, 'dir')
        mkdir(rel_path)
    end
end
% Move data for processing into the relevant folder (you may wish to delete the original file)
data_orig_path = [up.paths.root_folder, up.paths.data_load_filename, '.mat'];
data_new_path = [up.paths.data_load_folder, up.paths.data_load_filename, '.mat'];
if exist(data_orig_path, 'file')
    copyfile(data_orig_path, data_new_path);
end

% path of the folder in which to store results figures. This can be
% wherever you like.
up.paths.paper_figures_folder = up.paths.plots_save_folder;   % in which to save eps figures for direct import into publication.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% ANALYSIS PARAMETERS %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% (no editing required) %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Find out what subjects need processing
data_path = [up.paths.data_load_folder, up.paths.data_load_filename, '.mat'];
up.analysis.run_analysis = true;
if exist(data_path, 'file')
    load(data_path);
else
    warning('Note there was no raw data for this period...')
    path00 = [up.paths.root_folder, 'VORTAL_rest', up.paths.slash_direction];
    path0 = [path00, 'Analysis_files', up.paths.slash_direction, 'Data_for_Analysis', up.paths.slash_direction];
    path2 = 'VORTAL_rest_data';
    load([path0, path2]);
    up.analysis.run_analysis = false;
    up.paths.equipment_type = '';                                           % Possible equipment types: either '_clin' for clinical monitor, or empty, '', for raw signal acquisition
end
up.paramSet.subj_list = 1:length(data);
up.paramSet.groups = extractfield(data, 'group');

% window parameters
up.paramSet.winLeng = 32;                                                   % the duration of each window in secs
up.paramSet.winStep = 0;                                                    % the number of secs between each consecutive window
up.paramSet.buffer_period = 10;                                             % the number of secs to ignore at the start of each resp sig (since the filters might not have stabilised).

% Filter characteristics: Eliminate VHFs (above frequency content of signals)
up.paramSet.elim_vhf.ppg.Fpass = 38.5;  % in HZ
up.paramSet.elim_vhf.ppg.Fstop = 33.12;  % in HZ   (33.12 and 38.5 provide a -3 dB cutoff of 35 Hz)
up.paramSet.elim_vhf.ppg.Dpass = 0.05;
up.paramSet.elim_vhf.ppg.Dstop = 0.01;
up.paramSet.elim_vhf.ekg.Fpass = 103.7;  % in HZ
up.paramSet.elim_vhf.ekg.Fstop = 98;  % in HZ    (98 and 103.7 provide a -3 dB cutoff of 100 Hz)
up.paramSet.elim_vhf.ekg.Dpass = 0.05;
up.paramSet.elim_vhf.ekg.Dstop = 0.01;
% Filter characteristics: Eliminate Mains frequencies (50 Hz)
up.paramSet.elim_mains.fcuts = [41, 47.2, 52.8, 59];  % in HZ (41, 47.2, 52.8, 59 provide a -3 dB cutoff of 45-55 Hz)
up.paramSet.elim_mains.Dpass = 0.05;
up.paramSet.elim_mains.Dstop = 0.01;
% Filter characteristics: Eliminate HFs (above resp freqs)
up.paramSet.elim_hf.Fpass = 1.2;  % in Hz
up.paramSet.elim_hf.Fstop = 0.8899;  % in Hz     (1.2 and 0.8899 provide a -3dB cutoff of 1 Hz)
up.paramSet.elim_hf.Dpass = 0.05;
up.paramSet.elim_hf.Dstop = 0.01;
% Filter characteristics: Eliminate LFs (below cardiac freqs): For 30bpm cutoff
up.paramSet.elim_sub_cardiac.Fpass = 0.63;  % in Hz
up.paramSet.elim_sub_cardiac.Fstop = 0.43;  % in Hz     (0.63 and 0.43 provide a - 3dB cutoff of 0.5 Hz)
up.paramSet.elim_sub_cardiac.Dpass = 0.05;
up.paramSet.elim_sub_cardiac.Dstop = 0.01;
% Filter characteristics: Eliminate VLFs (below resp freqs): For 4bpm cutoff
up.paramSet.elim_vlf.Fpass = 0.157;  % in Hz
up.paramSet.elim_vlf.Fstop = 0.02;   % in Hz     (0.157 and 0.02 provide a - 3dB cutoff of 0.0665 Hz)
up.paramSet.elim_vlf.Dpass = 0.05;
up.paramSet.elim_vlf.Dstop = 0.01;
% Filter characteristics: duration of Tukey window taper in secs
up.paramSet.tukey_win_duration_taper = 2;

% frequency at which to resample feature-based respiratory signals
up.paramSet.resample_fs = 5;
% frequency at which to sample filter-based respiratory signals
up.paramSet.filt_resample_fs = 25;

% Plausible range of HRs [bpm]
up.paramSet.hr_range = [30, 220];

% Plausible range of plausible RRs [bpm]
up.paramSet.rr_range = [4 60];

% AR model order
up.paramSet.ar_model_order = 8;
up.paramSet.ar_resample_freq = 5;

% FFT downsample freq
up.paramSet.fft_resample_freq = 5;

% decide whether or not to redo stats
up.analysis.redo_stats = 1;

% QRS detection windows in secs
up.paramSet.qrs_detect.win_length = 10;
up.paramSet.qrs_detect.win_overlap = 5;

% PPG Peak Detector specifications
% CO
smooth_filt = load('parafilt'); % Load PPG smoothing filter
up.paramSet.CO_peak_det.PPG_smoother.coeffs = smooth_filt.parafiltPPG;
up.paramSet.CO_peak_det.PPG_smoother.details = smooth_filt.details;
up.paramSet.CO_peak_det.upctl = 0.9;
up.paramSet.CO_peak_det.lpctl = 0.1;

% Determine how to obtain reference RRs:
ref_fields = fieldnames(data(1).ref);
if sum(strcmp(ref_fields, 'resp_sig')) && sum(strcmp(fieldnames(data(1).ref.resp_sig), 'paw'))
    up.paramSet.ref_method = 'paw';
elseif sum(strcmp(ref_fields, 'resp_sig')) && sum(strcmp(fieldnames(data(1).ref.resp_sig), 'imp'))
    up.paramSet.ref_method = 'imp';
elseif sum(strcmp(ref_fields, 'breaths'))
    up.paramSet.ref_method = 'breaths';
elseif sum(strcmp(ref_fields, 'params')) && sum(strcmp(fieldnames(data(1).ref.params), 'rr'))
    up.paramSet.ref_method = 'rrs';
end

% Determine whether to calculate stats for impedance:
if sum(strcmp(fieldnames(data(1).ref), 'params')) && sum(strcmp(fieldnames(data(1).ref.params), 'rr')) ...
        && sum(strcmp(fieldnames(data(1).ref.params.rr), 'method')) && ~isempty(strfind(data(1).ref.params.rr.method, 'impedance'))
    up.analysis.imp_stats = true;
else
    up.analysis.imp_stats = false;
end
% Determine whether to calculate stats for synthetic data:
if ~isempty(strfind(up.paths.root_data_folder, 'synth'))
    up.analysis.calc_synth_stats = true;
else
    up.analysis.calc_synth_stats = false;
end
clear data

% Ref RR threshold
up.paramSet.resp_sig_thresh.paw = 0.42;                                     % thresh for use of paw (just a single paw, not 2*paw)
up.paramSet.resp_filter_roll_off = 1+(1/3);

% Paw SNR
up.paramSet.paw_snr_thresh = -8;

% Reference RR method
up.paramSet.ref_rr_method = 'paw_thresh';                                   % For Vortal: sum_both_thresh, imp_thresh, paw_thresh; For MIMIC: comb_imp; For CAPNOBASE: co2_zex

% Impedance RR methods
up.paramSet.imp_rr_method = 'num';   % num or thresh

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% SAVE UNIVERSAL PARAMETERS %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% (no editing required) %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

save(up_path, 'up');
end

function provide_licence_details

licence_details = ['\n\n RRest  Copyright (C) 2016 King''s College London', ...
    '\n This program comes with ABSOLUTELY NO WARRANTY', ...
    '\n and is available under the GNU public license.', ...
    '\n For details see the accompanying LICENSE.txt file.\n\n'];

fprintf(licence_details)

end