text:
function [fx,tt]=fxrapt(s,fs,mode)
%FXRAPT RAPT pitch tracker [FX,VUV]=(S,FS)
%
% Input:   s(ns)      Speech signal
%          fs         Sample frequency (Hz)
%          mode       'g' will plot a graph [default if no output arguments]
%                     'u' will include unvoiced fames (with fx=NaN)
%
% Outputs: fx(nframe)     Larynx frequency for each fram,e (or NaN for silent/unvoiced)
%          tt(nframe,3)  Start and end samples of each frame. tt(*,3)=1 at the start of each talk spurt
%
% Plots a graph if no outputs are specified showing lag candidates and selected path
%

% Bugs/Suggestions:
%   (1) Include backward DP pass and output the true cost for each candidate.
%   (2) Add an extra state to distinguish between voiceless and silent
%   (3) N-best DP to allow longer term penalties (e.g. for frequent pitch doubling/halving)

% The algorithm is taken from [1] with the following differences:
%
%      (a)  the factor AFACT which in the Talkin algorithm corresponds roughly
%           to the absolute level of harmonic noise in the correlation window. This value
%           is here calculated as the maximum of three figures:
%                   (i) an absolute floor set by PP.rapt_absnoise
%                  (ii) a multiple of the peak signal set by PP.rapt_signoise
%                 (iii) a multiple of the noise floor set by PP.rapt_relnoise
%      (b) The LPC used in calculating the Itakura distance uses a Hamming window rather than
%          a Hanning window.
%
% A C implementation of this algorithm by Derek Lin and David Talkin is included as  "get_f0.c"
% in the esps.zip package available from http://www.speech.kth.se/esps/esps.zip under the BSD
% license.
%
% Refs:
%      [1]   D. Talkin, "A Robust Algorithm for Pitch Tracking (RAPT)"
%            in "Speech Coding & Synthesis", W B Kleijn, K K Paliwal eds,
%            Elsevier ISBN 0444821694, 1995

%      Copyright (C) Mike Brookes 2006
%      Version: $Id: fxrapt.m 3601 2013-10-11 15:27:30Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s=s(:); % force s to be a column
if nargin<3
    mode=' ';
end
doback=0;   % don't do backwards DP for now

% read in parameters

PP=voicebox;
f0min=PP.rapt_f0min;            % Min F0 (Hz)                               [50]
f0max=PP.rapt_f0max;            % Max F0 (Hz)                               [500]
tframe=PP.rapt_tframe;          % frame size (s)                            [0.01]
tlpw=PP.rapt_tlpw;              % low pass filter window size (s)           [0.005]
tcorw=PP.rapt_tcorw;            % correlation window size (s)               [0.0075]
candtr=PP.rapt_candtr;          % minimum peak in NCCF                      [0.3]
lagwt=PP.rapt_lagwt;            % linear lag taper factor                   [0.3]
freqwt=PP.rapt_freqwt;          % cost factor for F0 change                 [0.02]
vtranc=PP.rapt_vtranc;          % fixed voice-state transition cost         [0.005]
vtrac=PP.rapt_vtrac;            % delta amplitude modulated transition cost [0.5]
vtrsc=PP.rapt_vtrsc;            % delta spectrum modulated transition cost  [0.5]
vobias=PP.rapt_vobias;          % bias to encourage voiced hypotheses       [0.0]
doublec=PP.rapt_doublec;        % cost of exact doubling or halving         [0.35]
absnoise=PP.rapt_absnoise;      % absolute rms noise level                  [0]
relnoise=PP.rapt_relnoise;      % rms noise level relative to noise floor   [2.0]
signoise=PP.rapt_signoise;      % ratio of peak signal rms to noise floor   [0.001]
ncands=PP.rapt_ncands;          % max hypotheses at each frame              [20]
trms=PP.rapt_trms;              % window length for rms measurement         [0.03]
dtrms=PP.rapt_dtrms;            % window spacing for rms measurement        [0.02]
preemph=PP.rapt_preemph;        % s-plane position of preemphasis zero      [-7000]
nfullag=PP.rapt_nfullag;        % number of full lags to try (must be odd)  [7]

% derived parameters (mostly dependent on sample rate fs)

krms=round(trms*fs);            % window length for rms measurement
kdrms=round(dtrms*fs);          % window spacing for rms measurement
rmswin=hanning(krms).^2;
kdsmp=round(0.25*fs/f0max);
hlpw=round(tlpw*fs/2);          % force window to be an odd length
blp=sinc((-hlpw:hlpw)/kdsmp).*hamming(2*hlpw+1).';
fsd=fs/kdsmp;
kframed=round(fsd*tframe);      % downsampled frame length
kframe=kframed*kdsmp;           % frame increment at full rate
rmsix=(1:krms)+floor((kdrms-kframe)/2); % rms index according to Talkin; better=(1:krms)+floor((kdrms-krms+1)/2)
minlag=ceil(fsd/f0max);
maxlag=round(fsd/f0min);        % use round() only because that is what Talkin does
kcorwd=round(fsd*tcorw);        % downsampled correlation window
kcorw=kcorwd*kdsmp;             % full rate correlation window
spoff=max(hlpw-floor(kdsmp/2),1+kdrms-rmsix(1)-kframe);  % offset for first speech frame at full rate
sfoff=spoff-hlpw+floor(kdsmp/2); % offset for downsampling filter
sfi=1:kcorwd;                   % initial decimated correlation window index array
sfhi=1:kcorw;                   % initial correlation window index array
sfj=1:kcorwd+maxlag;
sfmi=repmat((minlag:maxlag)',1,kcorwd)+repmat(sfi,maxlag-minlag+1,1);
lagoff=(minlag-1)*kdsmp;        % lag offset when converting to high sample rate
beta=lagwt*f0min/fs;            % bias towards low lags
log2=log(2);
lpcord=2+round(fs/1000);        % lpc order for itakura distance
hnfullag=floor(nfullag/2);
jumprat=exp((doublec+log2)/2);  % lag ratio at which octave jump cost is lowest
ssq=s.^2;
csssq=cumsum(ssq);
sqrt(min(csssq(kcorw+1:end)-csssq(1:end-kcorw))/kcorw);
afact=max([absnoise^2,max(ssq)*signoise^2,min(csssq(kcorw+1:end)-csssq(1:end-kcorw))*(relnoise/kcorw)^2])^2*kcorw^2;

% downsample signal to approx 2 kHz to speed up autocorrelation calculation
% kdsmp is the downsample factor

sf=filter(blp/sum(blp),1,s(sfoff+1:end));
sp=filter([1 exp(preemph/fs)],1,s); % preemphasised speech for LPC calculation
sf(1:length(blp)-1)=[];         % remove startup transient
sf=sf(1:kdsmp:end);             % downsample to =~2kHz
nsf=length(sf);                 % length of downsampled speech
ns=length(s);                   % length of full rate speech

% Calculate the frame limit to ensure we don't run off the end of the speech or decimated speech:
%   (a) For decimated autocorrelation when calculating sff():  (nframe-1)*kframed+kcorwd+maxlag <= nsf
%   (b) For full rate autocorrelation when calculating sfh():  max(fho)+kcorw+maxlag*kdsamp+hnfllag <= ns
%   (c) For rms ratio window when calculating rr            :  max(fho)+rmsix(end) <= ns
% where max(fho) = (nframe-1)*kframe + spoff

nframe=floor(1+min((nsf-kcorwd-maxlag)/kframed,(ns-spoff-max(kcorw-maxlag*kdsmp-hnfullag,rmsix(end)))/kframe));

% now search for autocorrelation peaks in the downsampled signal

cost=zeros(nframe,ncands);      % cumulative cost
prev=zeros(nframe,ncands);      % traceback pointer
mcands=zeros(nframe,1);         % number of actual candidates excluding voiceless
lagval=repmat(NaN,nframe,ncands-1);    % lag of each voiced candidate
tv=zeros(nframe,3);             % diagnostics: 1=voiceless cost, 2=min voiced cost, 3:cumulative voiceless-min voiced
if doback
    costms=cell(nframe,1);
end

% Main processing loop for each 10 ms frame

for iframe=1:nframe       % loop for each frame (~10 ms)

    % Find peaks in the normalized autocorrelation of subsampled (2Khz) speech
    % only keep peaks that are > 30% of highest peak

    fho=(iframe-1)*kframe+spoff;
    sff=sf((iframe-1)*kframed+sfj);
    sffdc=mean(sff(sfi));       % mean of initial correlation window length
    sff=sff-sffdc;              % subtract off the mean
    nccfd=normxcor(sff(1:kcorwd),sff(minlag+1:end));
    [ipkd,vpkd]=v_findpeaks(nccfd,'q');

    % Debugging: execute the line below to plot the autocorrelation peaks.
    % v_findpeaks(nccfd,'q'); xlabel(sprintf('Lag = (x+%d)*%g ms',minlag-1,1000*kdsmp/fs)); ylabel('Normalized Cross Correlation'); title (sprintf('Frame %d/%d',iframe,nframe));

    vipkd=[vpkd ipkd];
    vipkd(vpkd<max(vpkd)*candtr,:)=[];          % eliminate peaks that are small
    if size(vipkd,1)
        if size(vipkd,1)>ncands-1
            vipkd=sortrows(vipkd);
            vipkd(1:size(vipkd,1)-ncands+1,:)=[];   % eliminate lowest to leave only ncands-1
        end
        lagcan=round(vipkd(:,2)*kdsmp+lagoff);        % convert the lag candidate values to the full sample rate
        nlcan=length(lagcan);
    else
        nlcan=0;
    end

    % If there are any candidate lag values (nlcan>0) then refine their accuracy at the full sample rate

    if nlcan
        laglist=reshape(repmat(lagcan(:)',nfullag,1)+repmat((-hnfullag:hnfullag)',1,nlcan),nfullag*nlcan,1);
        sfh=s(fho+(1:kcorw+max(lagcan)+hnfullag));
        sfhdc=mean(sfh(sfhi));
        sfh=sfh-sfhdc;
        e0=sum(sfh(sfhi).^2);                     % energy of initial correlation window (only needed to store in tv(:,6)
        lagl2=repmat(lagcan(:)',nfullag+kcorw-1,1)+repmat((1-hnfullag:hnfullag+kcorw)',1,nlcan);
        nccf=normxcor(sfh(1:kcorw),sfh(lagl2),afact);

        [maxcc,maxcci]=max(nccf,[],1);
        vipk=[maxcc(:) lagcan(:)+maxcci(:)-hnfullag-1];
        vipk=vipk(:,[1 2 2]);
        maxccj=maxcci(:)'+nfullag*(0:nlcan-1);    % vector index into nccf array
        msk=mod(maxcci,nfullag-1)~=1 & 2*nccf(maxccj)-nccf(mod(maxccj-2,nfullag*nlcan)+1)-nccf(mod(maxccj,nfullag*nlcan)+1)>0;  % don't do quadratic interpolation for the end ones
        if any(msk)
            maxccj=maxccj(msk);
            vipk(msk,3)=vipk(msk,3)+(nccf(maxccj+1)-nccf(maxccj-1))'./(2*(2*nccf(maxccj)-nccf(maxccj-1)-nccf(maxccj+1)))';
        end
        vipk(maxcc<max(maxcc)*candtr,:)=[];          % eliminate peaks that are small
        if size(vipk,1)>ncands-1
            vipk=sortrows(vipk);
            vipk(1:size(vipk,1)-ncands+1,:)=[];   % eliminate lowest to leave only ncands-1
        end

        % vipk(:,1) has NCCF value, vipk(:,2) has integer peak position, vipk(:,3) has refined peak position

        mc=size(vipk,1);
    else
        mc=0;
    end

    % We now have mc lag candidates at the full sample rate

    mc1=mc+1;               % total number of candidates including "unvoiced" possibility
    mcands(iframe)=mc;      % save number of lag candidates (needed for pitch consistency cost calculation)
    if mc
        lagval(iframe,1:mc)=vipk(:,3)';
        cost(iframe,1)=vobias+max(vipk(:,1));   % voiceless cost
        cost(iframe,2:mc1)=1-vipk(:,1)'.*(1-beta*vipk(:,3)');   % local voiced costs
        tv(iframe,2)=min(cost(iframe,2:mc1));
    else
        cost(iframe,1)=vobias;          % if no lag candidates (mc=0), then the voiceless case is the only possibility
    end
    tv(iframe,1)=cost(iframe,1);
    if iframe>1                         % if it is not the first frame, then calculate pitch consistency and v/uv transition costs
        mcp=mcands(iframe-1);
        costm=zeros(mcp+1,mc1);         % cost matrix: rows and cols correspond to candidates in previous and current frames (incl voiceless)

        % if both frames have at least one lag candidate, then calculate a pitch consistency cost

        if mc*mcp
            lrat=abs(log(repmat(lagval(iframe,1:mc),mcp,1)./repmat(lagval(iframe-1,1:mcp)',1,mc)));
            costm(2:end,2:end)=freqwt*min(lrat,doublec+abs(lrat-log2));  % allow pitch doubling/halving
        end

        % if either frame has a lag candidate, then calculate the cost of voiced/voiceless transition and vice versa

        if mc+mcp
            rr=sqrt((rmswin'*s(fho+rmsix).^2)/(rmswin'*s(fho+rmsix-kdrms).^2)); % amplitude "gradient"
            ss=0.2/(distitar(lpcauto(sp(fho+rmsix),lpcord),lpcauto(sp(fho+rmsix-kdrms),lpcord),'e')-0.8);   % Spectral stationarity: note: Talkin uses Hanning instead of Hamming windows for LPC
            costm(1,2:end)= vtranc+vtrsc*ss+vtrac/rr;   % voiceless -> voiced cost
            costm(2:end,1)= vtranc+vtrsc*ss+vtrac*rr;
            tv(iframe,4:5)=[costm(1,mc1) costm(mcp+1,1)];
        end
        costm=costm+repmat(cost(iframe-1,1:mcp+1)',1,mc1);  % add in cumulative costs
        [costi,previ]=min(costm,[],1);
        cost(iframe,1:mc1)=cost(iframe,1:mc1)+costi;
        prev(iframe,1:mc1)=previ;
    else                            % first ever frame
        costm=zeros(1,mc1); % create a cost matrix in case doing a backward recursion
    end
    if mc
        tv(iframe,3)=cost(iframe,1)-min(cost(iframe,2:mc1));
        tv(iframe,6)=5*log10(e0*e0/afact);
    end
    if doback
        costms{iframe}=costm; % need to add repmatted cost into this
    end
end

% now do traceback

best=zeros(nframe,1);
[cbest,best(nframe)]=min(cost(nframe,1:mcands(nframe)+1));
for i=nframe:-1:2
    best(i-1)=prev(i,best(i));
end
vix=find(best>1);
fx=repmat(NaN,nframe,1);                        % unvoiced frames will be NaN
fx(vix)=fs*lagval(vix+nframe*(best(vix)-2)).^(-1); % leave as NaN if unvoiced
tt=zeros(nframe,3);
tt(:,1)=(1:nframe)'*kframe+spoff;       % find frame times
tt(:,2)=tt(:,1)+kframe-1;
jratm=(jumprat+1/jumprat)/2;
tt(2:end,3)=abs(fx(2:end)./fx(1:end-1)-jratm)>jumprat-jratm;    % new spurt if frequency ratio is outside (1/jumprat,jumprat)
tt(1,3)=1;           % first frame always starts a spurt
tt(1+find(isnan(fx(1:end-1))),3)=1; % NaN always forces a new spurt

% plot results if there are no output arguments of if the 'g' mode option is specified

if ~nargout || any(mode=='g')
    tf=spoff+(0:nframe-1)'*kframe;      % one sample before start of each frame
    blag=repmat(NaN,nframe,1);                        % unvoiced frames will be NaN
    blag(vix)=lagval(vix+nframe*(best(vix)-2)); % leave as NaN if unvoiced
    ts=(1:ns)/fs;                       % time scale for speech samples
    tsa=[1:tf(1) tf(end)+kframe+1:ns];  % indexes for unprocessed speech [-1 term is an error methinks]
    sup=repmat(NaN,ns,1);               % unprocessed speech - plot in black
    sup(tsa)=s(tsa);
    sv=reshape(s(tf(1)+1:tf(end)+kframe),kframe,nframe);               % processed speech
    su=sv;
    su(:,best>1)=NaN;                   % delete all voiced samples
    sv(:,best==1)=NaN;                  % delete all unvoiced samples
    tsuv=(tf(1)+1:tf(end)+kframe)/fs;
    su=su(:);
    sv=sv(:);
    ax=zeros(2,1);
    ax(1)=subplot(211);
    plot(ts,sup,'-k',tsuv,su,'r-',tsuv,sv,'b-');
    title('Speech');
    ax(2)=subplot(212);
    plot((tf+(kframe+1)/2)/fs,lagval*1000/fs,'xr',(tf+(kframe+1)/2)/fs,blag*1000/fs,'-b')
    xlabel('Time (s)');
    ylabel('Period (ms)');
    title('Lag Candidates');
    linkaxes(ax,'x');
end
if ~any(mode=='u')
tt(isnan(fx),:)=[];    % remove NaN spurts
fx(isnan(fx),:)=[];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function v=normxcor(x,y,d)
% Calculate the normalized cross correlation of column vectors x and y
% we can calculate this in two ways but fft is much faster even for nx small
% We must have nx<=ny and the output length is ny-nx+1
% note that this routine does not do mean subtraction even though this is normally a good idea
% if y is a matrix, we correlate with each column
% d is a constant added onto the normalization factor
% v(j)=x'*yj/sqrt(d + x'*x * yj'*yj) where yj=y(j:j+nx-1) for j=1:ny-nx+1

if nargin<3
    d=0;
end
nx=length(x);
[ny,my]=size(y);
nv=1+ny-nx;
if nx>ny
    error('second argument is shorter than the first');
end

nf=pow2(nextpow2(ny));
w=irfft(repmat(conj(rfft(x,nf,1)),1,my).*rfft(y,nf,1));
s=zeros(ny+1,my);
s(2:end,:)=cumsum(y.^2,1);
v=w(1:nv,:)./sqrt(d+(x'*x).*(s(nx+1:end,:)-s(1:end-nx,:)));

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
clear all
close all

tora_tp

tora_lmi

tora_sim

author: Sable

repo: mcbench-benchmarks


text:
function medit_io_test01 ( )

%*****************************************************************************80
%
%% MEDIT_IO_TEST01 creates a MESH dataset and writes it to a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 December 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MEDIT_IO_TEST01:\n' );
  fprintf ( 1, '  Create a hexahedral mesh and write it to a file.\n' );

  [ dim, vertices, edges, triangles, quadrilaterals, tetrahedrons, ...
    hexahedrons, vertex_coordinate, vertex_label, edge_vertex, edge_label, ...
    triangle_vertex, triangle_label, quadrilateral_vertex, ...
    quadrilateral_label, tetrahedron_vertex, tetrahedron_label, ...
    hexahedron_vertex, hexahedron_label ] = hexahexa_2x2x2 ( );

  filename = 'hexahexa_2x2x2.mesh';

  mesh_write ( filename, dim, vertices, edges, triangles, quadrilaterals, ...
    tetrahedrons, hexahedrons, vertex_coordinate, vertex_label, edge_vertex, ...
    edge_label, triangle_vertex, triangle_label, quadrilateral_vertex, ...
    quadrilateral_label, tetrahedron_vertex, tetrahedron_label, ...
    hexahedron_vertex, hexahedron_label );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Created the file "%s".\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef ResponseLossL2 < dagnn.Loss
%ResponseLossL2  layer
%  `ResponseLossL2.forward({r, x*})` computes L2 loss.
%
%  Here the Loss between two matrices is defined as:
%
%     Loss = |r - r_idea|.*|r - r_idea|
%
%  Input
%       - r    w x h x 1 x N.
%       - x0   N x 2
%   QiangWang, 2016
% -------------------------------------------------------------------------------------------------------------------------
    properties
        win_size = [3,3];
        sigma = 1;
    end
    properties (Transient)
        y = [];
    end
    methods
        function outputs = forward(obj, inputs, params)
            r1 = inputs{1}; % middle result of Z1
            r2 = inputs{2}; % middle result of Z2
            r3 = inputs{3}; % final label of X
            r4 = inputs{4}; % initial label of X

            % original loss
            res_minus = bsxfun(@minus,r3,r4);
            loss = res_minus.*res_minus;
            perSample_loss = sum(sum(sum(loss,1),2),3);
            [~, index] = sort(perSample_loss, 'ascend');
            sample_drop = zeros(1,1,1,size(r1,4));
            % drop unreliable samples
            max_num = round(0.9*size(r1,4));
            sample_drop(index(1:max_num)) = 1;

            % the motion between X and Z1
            delta1 = bsxfun(@minus,r1,r4);
            % the motion between Z1 and Z2
            delta2 = bsxfun(@minus,r1,r2);
            delta = delta1.*delta1 + delta2.*delta2;
            delta_sum = sum(sum(sum(delta, 1),2),3);
            delta_drop = bsxfun(@times,delta_sum,sample_drop);
            delta_norm = delta_drop./sum(delta_drop(:)) * size(r3,4);

            % spatial weight
            spatial = exp(r4);
            weight = bsxfun(@times, delta_norm, spatial);
            % add weight
            loss = loss.*weight;

            outputs{1} = sum(loss(:))/size(r3,4);

            n = obj.numAveraged ;
            m = n + 1 ;
            obj.average = (n * obj.average + gather(outputs{1})) / m ;
            obj.numAveraged = m ;

        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
            r1 = inputs{1};
            r2 = inputs{2};
            r3 = inputs{3};
            r4 = inputs{4};

            % original loss
            res_minus = bsxfun(@minus,r3,r4);
            loss = res_minus.*res_minus;
            perSample_loss = sum(sum(sum(loss,1),2),3);
            [~, index] = sort(perSample_loss, 'ascend');
            sample_drop = zeros(1,1,1,size(r1,4));
            % drop unreliable samples
            max_num = round(0.9*size(r1,4));
            sample_drop(index(1:max_num)) = 1;

            % the motion between X and Z1
            delta1 = bsxfun(@minus,r1,r4);
            % the motion between Z1 and Z2
            delta2 = bsxfun(@minus,r1,r2);
            delta = delta1.*delta1 + delta2.*delta2;
            delta_sum = sum(sum(sum(delta, 1),2),3);
            delta_drop = bsxfun(@times,delta_sum,sample_drop);
            delta_norm = delta_drop./sum(delta_drop(:)) * size(r3,4);

            % spatial weight
            spatial = exp(r4);
            weight = bsxfun(@times, delta_norm, spatial);

            derInputs{1} = {};
            derInputs{2} = {};
            derInputs{3} = (derOutputs{1}*2/size(r3,4)).*bsxfun(@minus,r3,r4).*weight ;
            derInputs{4} = {};
            derParams = {} ;
        end

        function obj = reset(obj)
            obj.y = [] ;
            obj.average = 0 ;
            obj.numAveraged = 0 ;
        end

        function obj = ResponseLossL2(varargin)
            obj.load(varargin);
            obj.win_size = obj.win_size;
            obj.sigma = obj.sigma ;
            obj.y = [];
        end

    end

end

repo: UDT


text:
function SimVaryResults = SimVary(obj, runs, OptTable, Opts)
% Performs Sensitivity analysis on Model object
% Each Parameter (obj.xnames) is varied in 10 steps between OptTable.lb and OptTable.ub.
% Noise is added N=runs times and simulated MR signals is fitted using obj.Sim_Single_Voxel_Curve
%
% USAGE: SimVaryResults = SimVary(obj, runs, OptTable, Opts)
%
%        obj: qMR Model
%        runs: Number of runs
%        OptTable: structure with vectors defining options for each fitting parameters obj.xnames
%               OptTable.fx: true/false. Vary this parameter?
%               OptTable.st: Nominal Values
%               OptTable.lb: Vary from lb...
%               OptTable.ub: ..to ub value
%        Opts.SNR
%

Nsteps = 10;

if ~exist('OptTable','var') || isempty(OptTable), OptTable = obj; end % use fitting boundaries
if ~exist('Opts','var') || isempty(Opts), Opts.SNR = 50; end
if isempty(getenv('ISCITEST')) || str2double(getenv('ISCITEST'))==0
    waitbarcreate = true;
else
    waitbarcreate = false;
    Nsteps=3;
end

if ~isempty(getenv('ISDOC')), Nsteps=10; end

fx = [OptTable.fx];
st = [OptTable.st];
lb = [OptTable.lb];
ub = [OptTable.ub];

for pp=1:length(fx)
    if ~fx(pp)
        Sens.x = linspace(lb(pp),ub(pp),Nsteps);
        % Create waitbar
        if waitbarcreate
            h = waitbar(0, sprintf('Data 0/%0.0f',length(Sens.x)), 'Name', sprintf('Simulating %s sensitivity data', obj.xnames{pp}),...
                'CreateCancelBtn', 'if ~strcmp(get(gcbf,''Name''),''canceling...''), setappdata(gcbf,''canceling'',1); set(gcbf,''Name'',''canceling...''); else delete(gcbf); end');
            setappdata(h,'canceling',0);
        end
        setappdata(0,'Cancel',0);

        for ii=1:length(Sens.x)
            x = st; x(pp)=Sens.x(ii);
            for N=1:runs
                if waitbarcreate && (~ishandle(h) || getappdata(h,'canceling')); setappdata(0,'Cancel',1); break; end
                Fittmp = obj.Sim_Single_Voxel_Curve(x, Opts,0);
                if ~isfield(Sens,'fit')
                    Sens.fit = Fittmp;
                else
                    for ff=fieldnames(Fittmp)'
                        Sens.fit.(ff{1})(ii,N) = Fittmp.(ff{1})(1);
                    end
                end
                %if ~ishandle(h), return; end
            end
            if waitbarcreate && (~ishandle(h) || getappdata(h,'canceling')); break; end
            if waitbarcreate && ishandle(h)
                waitbar(ii/length(Sens.x),h,sprintf('Data %0.0f/%0.0f',ii,length(Sens.x)));
            end
        end

        if waitbarcreate
            delete(h)       % DELETE the waitbar; don't try to CLOSE it.
        end

        for ff=fieldnames(Sens.fit)'
            Sens.(ff{1}).mean = mean(Sens.fit.(ff{1}),2);
            Sens.(ff{1}).std = std(Sens.fit.(ff{1}),0,2);
            Sens.(ff{1}).GroundTruth = st(strcmp(obj.xnames,ff{1}));
        end
        SimVaryResults.(obj.xnames{pp})=Sens;
    end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function x = conj_gradient_method(A, b)

%Parameters
n = 20;

%Initialization
clc
format short e
tol = 1e-2;

r = A*x-b;
p = -r;
k = 0;
abs_r = sqrt(r'*r);

while abs_r > tol

    alpha = (r'*r)/(p'*A*p);

    x = x+alpha*p;
    r_old = r;
    r = r + alpha*A*p;
    beta = (r'*r)/(r_old'*r_old);
    p = -r+beta*p;
    k=k+1;

    abs_r = sqrt(r'*r);
%    disp([k abs_r alpha]);
end

author: thomas-koehler

repo: SupER


text:
%RPY2TR Roll-pitch-yaw angles to SE(3) homogeneous transform
%
% T = RPY2TR(ROLL, PITCH, YAW, OPTIONS) is an SE(3) homogeneous
% transformation matrix (4x4) with zero translation and rotation equivalent
% to the specified roll, pitch, yaw angles angles. These correspond to
% rotations about the Z, Y, X axes respectively. If ROLL, PITCH, YAW are
% column vectors (Nx1) then they are assumed to represent a trajectory and
% R is a three-dimensional matrix (4x4xN), where the last index corresponds
% to rows of ROLL, PITCH, YAW.
%
% T = RPY2TR(RPY, OPTIONS) as above but the roll, pitch, yaw angles are
% taken from the vector (1x3) RPY=[ROLL,PITCH,YAW]. If RPY is a matrix
% (Nx3) then R is a three-dimensional matrix (4x4xN), where the last index
% corresponds to rows of RPY which are assumed to be
% ROLL,PITCH,YAW].
%
% Options::
%  'deg'      Compute angles in degrees (radians default)
%  'xyz'      Rotations about X, Y, Z axes (for a robot gripper)
%  'zyx'      Rotations about Z, Y, X axes (for a mobile robot, default)
%  'yxz'      Rotations about Y, X, Z axes (for a camera)
%  'arm'      Rotations about X, Y, Z axes (for a robot arm)
%  'vehicle'  Rotations about Z, Y, X axes (for a mobile robot)
%  'camera'   Rotations about Y, X, Z axes (for a camera)
%
% Note::
% - Toolbox rel 8-9 has the reverse angle sequence as default.
% - ZYX order is appropriate for vehicles with direction of travel in the X
%   direction.  XYZ order is appropriate if direction of travel is in the Z
%   direction.
% - 'arm', 'vehicle', 'camera' are synonyms for 'xyz', 'zyx' and 'yxz'
%   respectively.
%
% See also TR2RPY, RPY2R, EUL2TR.

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% 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.
%
% https://github.com/petercorke/spatial-math

function T = rpy2tr(roll, varargin)

    R = rpy2r(roll, varargin{:});
    T = r2t(R);

author: petercorke

repo: spatialmath-matlab


text:
function x = stretch_rs ( r, s, i )

%*****************************************************************************80
%
%% STRETCH_RS returns a data component given (R,S).
%
%  Discussion:
%
%    This routine shifts by (1,2) and stretches by (3,4).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, S, the coordinates of a point that lies on the
%    boundary of the square.
%
%    Input, integer I, the component of the data which is to be returned.
%
%    Output, real X, the I-th component of the data vector X, evaluated
%    at the point (R,S), which is on a corner, or edge, of the unit square.
%
  if ( i == 1 )
    x = 3.0 * r + 1.0;
  elseif ( i == 2 )
    x = 4.0 * s + 2.0;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'STRETCH_RS - Fatal error!\n' );
    fprintf ( 1, '  Illegal component index I = %d\n', i );
    error ( 'STRETCH_RS - Fatal error!' );
  end

  return
end


text:
filename='CantileverBeam_Triangle_Linear_Fine';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = 1;
%cost = {'compliance'};
%weights = [1];
constraint = {'volumeConstraint'};
%optimizer = 'DualNestedInPrimal';
optimizer = 'AlternatingPrimalDual';

optimizerUnconstrained = 'PROJECTED GRADIENT';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.4;
optimality_final =1e-3;
constr_final = 1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

plotting = false;
printing = false;
monitoring = false;

maxiter = 200;

monitoring_interval = 1;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [train_users, dev_users] = get_balanced_fold(BP4D_dir, users, au, prop_test)

    % Extracting the labels
    [labels, valid_ids, vid_ids, filenames] = extract_BP4D_labels(BP4D_dir, users, au);

    % the grouping should be done per person

    for f=1:numel(filenames)
        filenames{f} = filenames{f}(1:4);
    end

    counts = zeros(numel(users),1);
    for k=1:numel(users)
        counts(k) = sum(cat(1, labels{strcmp(filenames, users{k})}));
    end

    [sorted, inds] = sort(counts);

    dev_users = users(inds(1:round(1/prop_test):end));
    train_users = setdiff(users, dev_users);
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function r = ldivide(a,b)
%LDIVIDE      Taylor elementwise left division a .\ b
%

% written  05/21/09     S.M. Rump
%

  r = b ./ a;

repo: OpenMAS


text:
function fd1d_burgers_lax ( )

%*****************************************************************************80
%
%% FD1D_BURGERS_LAX solves the nonviscous Burgers equation using Lax-Wendroff.
%
%  Discussion:
%
%    The non-viscous time-dependent Burgers equation is:
%
%      du/dt + u du/dx = 0
%
%    which can be written in conservative form as
%
%      du/dt + 1/2 d/dx ( u^2 ) = 0
%
%    or
%
%      du/dt + dF/dx = 0
%
%    For the Burgers equation, we define
%
%      F(x,t) = 1/2 u^2,
%      A(x,t) = dF/dx = u
%
%    and then the Lax-Wendroff method approximates the solution
%    using the iteration:
%
%      u(x,t+dt) = u(t) - dt dF/dx + 1/2 dt^2 d/dx A dF/dx
%
%    which can be written:
%
%      u(x,t+dt) = u(x,t) - dt ( F(x+dx,t) - F(x-dx,t) ) / ( 2 * dx )
%        + 1/2 dt^2/dx^2 ( A(x+dx/2,t) * ( F(x+dx,t) - F(x,t) )
%                        - A(x-dx/2,t) * ( F(x,t) - F(x-dx,t) )
%
%    where we approximate:
%
%      A(x+dx/2,t) = 1/2 ( u(x+dx,t) + u(x,t) )
%      A(x-dx/2,t) = 1/2 ( u(x,t) + u(x-dx,t) )
%
%    There is a stability condition that applies here, which requires that
%
%      dt * max ( abs ( u ) ) / dx <= 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_BURGERS_LAX:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Solve the non-viscous time-dependent Burgers equation,\n' );
  fprintf ( 1, '  using the Lax-Wendroff method.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Equation to be solved:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    du/dt + u * du/dx = 0\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  for x in [ a, b ], for t in [t_init, t_last]\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  with initial conditions:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    u(x,o) = u_init\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  and boundary conditions:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    u(a,t) = u_a(t), u(b,t) = u_b(t)\n' );
%
%  Set and report the problem parameters.
%
  n = 41;
  a = -1.0;
  b = +1.0;
  dx = ( b - a ) / ( n - 1 );
  step_num = 80;
  t_init = 0.0;
  t_last = 1.0;
  dt = ( t_last - t_init ) / step_num;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  %f <= X <= %f\n', a, b );
  fprintf ( 1, '  Number of nodes = %d\n', n );
  fprintf ( 1, '  DX = %f\n', dx );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %f <= t <= %f\n', t_init, t_last );
  fprintf ( 1, '  Number of time steps = %d\n', step_num );
  fprintf ( 1, '  DT = %f\n', dt );

  x = r8vec_even ( n, a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  X:\n' );
  fprintf ( 1, '\n' );
  for ilo = 1 : 5 : n
    ihi = min ( ilo + 4, n );
    for i = ilo : ihi
      fprintf ( 1, '  %16f', x(i,1) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Set the initial condition,
%  and apply boundary conditions to first and last entries.
%
  step = 0;
  t = t_init;
  un(1:n,1) = u_init ( n, x, t );
  un(1,1) = u_a ( x(1,1), t );
  un(n,1) = u_b ( x(n,1), t );

  stability = ( dt / dx ) * max ( abs ( un(1:n,1) ) );
  report ( step, step_num, n, x, t, un, stability );

  if ( 1 )
    plot ( x, un );
    title ( sprintf ( 'Step %d, Time %f', step, t ) );
    pause
  end
%
%  March in time.
%
  c1 = - ( 0.5 * dt   / dx );
  c2 = - ( 0.5 * dt^2 / dx^2 );

  for step = 1 : step_num

    t = ( ( step_num - step ) * t_init   ...
        + (            step ) * t_last ) ...
        / ( step_num        );

    uo(1:n,1) = un(1:n,1);

    un(2:n-1,1) = uo(2:n-1,1) ...
      -       ( dt    / dx  ) * ( uo(3:n,1).^2 - uo(1:n-2,1).^2 ) ...
      + 0.5 * ( dt^2 / dx^2 ) * (  0.5 * ( uo(3:n,1)   + uo(2:n-1,1) ) ...
                                   .* ( uo(3:n,1).^2   - uo(2:n-1,1).^2 ) ...
                                 - 0.5 * ( uo(2:n-1,1) + uo(1:n-2,1) ) ...
                                   .* ( uo(2:n-1,1).^2 - uo(1:n-2,1).^2 ) );

    un(1,1) = u_a ( x(1,1), t );
    un(n,1) = u_b ( x(n,1), t );

    stability = ( dt / dx ) * max ( abs ( un(1:n,1) ) );
    report ( step, step_num, n, x, t, un, stability );

    if ( 1 )
      plot ( x, un );
      title ( sprintf ( 'Step %d, Time %f', step, t ) );
      pause
    end

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_BURGERS_LAX:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

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

  return
end
function a = r8vec_even ( n, alo, ahi )

%*****************************************************************************80
%
%% R8VEC_EVEN returns N real values, evenly spaced between ALO and AHI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of values.
%
%    Input, real ALO, AHI, the low and high values.
%
%    Output, real A(N), N evenly spaced values.
%    Normally, A(1) = ALO and A(N) = AHI.
%    However, if N = 1, then A(1) = 0.5*(ALO+AHI).
%
  if ( n == 1 )

    a(1,1) = 0.5 * ( alo + ahi );

  else

    a(1:n,1) = ( (n-1:-1:0) * alo + (0:n-1) * ahi ) / ( n - 1 );

  end

  return
end
function report ( step, step_num, n, x, t, u, stability )

%*****************************************************************************80
%
%% REPORT prints or plots or saves the data at the current time step.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer STEP, the index of the current step,
%    between 0 and STEP_NUM.
%
%    Input, integer STEP_NUM, the number of steps to take.
%
%    Input, integer N, the number of nodes.
%
%    Input, real X(N), the coordinates of the nodes.
%
%    Input, real T, the current time.
%
%    Input, real U(N), the initial values U(X,T).
%
%    Input, real STABILITY, the stability factor, which should be
%    no greater than 1.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  STEP = %d\n', step );
  fprintf ( 1, '  TIME = %f\n', t );
  fprintf ( 1, '  STABILTY = %f\n', stability )
  fprintf ( 1, '\n' );
  for ilo = 1 : 5 : n
    ihi = min ( ilo + 4, n );
    for i = ilo : ihi
      fprintf ( 1, '  %14g', u(i) );
    end
    fprintf ( 1, '\n' );
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
function ua = u_a ( x, t )

%*****************************************************************************80
%
%% U_A sets the boundary condition for U at A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, T, the position and time.
%
%    Output, real UA, the prescribed value of U(X,T).
%
  ua = + 0.5;

  return
end
function ub = u_b ( x, t )

%*****************************************************************************80
%
%% U_B sets the boundary condition for U at B.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, T, the position and time.
%
%    Output, real UB, the prescribed value of U(X,T).
%
  ub = - 0.5;

  return
end
function u = u_init ( n, x, t )

%*****************************************************************************80
%
%% U_INIT sets the initial condition for U.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of nodes.
%
%    Input, real X(N), the coordinates of the nodes.
%
%    Input, real T, the current time.
%
%    Output, real U(N), the initial values U(X,T).
%
  ua = u_a ( x(1,1), t );
  ub = u_b ( x(n,1), t );

  q = 2.0 * ( ua - ub ) / pi;
  r = ( ua + ub ) / 2.0;
%
%  S can be varied.  It is the slope of the initial condition at the midpoint.
%
  s = 1.0;
  u(1:n,1) = ( 2 * x(1:n,1) - x(n,1) - x(1,1) ) ...
           / (                x(n,1) - x(1,1) );

  u(1:n,1) = - q * atan ( s * u(1:n,1) ) + r;

  return
end


text:
function a = imag(a)
%IMAG         imaginary part of taylor
%
%   c = imag(a)
%
%

% written  05/21/09     S.M. Rump
%

  a.t = imag(a.t);


text:
function [ishappy, cutoff] = strictCheck(f, values, data, pref)
%STRICTCHECK   Attempt to trim trailing Chebyshev coefficients in a
%   CHEBTECH. [ISHAPPY, CUTOFF] = STRICTCHECK(F, VALUES, DATA) returns an
%   estimated location CUTOFF at which the CHEBTECH F could be truncated to
%   maintain an accuracy of the default CHEBTECH CHEBFUNEPS preference
%   relative to DATA.VSCALE and DATA.HSCALE. ISHAPPY is TRUE if
%   CUTOFF < MIN(LENGTH(F.COEFFS), 2) or VSCALE(F)=0, and FALSE otherwise.
%
%   [ISHAPPY, CUTOFF] = STRICTCHECK(F, VALUES, DATA, PREF) allows
%   additional preferences to be passed. In particular, one can adjust the
%   target accuracy with PREF.CHEBFUNEPS. The VALUES field is ignored, but
%   included for consistency with other happiness checks.
%
%   STRICTCHECK tests to see if the absolute values of the entries in the
%   tail of coeffs, i.e., f.coeffs(1:TESTLENGTH,:), where
%       TESTLENGTH = n,             for n = 1:4
%       TESTLENGTH = 5,             for n = 5:44
%       TESTLENGTH = round((n-1)/8) for n > 44
%   all lie below the value in PREF.CHEBFUNEPS.  CUTOFF is the location of
%   the first entry above PREF.CHEBFUNEPS in absolute value.
%
%   STRICKCHECK differs from CLASSICCHECK() in that the tolerance
%   PREF.CHEBFUNEPS is not relaxed by the length of the representation of F
%   or by any finite difference approximation of the gradient of F.
%
% See also STRICTCHECK, LOOSECHECK.

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

% Deal with special cases ------------------------------------------------------

% Determine n (the length of the input):
n = length(f);

% Assume we're not happy:
ishappy = false;

% Grab some preferences:
if ( nargin == 1 )
    pref = f.techPref();
    epslevel = pref.chebfuneps;
elseif ( isnumeric(pref) )
    epslevel = pref;
else
    epslevel = pref.chebfuneps;
end

% Convert scalar epslevel/tolerance inputs into vectors.
if ( isscalar(epslevel) )
    epslevel = repmat(epslevel, 1, size(f.coeffs, 2));
end

% Deal with the trivial case:
if ( n < 2 )  % (Can't be simpler than a constant.)
    cutoff = n;
    return
end

if ( max(data.vscale) == 0 )
    % This is the zero function, so we must be happy.
    ishappy = true;
    cutoff = 1;
    return
elseif ( any(isinf(data.vscale(:))) )
    % Inf located. No cutoff.
    cutoff = n;
    return
end

% NaNs are not allowed.
if ( any(isnan(f.coeffs(:))) )
    error('CHEBFUN:CHEBTECH:strictCheck:nanEval', ...
        'Function returned NaN when evaluated.')
end

% Check for convergence and chop location --------------------------------------
testLength = min(n, max(5, round((n-1)/8)));
ac = bsxfun(@rdivide, abs(f.coeffs), data.vscale);
f.coeffs(bsxfun(@le, ac, epslevel)) = 0;
tail = f.coeffs(end-testLength+1:end,:);
if ( ~any(tail(:)) )
    cutoff = find(max(f.coeffs, [], 2) > 0, 1, 'last');
    ishappy = true;
else
    cutoff = n;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function pers=categorical_persistence(S)
% CATEGORICAL_PERSISTENCE computes the persistence of an unordered multilayer partition
%
% Version: 2.2.0
% Date: Thu 11 Jul 2019 12:25:42 CEST
%
% pers = categorical_PERSISTENCE(S) with a single multilayer partition or a
% cell of multilayer partitions S (with S{i} the ith multilayer partition of
% cell S) computes the "persistence" of multilayer partitions. The output pers
% is a vector such that pers(i) is the persistence of multilayer partition
% S{i} (stored as an N*T matrix where N is the number of nodes in each layer
% and T is the number of layers). For a multilayer network with unordered
% layers, the value of persistence is the sum over all pairs of layers of
% the number of nodes that do not change community assignments. Categorical
% persistence varies between 0 and 1. Categorical persistence is related to
% the multilayer quality function developed in Mucha et al. 2010 when one
% uses categorical interlayer coupling.
%
%   References:
%
%     Mucha, Peter J., Thomas Richardson, Kevin Macon, Mason A. Porter, and
%     Jukka-Pekka Onnela. "Community Structure in Time-Dependent,
%     Multiscale, and Multiplex Networks," Science 328, 876-878 (2010).
%
%     Bazzi, Marya, Mason A. Porter, Stacy Williams, Mark McDonald, Daniel
%     J. Fenn, and Sam D. Howison. "Community Detection in Temporal
%     Multilayer Networks, with an Application to Correlation Networks",
%     MMS: A SIAM Interdisciplinary Journal 14, 1-41 (2016).

if ~iscell(S)
    S={S};
end

[N,T]=size(S{1});
all2all = N*[(-T+1):-1,1:(T-1)];
A=spdiags(ones(N*T,2*T-2),all2all,N*T,N*T);
pers=zeros(length(S),1);
for i=1:length(S)

    G=sparse(1:length(S{i}(:)),S{i}(:),1);
    pers(i)=trace(G'*A*G)/(N*T*(T-1));
end

end

author: GenLouvain

repo: GenLouvain


text:
%Script that explains how to combine computeGTOM output with a hierarchical
%clustering
%   Joaquin Goñi <jgoni@unav.es> & Iñigo Martincorena
%   <imartincore@alumni.unav.es>
%   University of Navarra - Dpt. of Physics and Applied Mathematics &
%   Centre for Applied Medical Research.  Pamplona (Spain).
%
%   November 22nd 2007
%
load adjExample.mat    %example of an adjacency matrix
step = 3;   %value to compute GTOMstep (GTOM0,GTOM1,GTOM2, etc)
[GTOM]=computeGTOM(adjExample,step);
dist = ones(size(adjExample)) - GTOM;  %GTOM dissimilarity matrix
distVector = squareform(dist,'tovector');   %conversion to vector taking upper triangular values
Z = linkage(distVector,'average');  %linkage with average criteria
subplot(2,1,1), [H,T,PERM] = dendrogram(Z,0,'colorthreshold',0.35); %dendrogram plot with colorthreshold set to 0.35
GTOMOrdered = GTOM(PERM,PERM); %GTOM ordered by dendrogram
subplot(2,1,2), imagesc(GTOMOrdered); axis square; %plot of GTOMordered to check modules detection
%Automatic figure resizing and better merging to be done


text:
%% Affine invariant feature-based image matching
%
% This sample is similar to |feature_homography_demo.m|, but uses the affine
% transformation space sampling technique, called
% <http://www.ipol.im/pub/algo/my_affine_sift/ ASIFT>. While the original
% implementation is based on SIFT, you can try to use SURF or ORB detectors
% instead. Homography RANSAC is used to reject outliers.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/asift.py>
%

function asift_demo()
    % load input grayscale images
    files = {
        fullfile(mexopencv.root(),'test','aero1.jpg')
        fullfile(mexopencv.root(),'test','aero3.jpg')
    };
    for i=1:numel(files)
        if exist(files{i}, 'file') ~= 2
            disp('Downloading...')
            url = 'https://cdn.rawgit.com/opencv/opencv/3.2.0/samples/data/';
            [~,fname,ext] = fileparts(files{i});
            urlwrite([url fname ext], files{i});
        end
    end
    img1 = cv.imread(files{1}, 'Grayscale',true);
    img2 = cv.imread(files{2}, 'Grayscale',true);

    % ASIFT
    [vis, H] = run_asift(img1, img2, 'BRISK', true);
    imshow(vis), title('ASIFT')

    % interactive exploration
    if ~mexopencv.isOctave() && mexopencv.require('images')
        % interactively select region in img1
        pos1 = select_poly(gca) - 1;
        if isempty(pos1), return; end

        % apply homography
        pos2 = cv.perspectiveTransform(pos1, H);
        pos2(:,1) = pos2(:,1) + size(img1,2);

        % draw region in img1 and transformed region in img2
        vis = cv.polylines(vis, pos1, 'Closed',true, ...
            'Color',[0 255 0], 'Thickness',2, 'LineType','AA');
        vis = cv.polylines(vis, pos2, 'Closed',true, ...
            'Color',[0 255 0], 'Thickness',2, 'LineType','AA');
        imshow(vis), title('ASIFT')
    end
end

function [vis, H] = run_asift(img1, img2, feature_name, use_flann)
    %RUN_ASIFT  Run ASIFT algorithm, estimate homography and visualize matches

    % create detector and matcher objects
    [detector, matcher] = init_feature(feature_name, use_flann);

    % detect features using ASIFT method
    disp('Detecting...')
    tic, [kp1, desc1] = affine_detect(detector, img1); toc
    tic, [kp2, desc2] = affine_detect(detector, img2); toc
    fprintf('img1: %d features, img2: %d features\n', numel(kp1), numel(kp2));

    % match ASIFT features
    disp('Matching...')
    tic, matches = matcher.knnMatch(desc1, desc2, 2); toc
    fprintf('%d matches\n', numel(matches));

    % filter matches
    [matches, p1, p2] = filter_matches(kp1, kp2, matches);
    fprintf('%d good matches\n', numel(matches));
    assert(numel(matches) >= 4, 'not enough matches for homography estimation');

    % estimate homography with RANSAC method
    [H,inliers] = cv.findHomography(p1, p2, 'Method','Ransac');
    assert(~isempty(H), 'homography estimation failed');
    inliers = logical(inliers);
    fprintf('%d inliers\n', nnz(inliers));

    % visualize good and inlier matches
    vis = cv.drawMatches(img1, kp1, img2, kp2, matches, ...
        'NotDrawSinglePoints',true, 'MatchesMask',inliers);
end

function [detector, matcher] = init_feature(feature_name, use_flann)
    %INIT_FEATURE  Create detector and matcher objects

    % detector
    switch upper(feature_name)
        case 'SURF'
            detector = cv.SURF('HessianThreshold',400);
        case 'SIFT'
            detector = cv.SIFT();
        case 'ORB'
            detector = cv.ORB();
        case 'BRISK'
            detector = cv.BRISK();
        case 'AKAZE'
            detector = cv.AKAZE();
        case 'KAZE'
            detector = cv.KAZE();
        otherwise
            error('unrecognized feature: %s', feature_name)
    end

    % matcher
    if use_flann
        if ~isempty(strfind(detector.defaultNorm(), 'Hamming'))
            opts = {'LSH', 'TableNumber',6, 'KeySize',12, 'MultiProbeLevel',1};
        else
            opts = {'KDTree', 'Trees',5};
        end
        matcher = cv.DescriptorMatcher('FlannBasedMatcher', 'Index',opts);
    else
        matcher = cv.DescriptorMatcher('BFMatcher', ...
            'NormType',detector.defaultNorm());
    end
end

function [kpts, descs] = affine_detect(detector, img, mask)
    %AFFINE_DETECT  Affine-SIFT (ASIFT) algorithm
    %
    %     [kpts, descs] = affine_detect(detector, img)
    %     [kpts, descs] = affine_detect(detector, img, mask)
    %
    % ## Input
    % * __detector__ detector object
    % * __img__ input image
    % * __mask__ optional mask, default all image included
    %
    % ## Output
    % * __kpts__ concatenated keypoints
    % * __descs__ corresponding concatenated descriptors
    %
    % Applies a set of affine transormations to the image, detect keypoints,
    % and reproject them into initial image coordinates.
    % See: http://www.ipol.im/pub/algo/my_affine_sift/ for the details.
    %

    % default mask if not specified
    if nargin < 3
        mask = 255 * ones(size(img,1), size(img,2), 'uint8');
    end

    % parameter sampling: Tilt and Phi
    % (simulates all possible affine distortions caused by the change of
    % camera optical axis orientation from a frontal position)
    params = [1 0];  % no tilt and no rotation, thus original image
    for t = sqrt(2).^(1:5)       % geometric series 1, a, a^2, .., a^n; a=sqrt(2), n=5
        for p = 0:72/t:180-72/t  % arithmetic series 0, b/t, .., k*b/t < 180; b=72
            params(end+1,:) = [t, p];
        end
    end

    % for each pair of distortion parameters
    % (this loop is candidate for parallelization, i.e parfor)
    N = size(params,1);
    kpts = cell(N, 1);
    descs = cell(N, 1);
    for i=1:N
        % transform image
        [timg, tmask, Ai] = affine_skew(params(i,1), params(i,2), img, mask);

        % detect features using a similarity invariant matching method (SIFT)
        [kp, feat] = detector.detectAndCompute(timg, 'Mask',tmask);

        %TODO: Remove keypoints close to the boundary of the transformed image

        % project keypoints from the coordinates of the rotated and tilted
        % image back to the original image corrdinates
        for j=1:numel(kp)
            kp(j).pt = [kp(j).pt, 1] * Ai.';
        end
        kpts{i} = kp(:);
        descs{i} = feat;
    end

    % concatenate all features from all simulated images
    kpts = cat(1, kpts{:});
    descs = cat(1, descs{:});
end

function [img, mask, Ai] = affine_skew(tilt, phi, img, mask)
    %AFFINE_SKEW  Transform image/mask by an affine distortion
    %
    %     [img, mask, Ai] = affine_skew(tilt, phi, img, mask)
    %
    % ## Input
    % * __tilt__ tilt
    % * __phi__ rotation angle in degrees
    % * __img__ image
    % * __mask__ mask
    %
    % ## Output
    % * __img__ transformed image
    % * __mask__ transformed mask
    % * __Ai__ affine transform matrix from output/skewed `img` to input `img`
    %

    % initialize affine transformation matrix (identity)
    A = eye(2,3);
    [h,w,~] = size(img);

    % in the case of (tilt=1, phi=0), then no transformation (identity)

    if phi ~= 0
        % rotate image
        A = [cosd(phi) -sind(phi); sind(phi) cosd(phi)];
        corners = [0 0; w 0; w h; 0 h];
        corners = round(corners * A.');
        rect = cv.boundingRect(corners);
        A(:,3) = -rect(1:2);
        img = cv.warpAffine(img, A, 'DSize',rect(3:4), ...
            'Interpolation','Linear', 'BorderType','Replicate');
    end

    if tilt ~= 1
        % anti-aliasing filtering in the vertical direction, then tilt image
        % (subsample in vertical direction by a factor of tilt)
        s = 0.8 * sqrt(tilt^2 - 1);
        img = cv.GaussianBlur(img, 'KSize',[0 0], 'SigmaX',s, 'SigmaY',0.01);
        img = cv.resize(img, 1/tilt, 1, 'Interpolation','Nearest');
        A(1) = A(1) / tilt;
    end

    if phi ~= 0 || tilt ~= 1
        % apply same transformation on mask
        [h,w,~] = size(img);
        mask = cv.warpAffine(mask, A, 'DSize',[w h], 'Interpolation','Nearest');
    end

    % inverse affine transformation
    Ai = cv.invertAffineTransform(A);
end

function [matches, p1, p2] = filter_matches(kp1, kp2, matches, ratio)
    %FILTER_MATCHES  Filter out 2-NN matches using Lowe's ratio test

    if nargin < 4, ratio = 0.75; end

    % good matches
    idx = cellfun(@(m) (numel(m) == 2) && ...
        (m(1).distance < ratio * m(2).distance), matches);
    matches = cellfun(@(m) m(1), matches(idx));

    % corresponding points
    p1 = cat(1, kp1([matches.queryIdx]+1).pt);
    p2 = cat(1, kp2([matches.trainIdx]+1).pt);
end

function pos = select_poly(ax)
    %SELECT_POLY  Select polygon area using the mouse
    %
    %     pos = select_poly()
    %     pos = select_poly(ax)
    %
    % ## Input
    % * __ax__ axes handle, default gca
    %
    % ## Output
    % * __pos__ Nx2 matrix of points
    %

    if nargin < 1, ax = gca; end
    hRoi = impoly(ax);
    api = iptgetapi(hRoi);
    try
        api.setColor('green');
        api.setPositionConstraintFcn(makeConstrainToRectFcn('impoly', ...
            get(ax,'XLim').*[1 0.5], get(ax,'YLim')));
        pos = wait(hRoi);
    catch
        pos = zeros(0,2);
    end
    delete(hRoi);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function plotColors(img)
%
%
%        plotColors(img)
%
%        This function visualizes colors of the image in its 3D color
%        space.
%
%        Input:
%           -img: an image
%
%     Copyright (C) 2015  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/>.
%

check3Color(img);

[r, c, col] = size(img);

x = img(:,:,1);
y = img(:,:,2);
z = img(:,:,3);

c = reshape(img, r * c, col);

scatter3(x(:), y(:), z(:), [], c, 'filled');

end

author: banterle


text:
% Gabor filters with ERBs between points from a scale
%
% Description
% Gabor filters are complex, mostly analytic filters that have a
% Gaussian envelope in both the time and frequency domains. They are
% defined as
%
% .. math::
%
%      f(t) = C \sigma^{-1/2} \pi^{-1/4}
%             e^{\frac{-t^2}{2\sigma^2} + i\xi t}
%
% in the time domain and
%
% .. math::
%
%      \widehat{f}(\omega) = C \sqrt{2\sigma} \pi^{1/4}
%                            e^{\frac{-\sigma^2(\xi - \omega)^2}{2}}
%
% in the frequency domain. Though Gaussians never truly reach 0, in
% either domain, they are effectively compactly supported. Gabor
% filters are optimal with respect to their time-bandwidth product.
%
% `scaling_function` is used to split up the frequencies between
% `high_hz` and `low_hz` into a series of filters. Every subsequent
% filter's width is scaled such that, if the filters are all of the
% same height, the intersection with the precedent filter's response
% matches the filter's Equivalent Rectangular Bandwidth (erb == True)
% or its 3dB bandwidths (erb == False). The ERB is the width of a
% rectangular filter with the same height as the filter's maximum
% frequency response that has the same L_2 norm.
%
% Properties
% centers_hz : [2-D Array]
% is_real : [bool]
% is_analytic : [bool]
% num_filts : [int]
% sampling_rate : [double]
% supports_hz : [double] 2-D Array
% supports : [double] 2-D Array
% supports_ms : [double] 2-D Array
% scaled_l2_norm : [bool]
% erb : [bool]
%
% See Also
% Config.EFFECTIVE_SUPPORT_THRESHOLD : the absolute
%     value below which counts as zero
%
% Copyright (c) 2018 Department of Computer Science,
%                    University of Toronto, Canada,
%                    Vector Institute, Canada
%
% License
% This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Yingxue Wang <yingxue@cs.toronto.edu>
%  Sean Robertson <sdrobert@cs.toronto.ca>
%

classdef GaborFilterBank < LinearFilterBank

    properties
        centers_hz
        is_real
        is_analytic = false
        num_filts = 40
        sampling_rate = 16000
        supports_hz
        supports
        supports_ms
    end

    properties (Access=private)
        scaled_l2_norm = false
        erb = false
        wrap_below
        centers_ang
        stds
        supports_ang
        wrap_supports_ang
        scale_l2_norm
    end

    methods

        function obj = GaborFilterBank(scaling_function,...
                num_filts,...
                low_hz,...
                sampling_rate,...
                scale_l2_norm,...
                erb,...
                high_hz)

            % Constructor for GaborFilterBank
            %
            % Input
            % scaling_function : [filterbanks/scales/ScalingFunction]
            %     Dictates the layout of filters in the Fourier domain. Can be
            %     a ScalingFunction
            % num_filts : [int]
            %     The number of filters in the bank
            % high_hz, low_hz : [double], optional
            %     The topmost and bottommost edge of the filters, respectively.
            %     The default for high_hz is the Nyquist
            % sampling_rate : [double], optional
            %     The sampling rate (cycles/sec) of the target recordings
            % scale_l2_norm : [bool]
            %     Whether to scale the l2 norm of each filter to 1. Otherwise the
            %     frequency response of each filter will max out at an absolute
            %     value of 1.
            % erb : [bool]
            %
            % Output
            % obj : [GaborFilterBank object]
            %
            % Example
            % >>scaling_function = MelScaling();
            % >>num_filts = 11;
            % >>low_hz = 0;
            % >>sampling_rate = 8000;
            % >>scale_l2_norm = false;
            % >>erb = false;
            % >>bank = GaborFilterBank(scaling_function,num_filts,...
            %                low_hz, sampling_rate, scale_l2_norm, erb);

            if nargin == 1
                num_filts = 40;
                low_hz = double(20);
                sampling_rate = double(16000);
                scale_l2_norm = false;
                erb = false;
                high_hz = double.empty;
            elseif nargin == 2
                low_hz = double(20);
                sampling_rate = double(16000);
                scale_l2_norm = false;
                erb = false;
                high_hz = double.empty;
            elseif nargin == 3
                sampling_rate = double(16000);
                scale_l2_norm = false;
                erb = false;
                high_hz = double.empty;
            elseif nargin == 4
                scale_l2_norm = false;
                erb = false;
                high_hz = double.empty;
            elseif nargin == 5
                erb = false;
                high_hz = double.empty;
            elseif nargin == 6
                high_hz = double.empty;
            end

            obj = obj@LinearFilterBank();
            obj.scaled_l2_norm = scale_l2_norm;
            obj.erb = erb;
            obj.sampling_rate = sampling_rate;
            if isempty(high_hz)
               high_hz = floor(sampling_rate/2);
            end

            % Input validation
            if low_hz < 0 || high_hz <= low_hz || high_hz > floor(sampling_rate/2)
                error('Invalid frequency range: (%d, %d)', low_hz, high_hz);
            end
            scale_low = scaling_function.hertz_to_scale(low_hz);
            scale_high = scaling_function.hertz_to_scale(high_hz);
            scale_delta = (scale_high - scale_low) / (num_filts + 1);
            % edges dictate the points where filters should intersect. We
            % make a pretend intersection halfway between low_hz and
            % the first filter center in the scaled domain. Likewise with
            % high_hz and the last filter center. Intersections are spaced
            % uniformly in the scaled domain
            edges = [];
            for idx= 1:num_filts+1
                edges = [edges ; scaling_function.scale_to_hertz(scale_low + scale_delta * ((idx-1) + .5))];
            end
            centers_hz = [];
            centers_ang = [];
            stds = [];
            supports_ang = [];
            supports = [];
            wrap_supports_ang = [];
            obj.wrap_below = false;
            log_2 = log(2);
            log_pi = log(pi);
            t_support_const = -2 * log(Config.EFFECTIVE_SUPPORT_THRESHOLD);
            f_support_const = t_support_const;

            if scale_l2_norm
               f_support_const = f_support_const + log_2 + .5 * log_pi;
               t_support_const = t_support_const - .5 * log_pi;
            else
               t_support_const = t_support_const - log_2 - log_pi;
            end
            if erb
               bandwidth_const = sqrt(pi) / 2;
            else
               bandwidth_const = sqrt(3 / 10 * log(10));
            end

            zip_edges = [edges(1:end-1) edges(2:end)];

            for ind = 1:size(zip_edges,1)  % 1st dimension of zip_edges
                left_intersect = zip_edges(ind);
                right_intersect = zip_edges(ind + size(zip_edges,1));
                center_hz = (left_intersect + right_intersect)/2;
                center_ang = Util.hertz_to_angular(center_hz, obj.sampling_rate);
                std = bandwidth_const / Util.hertz_to_angular(center_hz - left_intersect, obj.sampling_rate);
                log_std = log(std);

                if scale_l2_norm
                    diff_ang = sqrt(log_std + f_support_const) / std;
                    wrap_diff_ang = sqrt(log_std + f_support_const + log_2) / std;
                    diff_samps = ceil(std * sqrt(t_support_const - log_std));
                else
                    diff_ang = sqrt(f_support_const) / std;
                    wrap_diff_ang = sqrt(f_support_const + log_2) / std;
                    diff_samps = ceil(std * sqrt(t_support_const - 2 * log_std));
                end
                supp_ang_low = center_ang - diff_ang;
                if supp_ang_low < 0
                    obj.wrap_below = true;
                end
                centers_hz = [centers_hz; center_hz];
                centers_ang = [centers_ang; center_ang];
                supports_ang = [supports_ang; [center_ang - diff_ang center_ang + diff_ang]];
                wrap_supports_ang = [wrap_supports_ang; 2 * wrap_diff_ang];
                supports = [supports; [-diff_samps diff_samps]];
                stds = [stds; std];
            end

            obj.centers_ang = centers_ang;
            obj.centers_hz = centers_hz;
            obj.stds = stds;
            obj.supports_ang = supports_ang;
            obj.wrap_supports_ang = wrap_supports_ang;
            obj.num_filts = num_filts;

            supports_hz =[];
            for idx = 1:size(supports_ang,1) % 1st dimension of supports_ang
                supports_hz = [supports_hz; ...
                    [Util.angular_to_hertz(supports_ang(idx), obj.sampling_rate) ...
                    Util.angular_to_hertz(supports_ang(idx+size(supports_ang,1)), obj.sampling_rate)]];
            end
            obj.supports_hz = supports_hz;
            obj.supports = supports;
        end

    end

    methods
        function res = get_impulse_response (obj, filt_idx, width)
            % Helper function to get impulse response of a GaborFilterBank
            %
            % Inputs
            % filt_idx  : [int] filter index
            % width     : [int] width of the signal
            %
            % Outputs
            % res : [double]
            %
            % Example
            % >>scaling_function = MelScaling();
            % >>num_filts = 11;
            % >>low_hz = 0;
            % >>sampling_rate = 8000;
            % >>scale_l2_norm = false;
            % >>erb = false;
            % >>bank = GaborFilterBank(scaling_function,num_filts,...
            %                low_hz, sampling_rate, scale_l2_norm, erb);
            % >>x = bank.get_impulse_response(filt_idx, dft_size);
            %
            center_ang = obj.centers_ang(filt_idx);
            std = obj.stds(filt_idx);
            res = zeros(width,1); % we want it to store complex values
            if obj.scale_l2_norm
                const_term = -0.5 * log(std) - 0.25 * log(pi);
            else
                const_term = -0.5 * log(2 * pi) - log(std);
            end
            denom_term = 2 * std ^ 2;

            for t = 1: width+1
                val = complex(-(t-1) ^ 2 / denom_term + const_term, center_ang * (t-1));
                val = exp(val);
                if t ~= width+1
                    res(t) = res(t) + val;
                end
                if width-t+2 > 0 && t>1
                    res(end-t+2) = res(end-t+2) + conj(val);
                end
            end
        end

        function res = get_frequency_response(obj, filt_idx, width, half)
            % Helper function to get frequency response of a GaborFilterBank
            %
            % Inputs
            % filt_idx  : [int] filter index
            % width     : [int] width of the signal
            % half      : [bool] (def. false)
            %
            % Outputs
            % res : [hz] double
            %
            % Example
            % >>scaling_function = MelScaling();
            % >>num_filts = 11;
            % >>low_hz = 0;
            % >>sampling_rate = 8000;
            % >>scale_l2_norm = false;
            % >>erb = false;
            % >>bank = GaborFilterBank(scaling_function,num_filts,...
            %                low_hz, sampling_rate, scale_l2_norm, erb);
            % >>X = bank.get_frequency_response(filt_idx, dft_size, false);
            %
            if nargin == 3
                half = false;
            end
            center_ang = obj.centers_ang(filt_idx);
            lowest_ang = obj.supports_ang(filt_idx,1);
            highest_ang = obj.supports_ang(filt_idx,2);
            lowest_ang = filt_idx;
            std = obj.stds(filt_idx);
            dft_size = width;

            if half == true
                if mod(width,2) ~= 0
                    dft_size = floor((width + 1) / 2);
                else
                    dft_size = floor(width / 2) + 1;
                end
            end
            res = zeros(dft_size, 1);
            if obj.scale_l2_norm
                const_term = .5 * log(2 * std) + .25 * log(pi);
            else
                const_term = 0;
            end
            num_term = -(std ^ 2) / 2;

            for idx = 1: dft_size
                for period =  -1 - floor(max(-lowest_ang, 0) / (2 * pi)): 2 + floor(highest_ang / (2 * pi)) - 1
                    omega = ((double(idx)-1) / double(width) + double(period)) * 2 * pi;
                    val = num_term * (center_ang - omega) ^ 2 + const_term;
                    val = exp(val);
                    res(idx) = res(idx) + val;
                end
            end
        end

        function  [start_idx, truncated_response] = get_truncated_response(obj, filt_idx, width)
            % Helper function to get truncated response of a GaborFilterBank
            %
            % Desription
            % wrap_supports_ang contains the angular supports of each filter
            % if the effective support threshold were halved. If this
            % support exceeds the 2pi period, overlap from aliasing in the
            % periphery will exceed the effective support, meaning the
            % entire period lies in the support
            %
            % Inputs
            % filt_idx  : [int] filter index
            % width     : [int] width of the signal
            %
            % Outputs
            % start_idx : [int]
            % truncated_response : [1-D column vector]
            %
            % Example
            % >>scaling_function = MelScaling();
            % >>num_filts = 11;
            % >>low_hz = 0;
            % >>sampling_rate = 8000;
            % >>scale_l2_norm = false;
            % >>erb = false;
            % >>bank = GaborFilterBank(scaling_function,num_filts,...
            %                low_hz, sampling_rate, scale_l2_norm, erb);
            % [start_idx, truncated_response] = bank.get_truncated_response(filt_idx, dft_size);
            %
            % Please refer to test_fbank.test_truncated_matches_full
            % for more information
            %

            if obj.wrap_supports_ang(filt_idx) >= 2 * pi
                start_idx = 1;
                truncated_response = obj.get_frequency_response(filt_idx, width);
            else
                center_ang = obj.centers_ang(filt_idx);
                std = obj.stds(filt_idx);
                lowest_ang = obj.supports_ang(filt_idx, 1);
                highest_ang = obj.supports_ang(filt_idx, 2);
                left_idx = ceil(width * lowest_ang / (2 * pi));
                right_idx = round(width * highest_ang / (2 * pi));
                res = zeros(1 + right_idx - left_idx, 1);

                if obj.scale_l2_norm
                    const_term = .5 * log(2 * std) + .25 * log(pi);
                else
                    const_term = 0;
                end
                num_term = -(std ^2) / 2;

                for idx = left_idx: right_idx
                    for period = - round(max(-lowest_ang, 0) / (2 * pi)): 1 + round(highest_ang / (2 * pi)) - 1 % -1 because the syntax is different in matlab and python
                        omega = ((idx-1) / width + period) * 2 * pi;
                        val = num_term * (center_ang - omega) ^ 2 + const_term;
                        val = exp(val);
                        res(idx - left_idx + 1) = res(idx - left_idx + 1) + val; % Need to check if this is right
                    end
                end
                start_idx = mod(left_idx, width);
                truncated_response = res;
            end
        end
    end

end


text:
function x=ifmt(mellin,beta,M);
%IFMT	Inverse fast Mellin transform.
%	X=IFMT(MELLIN,BETA,M) computes the inverse fast Mellin
%	transform of MELLIN.
%	WARNING : the inverse of the Mellin transform is correct only
%	if the Mellin transform has been computed from FMIN to 0.5 Hz,
%	and if the original signal is analytic.
%
%	MELLIN : Mellin transform to be inverted. Mellin must have been
%	 obtained from FMT with frequency running from FMIN to 0.5 Hz.
%	BETA : Mellin variable issued from FMT.
%	M : number of points of the inverse Mellin transform.
%					(default : length(MELLIN)).
%	X : inverse Mellin transform with M points in time.
%
%	Example :
%	 sig=atoms(128,[64,0.25,32,1]);
%	 [MELLIN,BETA]=fmt(sig,0.05,0.5,256); clf;
%	 X=ifmt(MELLIN,BETA,128); plot(real(X)); hold;
%	 plot(real(sig),'g'); hold;
%
%	See also : fmt, fft, ifft.
%

%	P. Goncalves 9-95 - O. Lemoine, June 1996.
%	Copyright (c) 1995 Rice University - CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

N=length(mellin);

if nargin<2,
 error('At least 2 input parameters required');
elseif nargin==2,
 M=N;
end

No2 = (N+rem(N,2))/2;
q   = exp(1/(N*(beta(2)-beta(1))));
fmin = 0.5/(q^(No2-1));

% Inverse Mellin transform computation
p = 0:(N-1);
L = log(fmin)/log(q);
S = fft(fftshift(mellin.*exp(-j*2*pi*L*(p/N-1/2))/(N*log(q))));
S = S(1:No2);

% Inverse Fourier transform
k = (1:No2);
x = zeros(M,1);
t = (1:M)-(M+rem(M,2))/2-1;
geo_f = fmin*(exp((k-1).*log(q))) ;
itfmatx = zeros(M,No2);
itfmatx = exp(2*i*t'*geo_f*pi);
for k=1:M,
 x(k) = real(integ(itfmatx(k,:).*S,geo_f));
end;

x = hilbert(x);

% Normalization
SP = fft(x); fmax=0.5;
indmin = 1+round(fmin*(M-2));
indmax = 1+round(fmax*(M-2));
SPana = SP(indmin:indmax);
nu = (indmin:indmax)'/M;
SPp = SPana./nu;
Esm = SPp'*SPana;
x = x*norm(mellin)/sqrt(Esm);

author: HeLiangHIT


text:
function [datapt, unitStr] = lvd_ElectricalPowerSinkTasks(stateLogEntry, subTask, pwrSink)
%lvd_ElectricalPowerSinkTasks Summary of this function goes here
%   Detailed explanation goes here
    stageStates = stateLogEntry.stageStates;

    pwrSinkStates = AbstractLaunchVehicleElectricalPowerSnkState.empty(1,0);
    for(i=1:length(stageStates)) %#ok<*NO4LP>
        pwrSinkStates = horzcat(pwrSinkStates, stageStates(i).powerSinkStates); %#ok<AGROW>
    end

    pwrSinks = AbstractLaunchVehicleElectricalPowerSrcSnk.empty(1,0);
    for(i=1:length(pwrSinkStates))
        pwrSinks = horzcat(pwrSinks, pwrSinkStates(i).getEpsSinkComponent()); %#ok<AGROW>
    end

    pwrSinkInd = find(pwrSinks == pwrSink,1,'first');
    pwrSinkState = pwrSinkStates(pwrSinkInd);

    switch subTask
        case 'active'
            if(not(isempty(pwrSinkState)))
                pwrSinkState = pwrSinkState(1);
                datapt = pwrSinkState.getActiveState();

                if(isempty(datapt))
                    datapt = -1;
                end

                unitStr = '';
            else
                datapt = -1;
                unitStr = '';
            end

        case 'dischargeRate'
            if(not(isempty(pwrSinkState)))
                pwrSinkState = pwrSinkState(1);

                elemSet = stateLogEntry.getCartesianElementSetRepresentation();
                steeringModel = stateLogEntry.steeringModel;
                datapt = pwrSinkState.getElectricalPwrRate(elemSet, steeringModel, [], [], []);

                if(isempty(datapt))
                    datapt = -1;
                end

                unitStr = 'EC/s';
            else
                datapt = -1;
                unitStr = 'EC/s';
            end

        otherwise
            error('Unknow EPS sink task: %s', subTask);
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [groups, datasets] = hdf5dir(h5info, path, pos)
% [groups, datasets] = hdf5dir(h5info, path, pos)
%
% This function never returns an error - but will return the same thing for
% an empty group or failure to locate a node
%
% h5info - A Group in the hierarchy (send at least
% something.GroupHierarchy)
% path (optional) - the location of the hdf5 dataset or group
% pos (optional, mostly for recursion) - the path represented by h5info.

if nargin < 2
    path = '/';
end

if nargin < 3
    pos = '';
end

% The top level hdf5info doesn't drop you into the groups, so we check into
% that first
if isfield(h5info, 'GroupHierarchy')
    h5info = h5info.GroupHierarchy;
end

groups = {};
datasets = {};

[tok, path] = strtok(path, '/');

if isempty(tok)
    if length(h5info.Groups) > 0
        groups = h5info.Groups;
    end

    if length(h5info.Datasets) > 0
        datasets = h5info.Datasets;
    end
else
    pos = [pos, '/', tok];

    for i = 1:length(h5info.Groups)
        if strcmp(pos, h5info.Groups(i).Name)
            [groups, datasets] = hdf5dir(h5info.Groups(i), path, pos);
        end
    end

    for i = 1:length(h5info.Datasets)
        if strcmp(pos, h5info.Datasets(i).Name)
            datasets = {h5info.Datasets(i).Name};
        end
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% written by Ali A. Eftekhari
clc
L = 50;  % domain length
Nx = 20; % number of cells
m = createMesh3D(Nx, Nx, Nx, L, L, L);
BC = createBC(m); % all Neumann boundary condition structure
BC.left.a(:) = 0; BC.left.b(:)=1; BC.left.c(:)=1; % Dirichlet for the left boundary
BC.right.a(:) = 0; BC.right.b(:)=1; BC.right.c(:)=0; % right boundary
D_val = 1; % value of the diffusion coefficient
D = createCellVariable(m, D_val); % assign the diffusion coefficient to the cells
D_face = harmonicMean(D); % calculate harmonic average of the diffusion coef on the cell faces
Mdiff = diffusionTerm(D_face); % matrix of coefficients for the diffusion term
[Mbc, RHSbc] = boundaryCondition(BC); % matix of coefficients and RHS vector for the BC
M = Mdiff + Mbc; % matrix of cefficients for the PDE
c = solvePDE(m,M, RHSbc); % send M and RHS to the solver
visualizeCells(c); % visualize the results

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
%% Optimal Kernel Selection
%
%%
% In the section <DensityEstimation.html density estimation> we have seen
% that the correct choice of the kernel halfwidth is essential for creating a good
% match between the true density function and the reconstructed density
% function. If the halfwidth is set too small the reconstructed density
% function is usually oscillating and the indiviudual sampling points are
% visible as sharp peaks. If the halfwidth is too large the resulting
% density function is usually too smooth and does not reproduce the
% features of the original density function.
%
% Finding an optimal kernel halfwidth is a hard problem as the optimal
% kernel halfwidth depends not only on the number of sampling points but also
% on the smoothness of the true but unknown density function.
% MTEX offers several options set by flags during the kernel calculation operation.  A very
% conserative choice for the kernel halfwidth that takes into account only
% the number of sampling points is implemented in MTEX with the flag |'magicRule'|. The flag
% |'RuleOfThumb'| considers both the number of sampling
% points and the variance of the sampling points as an estimate of the
% smoothness of the true density function. The most advanced (and default)
% method for estimating the optimal kernel halfwidth is
% <orientation.KLCV.html Kullback Leibler cross validation>.
% This method tests different kernel halfwidths on a subset of the
% random sample and selects the halfwidth which best reproduces the
% ommited points of the random sample.
%
% In order to demonstrate this functionality let's start with the following
% orientation density function

% Define trigonal crystal symmetry using Enantiomorphic Point Group notation
cs = crystalSymmetry('32');

% Build a density function by combining a uniform texture with two pre-defined texture components
odf = 0.25*uniformODF(cs) + 0.25*unimodalODF(orientation.brass(cs)) + ...
  0.5*fibreODF(fibre.alpha(cs),'halfwidth',10*degree);

% plot the density function as six sigma sections
plot(odf,'sections',6,'silent','sigma')
mtexColorbar

%%
% and compute $10000$ random orientations representing this density function using the command
% |<orientation.discreteSample.html discreteSample>|

ori = odf.discreteSample(10000)

%%
% Next we estimate the optimal <ODFShapes.html kernel function> using the
% command |<orientation.calcKernel.html calcKernel>| with the default settings.

psi  = calcKernel(ori)

%%
% This kernel can now be used to reconstruct the original ODF from the sampled points using the command
% <DensityEsimation.html density estimation>

odf_rec = calcDensity(ori,'kernel',psi)

% plot the reconstructed ODF and compare it to the plot of the original function.  The results are similar but not identical.
figure;plot(odf_rec,'sections',6,'silent','sigma')
mtexColorbar

%% Exploration of the relationship between estimation error and number of single orientations
%
% In this section we want to compare the different methods for estimating
% the optimal kernel halfwidth. To this end we simulate 10, 100, ...,
% 1000000 single orientations from the model ODF |odf|, compute optimal
% kernels according to the |'magicRule'|, the |'RuleOfThumb'| and
% <orientation.KLCV.html Kullback Leibler cross validation> and then
% compute the fit between the reconstructed |odf_rec| and the original
% |odf|.

% define a variable to hold the calculated error values
e = [];
for i = 1:6

  % calculate a sample of orientations from the model ODF
  ori = discreteSample(odf,10^i,'silent');

  % calculate the kernel using the function defaults, reconstruct the odf, and calculate error between this and the original ODF
  psi1 = calcKernel(ori,'SamplingSize',10000,'silent');
  odf_rec = calcDensity(ori,'kernel',psi1,'silent');
  e(i,1) = calcError(odf_rec,odf,'resolution',2.5*degree);

  % calculate the kernel using the RuleOfThumb, reconstruct the odf, and calculate error between this and the original ODF
  psi2 = calcKernel(ori,'method','RuleOfThumb','silent');
  odf_rec = calcDensity(ori,'kernel',psi2,'silent');
  e(i,2) = calcError(odf_rec,odf,'resolution',2.5*degree);

  % calculate the kernel using the magicRule, reconstruct the odf, and calculate error between this and the original ODF
  psi3 = calcKernel(ori,'method','magicRule','silent');
  odf_rec = calcDensity(ori,'kernel',psi3,'silent');
  e(i,3) = calcError(odf_rec,odf,'resolution',2.5*degree);

  % generate text showing the kernel size calculated with each method in each loop
  disp(['RuleOfThumb: ' int2str(psi2.halfwidth/degree) mtexdegchar ...
    ' KLCV: ' int2str(psi1.halfwidth/degree) mtexdegchar ...
    ' magicRule: ' int2str(psi3.halfwidth/degree) mtexdegchar ...
    ]);

end

%%
% Plot the error to the number of single orientations sampled from the original ODF.

close all;
loglog(10.^(1:length(e)),e,'LineWidth',2)
legend('Default','RuleOfThumb','magicRule')
xlabel('Number of orientations (log scale)')
ylabel('Estimation Error in degrees')
title('Error between original ODF model and the reconstructed ODF','FontWeight','bold')

author: mtex-toolbox

repo: mtex

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


text:
function camera_info = readCameraInfo(filename)

file_words=file2cellarray(filename);

camera_info = [];

for iline=1:length(file_words)
    switch(file_words{iline}{1})
        case('camera_number:')
            fprintf('camera_number: %s\n',file_words{iline}{2});
        case('camera_index:')
            camera_info(end+1).camera_index = str2num(file_words{iline}{2});
        case('')
        otherwise
            value = str2double(file_words{iline}{2});
            if isnan(value)
                value = file_words{iline}{2};
            end
            camera_info(end).(file_words{iline}{1}(1:end-1)) = value;
    end
end

author: jianxiongxiao

repo: ProfXkit


text:
function c = sum(a,dim)
%SUM          Implements  sum(a,dim)  for slopes
%
%   u = sum(a,dim)
%
%parameter dim optional, functionality as Matlab function sum
%

% written  12/06/98     S.M. Rump
% modified 09/28/01     S.M. Rump  matrices and multi-dimensional arrays
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if nargin==1
    if a.size(1)==1
      dim = 2;
    else
      dim = 1;
    end
  end

  if dim==1
    c = ones(1,a.size(1))*a;
  else
    c = a*ones(a.size(2),1);
  end

  setround(rndold)


text:
%compute dratioarealength

function [data,units]=compute_dratioarealength(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
dratioarealength=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    dratioarealength{1,i}=(trx(larva).ratioarealength(2:end)-trx(larva).ratioarealength(1:end-1))./trx(larva).dt;
end

units=parseunits('mm/s');
data=dratioarealength;

author: kristinbranson

repo: JAABA


text:
function e = mgVcyclefracLap1d(A,r,pde,square,h,s,option)

global mA

%% Setting
if ~isfield('option','smootherType')
    smootherType = option.smootherType;
else
    smootherType = 3;
end
if ~isfield('option','smootherstep')
    smoothingstep = option.smoothingstep;
else
    smoothingstep = 3;
end
if ~exist('smootherType','var'), smootherType = 3; end
% mapping parameter
if s == 0.5
    gamma = 1;
else
    gamma = 3/(2*s)+0.1;
end

%% line smoothing
N = length(r);
rold = r;
x0 = square(1); x1 = square(2);
y0 = square(3); y1 = square(4);
nx = (x1-x0)/h+1;  % number of grid points in x-direction
ny = (y1-y0)/h+1;  % number of grid points in x-direction
% 1D vector to 2D matrix
r2d = reshape(r,ny,nx);
e = zeros(N,1);
e2d = reshape(e,ny,nx);
% index for Red-Black smoothing
nodeidx = reshape(1:N,ny,nx);
redidx = 2:2:nx-1;
blackidx = 3:2:nx-2;
rlinear = nodeidx(:,redidx);
blinear = nodeidx(:,blackidx);
for it = 1:smoothingstep
    switch smootherType
        case 0
            e = e + tril(A)\r;
            r = rold - A*e;
            e2d = reshape(e,ny,nx);
        case 1
            for k = 2:nx-1 % only update interiori lines
                % form residual
                idx = (k-1)*ny+1:k*ny;  % the end nodes is included?
                de = A(idx,idx)\r2d(:,k);
                e2d(:,k) = e2d(:,k) + de;
                r2d(:,k) = 0;  % exact solve
                if k<nx-1   % update right column
                    r2d(:,k+1) = r2d(:,k+1) - A(idx+ny,idx)*de;
                end
                if k>2    % update left column
                % for left-to-right ordering, the left r(idx-ny) is zero
                    r2d(:,k-1) = - A(idx-ny,idx)*de;
                end
            end
        case 2 % weighted Jacobi
            for k = 2:nx-1
                % form residual
                idx = (k-1)*ny+1:k*ny;
                de = A(idx,idx)\r(idx);
                e(idx) = e(idx) + 0.25*de;
            end
            r = rold - A*e;
        case 3 % red-black block Gauss-Seidel
            idx = ny+1:2*ny;
            % 1: red lines
            de = A(idx,idx)\r2d(:,redidx);
            e2d(:,redidx) = e2d(:,redidx) + de;
            % update residual
            r2d(:,redidx) = 0;
            r2d(:,blackidx) = r2d(:,blackidx)-reshape(A(blinear,rlinear)*de(:),ny,length(blackidx));
            % previous r2d(:,blackidx) = 0;
            % 2: black lines
            de = A(idx,idx)\r2d(:,blackidx);
            e2d(:,blackidx) = e2d(:,blackidx) + de;
            % update residual
            r2d(:,blackidx) = 0;
            r2d(:,redidx) = -reshape(A(rlinear,blinear)*de(:),ny,length(redidx));
            % previous r2d(:,redidx) = 0;
    end
end

%% Transfer operator
% prolongation and restriction in x-direction
Ix = prolongation1d(log2(nx-1)-1);
Rx = Ix';
% prolongation and restriction in y-direction
nxc = (nx - 1)/2 + 1;
nyc = (ny - 1)/2 + 1;
% geometric quantity
My = ny - 1;
Tyf = ((0:My)'/My).^gamma*(y1-y0) + y0;
hf = diff(Tyf);
Myc = nyc - 1;
Tyc = ((0:Myc)'/Myc).^gamma*(y1-y0) + y0;
hc = diff(Tyc);
jc = 2:nyc-1;  % interiori points in the coarse grid
j = 2*jc-1;  % index of coarse points in the fine grids
alpha = zeros(1,nyc);
beta = zeros(1,nyc);
alpha(jc) = hf(j)./hc(jc);
beta(jc) = hf(j-1)./hc(jc-1);
alpha(1) = 1-beta(2);
jc = 2:(nyc-1);
jf = 2*jc-1;  % index of coarse points in the fine grids
% due to the Neumann boundary condition, 1 is included
ii = [1 jf 2 jf+1 jf-1];
jj = [1 jc 1 jc jc];
ss = [ones(1,nyc-1) alpha(1:nyc-1) beta(2:nyc-1)];
Iy = sparse(ii,jj,ss,ny,nyc);
Ry = Iy';

%% Restriction
% 1D vector to 2D matrix
% r2d = reshape(r,ny,nx);
rc2d = Ry*r2d*Ix;
rc = rc2d(:);

%% Coarse grid correction
% option.solver = 'none';
% [u,eqn] = fracLap1d(square,2*h,pde,option);
level = -log2(h);
Ac = mA{level-1};
fixedNode = [1:nyc (2:nxc-1)*nyc (nxc-1)*nyc+(1:nyc)];
isBdNode = false(length(rc),1);
isBdNode(fixedNode) = true;
freeNode = find(~isBdNode);
if level <= 3
    ec = zeros(length(rc),1);
    ec(freeNode) = Ac(freeNode,freeNode)\rc(freeNode);
else
    ec = mgVcyclefracLap1d(Ac,rc,pde,square,2*h,s,option);
end

%% Prolongation
% 1D vector to 2D matrix
ec2d = reshape(ec,nyc,nxc);
e2d = e2d + Iy*ec2d*Rx;
e = e2d(:);
% update residual when e updated. remember we are solving Ae = rold.
r = rold - A*(e2d(:));
r2d = reshape(r,ny,nx);

%% Post-smoothing
for it = 1:smoothingstep
    switch smootherType
        case 0 % pointwise G-S smoothing
            e = e + triu(A)\r;
            r = rold - A*e;
            e2d = reshape(e,ny,nx);
        case 1
            for k = nx-1:-1:2
                % form residual
                idx = (k-1)*ny+1:k*ny;
                de = A(idx,idx)\r2d(:,k);
                e2d(:,k) = e2d(:,k) + de;
                r2d(:,k) = 0;
                if k < nx-1   % update right column
                % for right-to-left ordering, the right r(idx+ny) is zero
                    r2d(:,k+1) = - A(idx+ny,idx)*de;
                end
                if k>2    % update previous column
                    r2d(:,k-1) = r2d(:,k-1) - A(idx-ny,idx)*de;
                end
            end
        case 2 % weighted Jacobi
            for k = 1:nx
                % form residual
                idx = (k-1)*ny+1:k*ny;
                de = A(idx,idx)\r(idx);
                e(idx) = e(idx) + 0.25*de;
            end
            r = rold - A*e;
        case 3 % red-black block Gauss-Seidel
            idx = ny+1:2*ny;
            % 1: red lines
            de = A(idx,idx)\r2d(:,redidx);
            e2d(:,redidx) = e2d(:,redidx) + de;
            % update residual
            r2d(:,redidx) = 0;
            r2d(:,blackidx) = r2d(:,blackidx)-reshape(A(blinear,rlinear)*de(:),ny,length(blackidx));
            % previous r2d(:,blackidx) = 0;
            % 2: black lines
            de = A(idx,idx)\r2d(:,blackidx);
            e2d(:,blackidx) = e2d(:,blackidx) + de;
            % update residual
            r2d(:,blackidx) = 0;
            r2d(:,redidx) = -reshape(A(rlinear,blinear)*de(:),ny,length(redidx));
            % previous r2d(:,redidx) = 0;
    end
end
e = e2d(:);

repo: ifem