text:
function [gci,goi] = v_dypsa(s,fs)
%V_DYPSA   Derive glottal closure instances from speech [gci,goi] = (s,fs)
%   Note: Needs to be combined with a voiced-voiceless detector to eliminate
%   spurious closures in unvoiced and silent regions.
%
%   Inputs:
%   s     is the speech signal
%   fs    is the sampling frequncy
%
%   Outputs:
%   gci   is a vector of glottal closure sample numbers
%   gco   is a vector of glottal opening sample numbers derived from
%         an assumed constant closed-phase fraction
%
%   References:
%   [1]  P. A. Naylor, A. Kounoudes, J. Gudnason, and M. Brookes, "Estimation of Glottal Closure
%        Instants in Voiced Speech using the DYPSA Algorithm," IEEE Trans on Speech and Audio
%        Processing, vol. 15, pp. 34-43, Jan. 2007.
%   [2]  M. Brookes, P. A. Naylor, and J. Gudnason, "A Quantitative Assessment of Group Delay Methods
%        for Identifying Glottal Closures in Voiced Speech," IEEE Trans on Speech & Audio Processing,
%        vol. 14, no. 2, pp. 456-466, Mar. 2006.
%   [3]  A. Kounoudes, P. A. Naylor, and M. Brookes, "The DYPSA algorithm for estimation of glottal
%        closure instants in voiced speech," in Proc ICASSP 2002, vol. 1, Orlando, 2002, pp. 349-352.
%   [4]  C. Ma, Y. Kamp, and L. F. Willems, "A Frobenius norm approach to glottal closure detection
%        from the speech signal," IEEE Trans. Speech Audio Processing, vol. 2, pp. 258-265, Apr. 1994.
%   [5]  A. Kounoudes, "Epoch Estimation for Closed-Phase Analysis of Speech," PhD Thesis,
%        Imperial College, 2001.

% Algorithm Parameters
%       The following parameters are defined in v_voicebox()
%
%   dy_cpfrac=0.3;           % presumed closed phase fraction of larynx cycle
%   dy_cproj=0.2;            % cost of projected candidate
%   dy_cspurt=-0.45;         % cost of a talkspurt
%   dy_dopsp=1;              % Use phase slope projection (1) or not (0)?
%   dy_ewdly=0.0008;         % window delay for energy cost function term [~ energy peak delay from closure] (sec)
%   dy_ewlen=0.003;          % window length for energy cost function term (sec)
%   dy_ewtaper=0.001;        % taper length for energy cost function window (sec)
%   dy_fwlen=0.00045;        % window length used to smooth group delay (sec)
%   dy_fxmax=500;            % max larynx frequency (Hz)
%   dy_fxmin=50;             % min larynx frequency (Hz)
%   dy_fxminf=60;            % min larynx frequency (Hz) [used for Frobenius norm only]
%   dy_gwlen=0.0030;         % group delay evaluation window length (sec)
%   dy_lpcdur=0.020;         % lpc analysis frame length (sec)
%   dy_lpcn=2;               % lpc additional poles
%   dy_lpcnf=0.001;          % lpc poles per Hz (1/Hz)
%   dy_lpcstep=0.010;        % lpc analysis step (sec)
%   dy_nbest=5;              % Number of NBest paths to keep
%   dy_preemph=50;           % pre-emphasis filter frequency (Hz) (to avoid preemphasis, make this very large)
%   dy_spitch=0.2;           % scale factor for pitch deviation cost
%   dy_wener=0.3;            % DP energy weighting
%   dy_wpitch=0.5;           % DP pitch weighting
%   dy_wslope=0.1;           % DP group delay slope weighting
%   dy_wxcorr=0.8;           % DP cross correlation weighting
%   dy_xwlen=0.01;           % cross-correlation length for waveform similarity (sec)

%   Revision History:
%
%   3.0 - 29 Jun 2006  - Rewrote DP function to improve speed
%   2.6 - 29 Jun 2006  - Tidied up algorithm parameters
%   2.4 - 10 Jun 2006  - Made into a single file aand put into VOICEBOX
%   2.3 - 18 Mar 2005  - Removed 4kHz filtering of phase-slope function
%   2.2 - 05 Oct 2004  -  dpgci uses the slopes returned from xewgrdel
%                      -  gdwav from speech with fs<9000 is not filtered
%                      -  Various outputs and inputs of functions have been
%                         removed since now there is no plotting
%   1.0 - 30 Jan 2001  - Initial version [5]

%   Bugs:
%         1. Allow the projections only to extend to the end of the larynx cycle
%         2. Compensate for false pitch period cost at the start of a voicespurt
%         3. Should include energy and phase-slope costs for the first closure of a voicespurt
%         4. should delete candidates that are too close to the beginning or end of speech for the cost measures
%            currently this is 200 samples fixed in the main routine but it should adapt to window lengths of
%            cross-correlation, lpc and energy measures.
%         5. Should have an integrated voiced/voiceless detector
%         6. Allow v_dypsa to be called in chunks for a long speech file
%         7. Do forward & backward passes to allow for gradient descent and/or discriminative training
%         8. Glottal opening approximation does not really belong in this function
%         9. The cross correlation window is asymmetric (and overcomplex) for no very good reason
%        10. Incorporate -0.5 factor into dy_wxcorr and abolish erroneous (nx2-1)/(nx2-2) factor
%        11. Add in talkspurt cost at the beginning rather than the end of a spurt (more efficient)
%        12. Remove qmin>2 condition from voicespurt start detection (DYPSA 2 compatibility) in two places
%        13. Include energy and phase-slope costs at the start of a voicespurt
%        14. Single-closure voicespurt are only allowed if nbest=1 (should always be forbidden)
%        15. Penultimate closure candidate is always acceptd
%        16. Final element of gcic, Cfn and Ch is unused
%        17. Needs to cope better with irregular voicing (e.g. creaky voice)
%        18. Should give non-integer GCI positions for improved accuracy
%        19. Remove constraint that first voicespurt cannot begin until qrmax after the first candidate

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

% Extract algorithm constants from VOICEBOX

dy_preemph=v_voicebox('dy_preemph');
dy_lpcstep=v_voicebox('dy_lpcstep');
dy_lpcdur=v_voicebox('dy_lpcdur');
dy_dopsp=v_voicebox('dy_dopsp');              % Use phase slope projection (1) or not (0)?
dy_ewtaper=v_voicebox('dy_ewtaper');        % Prediction order of FrobNorm method  in seconds
dy_ewlen=v_voicebox('dy_ewlen');        % windowlength of FrobNorm method  in seconds
dy_ewdly=v_voicebox('dy_ewdly');        % shift for assymetric speech shape at start of voiced cycle
dy_cpfrac=v_voicebox('dy_cpfrac');        % presumed ratio of larynx cycle that is closed
dy_lpcnf=v_voicebox('dy_lpcnf');          % lpc poles per Hz (1/Hz)
dy_lpcn=v_voicebox('dy_lpcn');            % lpc additional poles
dy_xwlen=v_voicebox('dy_xwlen');            % cross-correlation length for waveform similarity (sec)
dy_fxminf=v_voicebox('dy_fxminf');            % minimum pitch for Frobenius norm calculation

lpcord=ceil(fs*dy_lpcnf+dy_lpcn);       % lpc poles

%PreEmphasise input speech
s_used=filter([1 -exp(-2*pi*dy_preemph/fs)],1,s);

% perform LPC analysis, AC method with Hamming windowing
[ar, e, k] = v_lpcauto(s_used,lpcord,floor([dy_lpcstep dy_lpcdur]*fs));

if any(any(isinf(ar)))    % if the data is bad and gives infinite prediction coefficients we return with a warning
    warning('No GCIs returned');
    gci=[];
    return;
end;

% compute the prediction residual
r = v_lpcifilt(s_used,ar,k);

% compute the group delay function:  EW method from reference [2] above
[zcr_cand,sew,gdwav,toff]=xewgrdel(r,fs);
gdwav=-[zeros(toff,1); gdwav(1:end-toff)];
zcr_cand=[round(zcr_cand), ones(size(zcr_cand))];   %flag zero crossing candidates with ones

sew=0.5+sew';  %the phase slope cost of each candidate

pro_cand=[];
if dy_dopsp ~= 0
    pro_cand = psp(gdwav,fs);
    pro_cand = [pro_cand, zeros(length(pro_cand),1)]; %flag projected candidates with zeros
    sew =      [sew zeros(1,size(pro_cand,1))];      %the phase slope cost of a projected candidate is zero
end;

%Sort the zero crossing and projected candidates together and remove any candidates that
%are to close to the start or end of the speech signal because the cost functions
%need room either side.

[gcic,sin] = sortrows([zcr_cand; pro_cand],1);
sew=sew(sin);
saf=max([200,dy_xwlen*fs/2+1,fs/dy_fxminf]);
sin=find(and(saf<gcic,gcic<length(gdwav)-saf));
gcic=gcic(sin,:);
sew=sew(sin);

% compute the frobenious norm function used for a cost in the DP
fnwav=frobfun(s_used,dy_ewtaper*fs,dy_ewlen*fs,dy_ewdly*fs);

%Dynamic programming, picks the most likely candidates based on the pitch consistency, energy etc.
[gci] = dpgci(gcic, s_used(:), sew, fnwav, fs);

%Evaluate goi ... determined to be dy_cpfrac percentage of the larynx cylce away from the last gci
goi=simplegci2goi(gci,dy_cpfrac);

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

function z = psp(g,fs)
%PSP  Calculates the phase slope projections of the group delay function
%   Z = PSP(G) computes the

%   Revision History:
%       V1.0 July 12th 2002:
%            Nov. 6th 2002 : if statement added to remove "last" midpoint

g = g(:);

gdot = [diff(g);0];
gdotdot = [diff(gdot);0];

% find the turning points  as follows: [tp_number, index_of_tp, min(1) or max(-1), g(index_of_tp)]
turningPoints = zcr(gdot);
turningPoints = [[1:length(turningPoints)]', turningPoints, sign(gdotdot(turningPoints)), g(turningPoints)];

% useful for debug/plotting
%tplot = zeros(length(g),1);
%tplot(turningPoints(:,1)) = turningPoints(:,2);

% find any maxima which are < 0
negmaxima = turningPoints(find(turningPoints(:,3) == -1 & turningPoints(:,4) < 0 & turningPoints(:,1)~=1),:);  %Change 01.05.2003 JG: The first row can't be included

% find the midpoint between the preceding min and the negative max
nmi = negmaxima(:,1);
midPointIndex = turningPoints(nmi-1,2) + round(0.5*(turningPoints(nmi,2) - turningPoints(nmi-1,2)));
midPointValue = g(midPointIndex);

% project a zero crossing with unit slope
nz = midPointIndex - round(midPointValue);

% find any minima which are > 0
posminima = turningPoints(find(turningPoints(:,3) == 1 & turningPoints(:,4) > 0),:);

% find the midpoint between the positive min and the following max
pmi = posminima(:,1);

%Remove last midpoint if it is the last sample
if ~isempty(pmi), if pmi(end)==size(turningPoints,1), pmi=pmi(1:end-1); end; end;

midPointIndex = turningPoints(pmi,2) + round(0.5*(turningPoints(pmi+1,2) - turningPoints(pmi,2)));
midPointValue = g(midPointIndex);

% project a zero crossing with unit slope
pz = midPointIndex - round(midPointValue);

z = sort([nz;pz]);

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

function i = zcr(x, p)
%ZCR  Finds the indices in a vector to  zero crossings
%   I = ZCR(X) finds the indices of vector X which are closest to zero-crossings.
%   I = ZCR(X, P) finds indices for positive-going zeros-crossings for P=1 and
%   negative-going zero-crossings for P=0.

x = x(:);

if (nargin==2)
    if (p==0)
        z1 = zcrp(x);   % find positive going zero-crossings
    elseif (p==1)
        z1 = zcrp(-x);  % find negative going zero-crossings
    else
        error('ZCR: invalid input parameter 2: must be 0 or 1');
    end
else
    z1 = [zcrp(x); zcrp(-x)];
end

% find crossings when x==0 exactly
z0 = find( (x(1:length(x))==0) & ([x(2:length(x));0] ~= 0));

% concatenate and sort the two types of zero-crossings
i = sort([z0; z1]);

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

function zz = zcrp(xx)  %only used in zcr
% find positive-going zero-crossing
z1 = find(diff(sign(xx)) == -2);
% find which out of current sample or next sample is closer to zero
[m, z2] = min([abs(xx(z1)), abs(xx(z1+1))], [], 2);
zz =  z1 -1 + z2;

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

function [frob]=frobfun(sp,p,m,offset)

% [frob]=frobfun(sp,p,m)
%
% sp is the speech signal assumed to be preemphasised
% p  is the prediction order  : recomended to be 1 ms in above paper
% m  is the window length     : recomended to be 1 ms in above paper
% offset is shift for assymetric speech shape at start of voiced cycle -
% default 1.5ms.
%
% This function implements the frobenius norm based measure C defined in [4] below.
% It equals the square of the Frobenius norm of the m by p+1 data matrix divided by p+1
%
% Reference:
%   [4]  C. Ma, Y. Kamp, and L. F. Willems, "A Frobenius norm approach to glottal closure detection
%        from the speech signal," IEEE Trans. Speech Audio Processing, vol. 2, pp. 258"265, Apr. 1994.

%   Revision History:
%       V1.0 July 12th 2002:
%            Nov. 6th 2002 : if statement added to remove "last" midpoint

%force p m and offset to be integers
p=round(p);
m=round(m);
offset=round(offset);

w=(p+1)*ones(1,m+p);
w(1:p)=1:p;
w(m+1:p+m)=p:-1:1;

w=w./(p+1);
frob=filter(w,1,sp.^2);
frob(1:(round((p+m-1)/2) + offset))=[];

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

function goi=simplegci2goi(gci,pr)

% Estimate glottal opening instants by assuming a fixed closed-phase fraction

gci=round(gci);
maxpitch=max(medfilt1(diff(gci),7));

% calculate opening instants
for kg=1:length(gci)-1
    goi(kg)=gci(kg)+min(pr*(gci(kg+1)-gci(kg)),pr*maxpitch);
end;
kg=kg+1;
goi(kg)=round(gci(kg)+pr*(gci(kg)-gci(kg-1)));  %use the previous pitch period instead
goi=round(goi);

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

function [tew,sew,y,toff]=xewgrdel(u,fs)

% implement EW group delay epoch extraction

dy_gwlen=v_voicebox('dy_gwlen');          % group delay evaluation window length
dy_fwlen=v_voicebox('dy_fwlen');          % window length used to smooth group delay

% perform group delay calculation

gw=2*floor(dy_gwlen*fs/2)+1;            % force window length to be odd
ghw=window('hamming',gw,'s');
ghw = ghw(:);                           % force to be a column (dmb thinks window gives a row - and he should know as he wrote it!)
ghwn=ghw'.*(gw-1:-2:1-gw)/2;            % weighted window: zero in middle

u2=u.^2;
yn=filter(ghwn,1,u2);
yd=filter(ghw,1,u2);
yd(abs(yd)<eps)=10*eps;                 % prevent infinities
y=yn(gw:end)./yd(gw:end);               % delete filter startup transient
toff=(gw-1)/2;
fw=2*floor(dy_fwlen*fs/2)+1;            % force window length to be odd
if fw>1
    daw=window('hamming',fw,'s');
    y=filter(daw,1,y)/sum(daw);         % low pass filter
    toff=toff-(fw-1)/2;
end
[tew,sew]=v_zerocros(y,'n');              % find zero crossings

tew=tew+toff;                           % compensate for filter delay and frame advance

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

function Cfn=fnrg(gcic,frob,fs)

%Frobenious Energy Cost

dy_fxminf=v_voicebox('dy_fxminf');
frob=frob(:)';
mm=round(fs/dy_fxminf);
mfrob=v_maxfilt(frob,1,mm);
mfrob=[mfrob(floor(mm/2)+1:end) max(frob(end-ceil(mm/2):end))*ones(1,floor(mm/2))];
rfr=frob./mfrob;
Cfn=0.5-rfr(round(gcic));

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

function gci=dpgci(gcic, s, Ch, fnwav, fs)

%DPGCI   Choose the best Glottal Closure Instances with Dynamic Programming
%   gci=dpgci(gcic, s(:), Ch, fnwav, fs) returns vectors of sample indices corresponding
%   to the instants of glottal closure in the speech signal s at sampling frequency fs Hz.
%
%   Inputs:
%   gcic    is a matrix whos first column are the glottal closure instance candidates and
%           the second column is 1 if the corresponding gci is derived from a zero crossing
%           but zero if the gci is from a a projected zero crossing
%   s       is the speech signal - MUST be a column vector
%   Ch      the phase slope cost of every candidate
%   fnwav   is the frobenious norm function of s
%   fs      is the sampling frequncy
%
%   Outputs:
%   gci     is a vector of glottal closure instances chosen by the DP

%   Revision History:
%   Bugs:  Constants are hardwired but defined in a structure like pv (defined in grpdelpv)
%

% get algorithm parameters from v_voicebox()

dy_fxmin=v_voicebox('dy_fxmin');        % min larynx frequency (Hz)
dy_fxmax=v_voicebox('dy_fxmax');        % min larynx frequency (Hz)
dy_xwlen=v_voicebox('dy_xwlen');        % cross-correlation length for waveform similarity (sec)
dy_nbest=v_voicebox('dy_nbest');        % Number of NBest paths to keep
dy_spitch=v_voicebox('dy_spitch');              % scale factor for pitch deviation cost
wproj=v_voicebox('dy_cproj');           % cost of projected candidate
dy_cspurt=v_voicebox('dy_cspurt');           % cost of a talkspurt
dy_wpitch=v_voicebox('dy_wpitch');           % DP pitch weighting
dy_wener=v_voicebox('dy_wener');           % DP energy weighting
dy_wslope=v_voicebox('dy_wslope');           % DP group delay slope weighting
dy_wxcorr=v_voicebox('dy_wxcorr');           % DP cross correlation weighting

%Constants
Ncand=length(gcic);
sv2i=-(2*dy_spitch^2)^(-1);                 % scale factor for pitch deviation cost
nxc=ceil(dy_xwlen*fs);                      % cross correlation window length in samples
% === should delete any gci's that are within this of the end.
% === and for energy window

%Limit the search:
qrmin=ceil(fs/dy_fxmax);
qrmax=floor(fs/dy_fxmin);

%Cost and tracking r = current, q = previous, p = preprevious
cost=zeros(Ncand, dy_nbest); cost(:,:)=inf;    %Cost matrix, one row for each candidate
maxcost=zeros(Ncand,1); maxcost(:,:)=inf;   %Maximum cost in each row
imaxcost=ones(Ncand,1);                     %Index of maximum cost

prev = ones(Ncand, dy_nbest);                  %index of previous, q candidates
ind = ones(Ncand, dy_nbest);                   %index of p in row q (from prev)
qbest = [zeros(Ncand,1), ones(Ncand,2)]; % the minimum cost in any previous q [cost,q,i]

Cfn=fnrg(gcic(:,1),fnwav,fs);  %Frob.Energy Cost

%Add start and end state
% === should probably delete candidates that are too close to either end of the input
% === why do we ever need the additional one at the tail end ?
gcic=[[gcic(1,1)-qrmax-2 0];gcic;[gcic(end,1)+qrmax+2 0]];
Cfn=[0 Cfn 0];
Ch = [0 Ch 0];

% first do parallelized version

% rather complicated window specification is for compatibility with DYPSA 2
% === +1 below is for compatibility - probably a bug
wavix=(-floor(nxc/2):floor(nxc/2)+1)';                 % indexes for segments [nx2,1]
nx2=length(wavix);
sqnx2=sqrt(nx2);
g_cr=dy_wener*Cfn+dy_wslope*Ch+wproj*(1-gcic(:,2))';  % fixed costs

g_n=gcic(:,1)';                  % gci sample number [1,Ncand+2]
g_pr=gcic(:,2)';                 % unprojected flag [1,Ncand+2]
g_sqm=zeros(1,Ncand+1);         % stores: sqrt(nx2) * mean for speech similarity waveform
g_sd=zeros(1,Ncand+1);         % stores: 1/(Std deviation * sqrt(nx2)) for speech similarity waveform
f_pq=zeros((Ncand+1)*dy_nbest,1);   % (q-p) period for each node
f_c=repmat(Inf,(Ncand+1)*dy_nbest,1);    % cumulative cost for each node - initialise to inf
f_c(1)=0;                       % initial cost of zero for starting node
% f_costs=zeros(Ncand*dy_nbest,6);   % === debugging only remember costs of candidate
f_f=ones((Ncand+1)*dy_nbest,1);    % previous node in path
f_fb=ones((Ncand+1),1);    % points back to best end-of-spurt node
fbestc=0;                       % cost of best end-of-spurt node

qmin=2;
for r=2:Ncand+1
%     if r==86
%         r;
%     end
    r_n=g_n(r);             % sample number of r = current candidate
    rix=dy_nbest*(r-1)+(1:dy_nbest);    % index range within node variables

    % determine the range of feasible q candidates
    qmin0=qmin;
    qmin=find(g_n(qmin0-1:r-1)<r_n-qrmax);      % qmin is the nearest candidate that is >qrmax away
    qmin=qmin(end)+qmin0-1;             % convert to absolute index of first viable candidate
    qmax=find(g_n(qmin-1:r-1)<=r_n-qrmin);      % qmax is the nearest candidate that is >=qrmin away
    qmax=qmax(end)+qmin-2;

    % calculate waveform similarity cost measure statistics

    sr=s(r_n+wavix);        % note s MUST be a column vector so sr is also
    wsum=sum(sr);
    g_sqm(r)=wsum/sqnx2;                % mean * sqrt(nx2)
    g_sd(r)=1/sqrt(sr.'*sr-wsum^2/nx2);   % 1/(Std deviation * sqrt(nx2))

    % now process the candidates

    if qmin<=qmax
        qix=qmin:qmax;      % q index
        nq=length(qix);
        % === should integrate the -0.5 into dy_wxcorr
        % === the factor (nx2-1)/(nx2-2) is to compensate for a bug in swsc()
        q_cas=-0.5*(nx2-1)/(nx2-2)*dy_wxcorr*(sum(s(repmat(g_n(qix),nx2,1)+repmat(wavix,1,nq)).*repmat(sr,1,nq),1)-g_sqm(qix)*g_sqm(r)).*g_sd(qix)*g_sd(r);
        % compare: i=35; Ca=swsc(g_n(qix(i)),g_n(r),s,fs); [i qix(i) r  g_n(qix(i)) g_n(r) dy_wxcorr*Ca q_cas(i)]

        % now calculate pitch deviation cost

        fix = 1+(qmin-1)*dy_nbest:qmax*dy_nbest;    % node index range
        f_qr=repmat(r_n-g_n(qix),dy_nbest,1);    % (r-p) period for each node
        f_pr=f_qr(:)+f_pq(fix);
        % === could absorb the 2 into sv2i
        f_nx=2-2*f_pr./(f_pr+abs(f_qr(:)-f_pq(fix)));
        f_cp=dy_wpitch*(0.5-exp(sv2i*f_nx.^2));
        % === fudge to match dypsa2.4 - could more efficiently be added
        % === onto the cost of a talkspurt end
        % === should be a v_voicebox parameter anyway
        f_cp(f_pq(fix)==0)=dy_cspurt*dy_wpitch;

        % now find the N-best paths

        [r_cnb,nbix]=sort(f_c(fix)+f_cp+reshape(repmat(q_cas,dy_nbest,1),nq*dy_nbest,1));
        f_c(rix)=r_cnb(1:dy_nbest)+g_cr(r);     % costs
        f_f(rix)=nbix(1:dy_nbest)+(qmin-1)*dy_nbest;       % traceback nodes
        f_pq(rix)=f_qr(nbix(1:dy_nbest));       % previous period
        % === f_costs is only for debugging
%         r;
%         f_costs(rix,1)=f_c(fix(nbix(1:dy_nbest)));
%         f_costs(rix,2)=wproj*(1-gcic(r,2));
%         f_costs(rix,3)=f_cp(nbix(1:dy_nbest));
%         f_costs(rix,4)=dy_wener*Cfn(r);
%         f_costs(rix,5)=dy_wslope*Ch(r);
%         f_costs(rix,6)=reshape(q_cas(1+floor((nbix(1:dy_nbest)-1)/dy_nbest)),dy_nbest,1);

        % check cost of using this candidate as the start of a new spurt
        % ==== the qmin>2 condition is for compatibility with v_dypsa 2 and
        % prevents any spurts starting until at least qrmax past the first
        % gci. This is probably a bug (see again below)
        iNb=rix(end);
        if (qmin>2) && (f_c(f_fb(qmin-1))+wproj*(1-gcic(r,2))<f_c(iNb))        % compare with worst of Nbest paths
            f_f(iNb)=f_fb(qmin-1);
            % === for now we exclude the energy and phase-slope costs for compatibility with dypsa2
            % === this is probably a bug
            f_c(iNb)=f_c(f_fb(qmin-1))+wproj*(1-gcic(r,2));     % replace worst of the costs with start voicespurt cost
            f_pq(iNb)=0;                    % false pq period
        end
        if f_c(rix(1))<fbestc
            f_fb(r)=rix(1);                          % points to the node with lowest end-of-spurt cost
            % === should compensate for the pitch period cost incurred at the start of the next spurt
            % === note that a node can never be a one-node voicespurt on its own unless dy_nbest=1
            % since the start voices[purt option replaced the worst Nbest cost. This is probably good but
            % is a bit inconsistent.
            fbestc=f_c(rix(1));
        else
            f_fb(r)=f_fb(r-1);
        end
    else            % no viable candidates - must be the start of a voicespurt if anything
        % === for now we exclude the energy and phase-slope costs for compatibility with dypsa2
        % === this is probably a bug
        % ==== the qmin>2 condition is for compatibility with v_dypsa 2 and
        % prevents any spurts starting until at least qrmax past the first
        % gci. This is probably a bug (see again above)
        if (qmin>2)
            f_c(rix(1))=f_c(f_fb(qmin-1))+wproj*(1-gcic(r,2));  % cost of new voicespurt
            f_f(rix)=f_fb(qmin-1);                              % traceback to previous talkspurt end
            f_pq(rix)=0;                                        % previous period
        end
        f_fb(r)=f_fb(r-1);                                  % cannot be the end of a voicespurt
    end
end

% now do the traceback

gci = zeros(1,Ncand+1);

% === for compatibility with dypsa2, we force the penultimate candidate to be accepted
% === should be: i=f_fb(Ncand+1) but instead we pick the best of the penultimate candidate
i=rix(1)-dy_nbest;
if f_c(i-dy_nbest+1)<f_c(i)     % check if start of a talkspurt
    i=i-dy_nbest+1;
end
k=1;
while i>1
    j=1+floor((i-1)/dy_nbest);          % convert node number to candidate number
    gci(k)=g_n(j);
    i=f_f(i);
    k=k+1;
end
gci=gci(k-1:-1:1);           % put into ascending order

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

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function bessel_y1_spherical_values_test ( )

%*****************************************************************************80
%
%% BESSEL_Y1_SPHERICAL_VALUES_TEST tests BESSEL_Y1_SPHERICAL_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BESSEL_Y1_SPHERICAL_VALUES_TEST:\n' );
  fprintf ( 1, '  BESSEL_Y1_SPHERICAL_VALUES stores values of\n' );
  fprintf ( 1, '  the y1 spherical Bessel function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X            FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = bessel_y1_spherical_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ y, m, d, f ] = jed_to_ymdf_khwarizmian ( jed )

%*****************************************************************************80
%
%% JED_TO_YMDF_KHWARIZMIAN converts a JED to a Khwarizmian YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Algorithm F,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, pages 324-325.
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, integer Y, M, D, real F,
%    the YMDF date.
%

%
%  Determine the computational date (Y'/M'/D').
%
  j = floor ( jed + 0.5 );
  f = ( jed + 0.5 ) - j;

  j_prime = j + 317;

  y_prime = floor ( j_prime / 365 );
  t_prime = mod ( j_prime, 365 );
  m_prime = floor ( t_prime / 30 );
  d_prime = mod ( t_prime, 30 );
%
%  Convert the computational date to a calendar date.
%
  d = d_prime + 1;
  m = mod ( m_prime, 13 ) + 1;
  y = y_prime - 5348 + floor ( ( 13 - m ) / 13 );

  return
end


text:
%GAUSS Generation of a multivariate Gaussian dataset
%
% 	A = GAUSS(N,U,G,LABTYPE)
%
% INPUT (in case of generation a 1-class dataset in K dimensions)
%   N		    Number of objects to be generated (default 50).
%   U		    Desired mean (vector of length K).
%   G       K x K covariance matrix. Default eye(K).
%   LABTYPE Label type (default 'crisp')
%
% INPUT (in case of generation a C-class dataset in K dimensions)
%   N       Vector of length C with numbers of objects per class.
%   U       C x K matrix with class means, or
%           Dataset with means, labels and priors of classes
%           (default: zeros(C,K))
%   G       K x K x C covariance matrix of right size.
%           Default eye(K);
%   LABTYPE	Label type (default 'crisp')
%
% OUTPUT
%   A       Dataset containing multivariate Gaussian data
%
% DESCRIPTION
% Generation of N K-dimensional Gaussian distributed samples for C classes.
% The covariance matrices should be specified in G (size K*K*C) and the
% means, labels and prior probabilities can be defined by the dataset U with
% size (C*K). If U is not a dataset, it should be a C*K matrix and A will
% be a dataset with C classes.
%
% If N is a vector, exactly N(I) objects are generated for class I, I = 1..C.
%
% EXAMPLES
% 1. Generation of 100 points in 2D with mean [1 1] and default covariance
%    matrix:
%
%        GAUSS(100,[1 1])
%
% 2. Generation of 50 points for each of two 1-dimensional distributions with
%    mean -1 and 1 and with variances 1 and 2:
%
%	     GAUSS([50 50],[-1;1],CAT(3,1,2))
%
%   Note that the two 1-dimensional class means should be given as a column
%   vector [1;-1], as [1 -1] defines a single 2-dimensional mean. Note that
%   the 1-dimensional covariance matrices degenerate to scalar variances,
%   but have still to be combined into a collection of square matrices using
%   the CAT(3,....) function.
%
% 3. Generation of 300 points for 3 classes with means [0 0], [0 1] and
%    [1 1] and covariance matrices [2 1; 1 4], EYE(2) and EYE(2):
%
%      GAUSS(300,[0 0; 0 1; 1 1]*3,CAT(3,[2 1; 1 4],EYE(2),EYE(2)))
%
% SEE ALSO (PRTools Guide)
% DATASETS, PRDATASETS

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function a = gauss(n,u,g,labtype)

		if (nargin < 1)
		prwarning (2,'number of samples not specified, assuming N = 50');
		n = 50;
		end
		cn = length(n);
	if (nargin < 2)
		prwarning (2,'means not specified; assuming one dimension, mean zero');
		u = zeros(cn,1);
	end;
	if (nargin < 3)
		prwarning (2,'covariances not specified, assuming unity');
	 	g = eye(size(u,2));
	end
	if (nargin < 4)
		prwarning (3,'label type not specified, assuming crisp');
		labtype = 'crisp';
	end

	% Return an empty dataset if the number of samples requested is 0.

	if (length(n) == 1) & (n == 0)
		a = prdataset([]);
		return
	end

	% Find C, desired number of classes based on U and K, the number of
	% dimensions. Make sure U is a dataset containing the means.

	if (isa(u,'prdataset'))
		[m,k,c] = getsize(u);
		lablist = getlablist(u);
		p = getprior(u);
		if c == 0
			u = double(u);
		end
	end
	if isa(u,'double')
		[m,k] = size(u);
		c = m;
		lablist = genlab(ones(c,1));
		u = prdataset(u,lablist);
		p = ones(1,c)/c;
	end

	if (cn ~= c) & (cn ~= 1)
		error('The number of classes specified by N and U does not match');
	end

 	% Generate a class frequency distribution according to the desired priors.
	n = genclass(n,p);

	% Find CG, the number of classes according to G.
	% Make sure G is not a dataset.

	if (isempty(g))
		g = eye(k);
		cg = 1;
	else
		g = real(+g); [k1,k2,cg] = size(g);
		if (k1 ~= k) | (k2 ~= k)
			error('The number of dimensions of the means U and covariance matrices G do not match');
		end
		if (cg ~= m & cg ~= 1)
			error('The number of classes specified by the means U and covariance matrices G do not match');
		end
	end

	% Create the data A by rotating and scaling standard normal distributions
	% using the eigenvectors of the specified covariance matrices, and adding
	% the means.

	%a = [];
  a = zeros(sum(n),k);
  nn = 0;
	for i = 1:m
		j = min(i,cg);						% Just in case CG = 1 (if G was not specified).

		% Sanity check: user can pass non-positive definite G.
		[V,D] = preig(g(:,:,j)); V = real(V); D = real(D); D = max(D,0);
    a(nn+1:nn+n(i),:) = randn(n(i),k)*sqrt(D)*V' + repmat(+u(i,:),n(i),1);
		%a = [a; randn(n(i),k)*sqrt(D)*V' + repmat(+u(i,:),n(i),1)];
    nn = nn+n(i);
	end

	% Convert A to dataset by adding labels and priors.

	labels = genlab(n,lablist);
	a = prdataset(a,labels,'lablist',lablist,'prior',p);

	% If non-crisp labels are requested, use output of Bayes classifier.
	switch (labtype)
		case 'crisp'
			;
		case 'soft'
			w = nbayesc(u,g);
			targets = a*w*classc;
			a = setlabtype(a,'soft',targets);
		otherwise
			error(['Label type ' labtype ' not supported'])
	end

	a = setname(a,'Gaussian Data');

return

author: marianux

repo: ecg-kit

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


text:
function deleteOptional(fname)
    if(exist(fname,'file'))
        delete(fname);
    elseif(exist(fname,'dir'))
        rmdir(fname);
    else
        return;
    end
end

repo: CategoryShapes


text:
classdef TestResize
    %TestResize

    methods (Static)
        function test_1
            img = imread(fullfile(mexopencv.root(),'test','img001.jpg'));
            ref = img(1:2:end,1:2:end,:);
            dst = cv.resize(img, 0.5, 0.5, 'Interpolation','Nearest');
            validateattributes(dst, {class(img)}, {'size',size(ref)});
            assert(isequal(ref, dst));
        end

        function test_2
            img = imread(fullfile(mexopencv.root(),'test','img001.jpg'));
            ref = img(1:2:end,1:2:end,:);
            dst = cv.resize(img, [256,256], 'Interpolation','Nearest');
            validateattributes(dst, {class(img)}, {'size',size(ref)});
            assert(all(abs(ref(:) - dst(:)) < 1e-5));
        end

        function test_3
            img = imread(fullfile(mexopencv.root(),'test','img001.jpg'));
            interps = {'Nearest', 'Linear', 'Cubic', 'Area', 'Lanczos4'};
            for i=1:numel(interps)
                dst = cv.resize(img, [700 800], 'Interpolation',interps{i});
                validateattributes(dst, {class(img)}, {'size',[800 700 size(img,3)]});
            end
        end

        function test_error_argnum
            try
                cv.resize();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function root_num = p16_root_num ( )

%*****************************************************************************80
%
%% P16_ROOT_NUM returns the number of known roots for problem 16.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer ROOT_NUM, the number of known roots.
%
  root_num = 0;

  return
end


text:
function [y,deriv] = square_distance_mv2df(w,input_data,new_dim)
% This is an MV2DF. See MV2DF_API_DEFINITION.readme.
%
% The function computes the square distance of the vectors for each trial.
%   y.' = sum((W(:,1:end-1).'*input_data + W(:,end)).^2,1)
%
%      W is the augmented matrix [M c] where M maps a score vector
%      to a lower dimensional space and c is an offset vector in
%      the lower dimensional space.
%
% Parameters:
%     w: the vectorized version of the W matrix
%     input_data: is an M-by-T matrix of input vectors of length M, for each of T
%             trials.
%     new_dim: the dimension of vectors in the lower dimensional space.
%

if nargin==0
    test_this();
    return;
end

if isempty(w)
    [dim, num_trials] = size(input_data);
    map = @(w) map_this(w,input_data,dim,new_dim);
    transmap = @(w) transmap_this(w,input_data,num_trials,new_dim);
    delta = linTrans(w,map,transmap);
    y = sums_of_squares(delta,new_dim);
    return;
end

if isa(w,'function_handle')
    f = square_distance_mv2df([],input_data,new_dim);
    y = compose_mv(f,w,[]);
    return;
end

f = square_distance_mv2df([],input_data,new_dim);
if nargout==1
    y = f(w);
else
    [y,deriv] = f(w);
end

function y = map_this(w,input_data,dim,new_dim)
% V = [input_data; ones(1,num_trials)];
W = reshape(w,new_dim,dim+1);
y = bsxfun(@minus,W(:,1:end-1)*input_data,W(:,end));
y = y(:);

function dx = transmap_this(dy,input_data,num_trials,new_dim)
dY = reshape(dy,new_dim,num_trials);
% V = [input_data; ones(1,num_trials)];
% Vt = V.';
% dX = dY*Vt;
dYt = dY.';
dYtSum = sum(dYt,1);
dX = [input_data*dYt;-dYtSum].';
dx = dX(:);

function test_this()
K = 5;
N = 10;
P = 3;
M = randn(P,N);
c = randn(P,1);
W = [M c];
w = W(:);
input_data = randn(N,K);

f = square_distance_mv2df([],input_data,P);
test_MV2DF(f,w);

author: nesl


text:
function res = imSkeleton(img, varargin)
% Homothopic skeleton of a binary image.
%
%   SK = imSkeleton(IMG)
%   Computes the skeleton of the binary image IMG. The skeleton is computed
%   by iterative homothopic thinning of the input image.
%   Basically, this is just a wrapper for the "bwmorph" function. If the
%   "bwskel" function (introduced in 2018) is found, it is used instead.
%
%   Example
%     img = imread('circles.png');
%     skel = imSkeleton(img);
%     ovr = imOverlay(img, skel);
%     imshow(ovr)
%
%   See also
%     bwmorph, bwskel, imLabelSkeleton, imBoundary, imOverlay
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2012-05-16,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2012 INRA - Cepia Software Platform.

if exist('bwskel', 'file')
    res = bwskel(img);
else
    res = bwmorph(img, 'thin', Inf);
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
%This function develop by Renoald Tang
%University Teknologi Malaysia, Photogrammetry and Laser scanning group
%for academic purpose
%Email:renoald@live.com
%This function read Ply file

function readPly(file)
nv=numvert(file);
fprintf('Number_of_point:%d\n',nv);
nh=numhead(file);
fprintf('Number of end line:%d\n',nh);
nf=numface(file,nh);
fprintf('Number_of_face:%d\n',nf);
value(nh,file,nv)

%return number of vertext
function nv=numvert(file)
C=0;
fid=fopen(file);
while 1
    tline = fgetl(fid);
    C=C+1;
    if C==4

      num=textscan(tline,'%s %s %d');
        break;
    end
    if ~ischar(tline),   break,   end

end
fclose(fid);
nv=num{1,3};
%return number of face
function nf=numface(file,nh)
nh=nh-2;
C=0;
fid=fopen(file);
while 1
    tline = fgetl(fid);
    C=C+1;
     if ~ischar(tline),   break,   end

    if C==nh

      num=textscan(tline,'%s %s %d');
        break;
    end

end
fclose(fid);
nf=num{1,3};
%return the number of end_hearder
function nh=numhead(file)
C=0;
fid=fopen(file);
while 1
    tline = fgetl(fid);
    C=C+1;
    map1=findstr(tline,'end_header');
  b1=isempty(map1);
  if b1==0
      nh=C;

      break;
  end

    if ~ischar(tline),   break,   end

end
fclose(fid);
function value(nh,file,n)
%[K1,R]=check(file);
C=0;
K=0;
fid=fopen(file);
fid1=fopen('point','w');
fid2=fopen('face','w');
while 1
    tline=fgetl(fid);
    C=C+1;
     if ~ischar(tline) , break , end
    if C > nh
       K=K+1;
       if K <=n
           fprintf(fid1,'%s\n',tline);
       else
           fprintf(fid2,'%s\n',tline);
       end

    end

end
fclose(fid);
fclose(fid1);
fclose(fid2);%return and extract face and point value
%disp(K)

author: Sable

repo: mcbench-benchmarks


text:
function c = tapas_hgf_ar1_binary_mab_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the Hierarchical Gaussian Filter (HGF) for AR(1) processes in a
% multi-armded bandit situation for binary inputs in the absence of perceptual uncertainty.
%
% The HGF is the model introduced in
%
% Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation
% for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39.
%
% The binary HGF model has since been augmented with a positive factor kappa1 which
% scales the second level with respect to the first, i.e., the relation between the
% first and second level is
%
% p(x1=1|x2) = s(kappa1*x2), where s(.) is the logistic sigmoid.
%
% By default, kappa1 is fixed to 1, leading exactly to the model introduced in
% Mathys et al. (2011).
%
% This file refers to BINARY inputs (Eqs 1-3 in Mathys et al., (2011));
% for continuous inputs, refer to tapas_hgf_config.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The HGF configuration consists of the priors of parameters and initial values. All priors are
% Gaussian in the space where the quantity they refer to is estimated. They are specified by their
% sufficient statistics: mean and variance (NOT standard deviation).
%
% Quantities are estimated in their native space if they are unbounded (e.g., the omegas). They are
% estimated in log-space if they have a natural lower bound at zero (e.g., the sigmas).
%
% The phis are estimated in 'logit space' because they are confined to the interval from 0 to 1.
% 'Logit-space' is a logistic sigmoid transformation of native space with a variable upper bound
% a>0:
%
% tapas_logit(x) = ln(x/(a-x)); x = a/(1+exp(-tapas_logit(x)))
%
% Parameters can be fixed (i.e., set to a fixed value) by setting the variance of their prior to
% zero. Aside from being useful for model comparison, the need for this arises whenever the scale
% and origin at the j-th level are arbitrary. This is the case if the observation model does not
% contain the representations mu_j and sigma_j. A choice of scale and origin is then implied by
% fixing the initial value mu_j_0 of mu_j and either kappa_j-1 or omega_j-1.
%
% Fitted trajectories can be plotted by using the command
%
% >> tapas_hgf_binary_mab_plotTraj(est)
%
% where est is the stucture returned by tapas_fitModel. This structure contains the estimated
% perceptual parameters in est.p_prc and the estimated trajectories of the agent's
% representations (cf. Mathys et al., 2011). Their meanings are:
%
%         est.p_prc.mu_0       row vector of initial values of mu (in ascending order of levels)
%         est.p_prc.sa_0       row vector of initial values of sigma (in ascending order of levels)
%         est.p_prc.phi        row vector of phis (representing reversion slope to attractor; in ascending order of levels)
%         est.p_prc.m          row vector of ms (representing attractors; in ascending order of levels)
%         est.p_prc.ka         row vector of kappas (in ascending order of levels)
%         est.p_prc.om         row vector of omegas (in ascending order of levels)
%
% Note that the first entry in all of the row vectors will be NaN because, at the first level,
% these parameters are either determined by the second level (mu_0 and sa_0) or undefined (rho,
% kappa, and omega).
%
%         est.traj.mu          mu (rows: trials, columns: levels, 3rd dim: bandits)
%         est.traj.sa          sigma (rows: trials, columns: levels, 3rd dim: bandits)
%         est.traj.muhat       prediction of mu (rows: trials, columns: levels, 3rd dim: bandits)
%         est.traj.sahat       precisions of predictions (rows: trials, columns: levels, 3rd dim: bandits)
%         est.traj.v           inferred variance of random walk (rows: trials, columns: levels)
%         est.traj.w           weighting factors (rows: trials, columns: levels)
%         est.traj.da          volatility prediction errors  (rows: trials, columns: levels)
%         est.traj.ud          updates with respect to prediction  (rows: trials, columns: levels)
%         est.traj.psi         precision weights on prediction errors  (rows: trials, columns: levels)
%         est.traj.epsi        precision-weighted prediction errors  (rows: trials, columns: levels)
%         est.traj.wt          full weights on prediction errors (at the first level,
%                                  this is the learning rate) (rows: trials, columns: levels)
%
% Note that in the absence of sensory uncertainty (which is the assumption here), the first
% column of mu, corresponding to the first level, will be equal to the inputs. Likewise, the
% first column of sa will be 0 always.
%
% Tips:
% - When analyzing a new dataset, take your inputs u and responses y and use
%
%   >> est = tapas_fitModel(y, u, 'tapas_hgf_ar1_binary_mab_config', 'tapas_bayes_optimal_binary_config');
%
%   to determine the Bayes optimal perceptual parameters (given your current priors as defined in
%   this file here, so choose them wide and loose to let the inputs influence the result). You can
%   then use the optimal parameters as your new prior means for the perceptual parameters.
%
% - If you get an error saying that the prior means are in a region where model assumptions are
%   violated, lower the prior means of the omegas, starting with the highest level and proceeding
%   downwards.
%
% - Alternatives are lowering the prior means of the kappas, if they are not fixed, or adjusting
%   the values of the kappas or omegas, if any of them are fixed.
%
% - If the log-model evidence cannot be calculated because the Hessian poses problems, look at
%   est.optim.H and fix the parameters that lead to NaNs.
%
% - Your guide to all these adjustments is the log-model evidence (LME). Whenever the LME increases
%   by at least 3 across datasets, the adjustment was a good idea and can be justified by just this:
%   the LME increased, so you had a better model.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013-2017 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'hgf_ar1_binary_mab';

% Number of levels (minimum: 3)
c.n_levels = 3;

% Number of bandits
c.n_bandits = 3;

% Coupling
% This may only be set to true if c.n_bandits is set to 2 above. If
% true, it means that the two bandits' winning probabilities are
% coupled in the sense that they add to 1 and are both updated on
% each trial even though only the outcome for one of them is observed.
c.coupled = false;

% Input intervals
% If input intervals are irregular, the last column of the input
% matrix u has to contain the interval between inputs k-1 and k
% in the k-th row, and this flag has to be set to true
c.irregular_intervals = false;

% Sufficient statistics of Gaussian parameter priors

% Initial mus and sigmas
% Format: row vectors of length n_levels
% For all but the first two levels, this is usually best
% kept fixed to 1 (determines origin on x_i-scale). The
% first level is NaN because it is determined by the second,
% and the second implies neutrality between outcomes when it
% is centered at 0.
c.mu_0mu = [NaN, 0, 1];
c.mu_0sa = [NaN, 1, 1];

c.logsa_0mu = [NaN,   log(0.1), log(1)];
c.logsa_0sa = [NaN,          1,      1];

% Phis
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
% Fix this to zero (-Inf in logit space) to set to zero.
c.logitphimu = [NaN, tapas_logit(0.4,1), tapas_logit(0.2,1)];
c.logitphisa = [NaN,                  1,                   1];

% ms
% Format: row vector of length n_levels.
% This should be fixed for all levels where the omega of
% the next lowest level is not fixed because that offers
% an alternative parametrization of the same model.
c.mmu = [NaN, c.mu_0mu(2), c.mu_0mu(3)];
c.msa = [NaN,           0,           0];

% Kappas
% Format: row vector of length n_levels-1.
% Fixing log(kappa1) to log(1) leads to the original HGF model.
% Higher log(kappas) should be fixed (preferably to log(1)) if the
% observation model does not use mu_i+1 (kappa then determines the
% scaling of x_i+1).
c.logkamu = [log(1), log(1)];
c.logkasa = [     0,    0.1];

% Omegas
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
c.ommu = [NaN,   -2,   -2];
c.omsa = [NaN,    1,    1];

% Gather prior settings in vectors
c.priormus = [
    c.mu_0mu,...
    c.logsa_0mu,...
    c.logitphimu,...
    c.mmu,...
    c.logkamu,...
    c.ommu,...
         ];

c.priorsas = [
    c.mu_0sa,...
    c.logsa_0sa,...
    c.logitphisa,...
    c.msa,...
    c.logkasa,...
    c.omsa,...
         ];

% Check whether we have the right number of priors
expectedLength = 5*c.n_levels+(c.n_levels-1);
if length([c.priormus, c.priorsas]) ~= 2*expectedLength;
    error('tapas:hgf:PriorDefNotMatchingLevels', 'Prior definition does not match number of levels.')
end

% Model function handle
c.prc_fun = @tapas_hgf_ar1_binary_mab;

% Handle to function that transforms perceptual parameters to their native space
% from the space they are estimated in
c.transp_prc_fun = @tapas_hgf_ar1_binary_mab_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
domain = [-1,1,-1,1,-1,1];
bc = [1,1,1,1,1,1]; % Dirichelet BC

load('BoxTorusMesh.mat')
% epsm = 8.85*10^(-5)*5;
% epsp = 8.85*10^(-5);
% sigm = 100;
% sigp = 1;
% mum = (4*pi);
% mup = (4*pi)*2;
% x0=0; y0=0; z0=-0.3; r1=0.2; r2=pi/5; omega = 20; stre = 1;
% pde = TorusTime1(mum,mup,sigm,sigp,epsm,epsp,omega,stre,x0,y0,z0,r1,r2);
% TEND = 0.1;
epsm = 8.854*10^(-2);
epsp = 8.854*10^(-3);
sigm = 10;%100;
sigp = 1;%1;
mum = 4*pi*5;
mup = 4*pi;
x0=0; y0=0; z0=-0.3; r1=0.2; r2=pi/5; omega = 1; a = omega*sqrt(epsp*mup); b= 120; intPt = -1;
pde = TorusTimeInitial3(mum,mup,sigm,sigp,epsm,epsp,omega,x0,y0,z0,r1,r2,a,b,intPt);
TEND = 0.8;
InitCond = 1;
BDCond = 1;

nx = 64;  h=(domain(2) - domain(1))/nx;
ny = nx;
nz = nx;

Ntime = ceil(TEND/sqrt(8.854*10^(-3)*4*pi*2))*5*nx;%ceil(sqrt(nx));
deltaT = TEND/Ntime;
mesh = enrichMesh3D(mesh,0);
fem = genNedFEM3D(mesh,bc);
bm = epsm/deltaT^2 + sigm/(2*deltaT);
bp = epsp/deltaT^2 + sigp/(2*deltaT);
am = mum^(-1); ap = mup^(-1);

S = globMatrixNedFit3D(am,ap,1,mesh,fem,fem);
Me = globMatrixNedFit3D(epsm,epsp,0,mesh,fem,fem);
Ms = globMatrixNedFit3D(sigm,sigp,0,mesh,fem,fem);
Atotal = Me/deltaT^2 + Ms/(2*deltaT) + S/2;

Ndof = size(mesh.e,1);
bdidx = zeros(Ndof,1);
isBdEdge = true(Ndof,1);
isBdEdge(fem.mapper) = false;
bdidx(isBdEdge) = 1;
Tbd = spdiags(bdidx,0,Ndof,Ndof);
T = spdiags(1-bdidx,0,Ndof,Ndof);
A = T*Atotal*T + Tbd;

if InitCond == 0
    uppre = zeros(size(mesh.e,1),1);
    upre = zeros(size(mesh.e,1),1);
elseif InitCond ~=0
    tu1 = sum(feval(pde.E1,fem.gex,fem.gey,fem.gez,0).*fem.gew,2);
    tu2 = sum(feval(pde.E2,fem.gex,fem.gey,fem.gez,0).*fem.gew,2);
    tu3 = sum(feval(pde.E3,fem.gex,fem.gey,fem.gez,0).*fem.gew,2);
    tgt = mesh.p(mesh.e(:,2),:) - mesh.p(mesh.e(:,1),:);
    tgt = tgt./sum(tgt.^2,2).^(1/2);
    uppre = tu1.*tgt(:,1) + tu2.*tgt(:,2) + tu3.*tgt(:,3);
    tu1 = sum(feval(pde.Et1,fem.gex,fem.gey,fem.gez,0).*fem.gew,2);
    tu2 = sum(feval(pde.Et2,fem.gex,fem.gey,fem.gez,0).*fem.gew,2);
    tu3 = sum(feval(pde.Et3,fem.gex,fem.gey,fem.gez,0).*fem.gew,2);
    upretmp = tu1.*tgt(:,1) + tu2.*tgt(:,2) + tu3.*tgt(:,3);
    upre = uppre + deltaT*upretmp;
end

UHfit =zeros(size(mesh.e,1),ceil(Ntime/10)+2);
NumStep = 0;
CurrentT = 0;
SaveCount = 0;

tID1 = (mesh.tLoc == 1);
tID2 = (mesh.tLoc == 2);
e1tmp = unique(reshape(mesh.t_e(tID1,:),[],1));
e2tmp = unique(reshape(mesh.t_e(tID2,:),[],1));
eInt = intersect(e1tmp,e2tmp);
eid1 = setdiff(e1tmp,eInt);
eid2 = setdiff(e2tmp,eInt);
alpha = am*ones(size(mesh.e,1),1);
alpha(eid2) = ap;
alpha(eInt) = (am+ap)/2;
beta = bm*ones(size(mesh.e,1),1);
beta(eid2) = bp;
beta(eInt) = (bm+bp)/2;
NEdof = size(A,1);

while CurrentT < TEND

    NumStep = NumStep + 1;
    CurrentT = CurrentT + deltaT;
    fm1Current = @(x,y,z) pde.fm1(x,y,z,CurrentT);
    fm2Current = @(x,y,z) pde.fm2(x,y,z,CurrentT);
    fm3Current = @(x,y,z) pde.fm3(x,y,z,CurrentT);
    fp1Current = @(x,y,z) pde.fp1(x,y,z,CurrentT);
    fp2Current = @(x,y,z) pde.fp2(x,y,z,CurrentT);
    fp3Current = @(x,y,z) pde.fp3(x,y,z,CurrentT);
    rhsF1 = globNedFitRHS3D(fm1Current,fp1Current, mesh, fem, 0, 1);
    rhsF2 = globNedFitRHS3D(fm2Current,fp2Current, mesh, fem, 0, 2);
    rhsF3 = globNedFitRHS3D(fm3Current,fp3Current, mesh, fem, 0, 3);
    rhsFCurrent = rhsF1 + rhsF2 + rhsF3;

    if BDCond == 0
        tuCurrent = zeros(size(mesh.e,1),1);
    elseif BDCond ~= 0
        exactu1Current = @(x,y,z) pde.E1(x,y,z,CurrentT);
        exactu2Current = @(x,y,z) pde.E2(x,y,z,CurrentT);
        exactu3Current = @(x,y,z) pde.E3(x,y,z,CurrentT);
        tu1 = sum(feval(exactu1Current,fem.gex,fem.gey,fem.gez).*fem.gew,2);
        tu2 = sum(feval(exactu2Current,fem.gex,fem.gey,fem.gez).*fem.gew,2);
        tu3 = sum(feval(exactu3Current,fem.gex,fem.gey,fem.gez).*fem.gew,2);
        tgt = mesh.p(mesh.e(:,2),:) - mesh.p(mesh.e(:,1),:);
        tgt = tgt./sum(tgt.^2,2).^(1/2);
        tuCurrent = tu1.*tgt(:,1) + tu2.*tgt(:,2) + tu3.*tgt(:,3);
    end

    ub = tuCurrent;
    ub(fem.mapper) = 0;
    rhsB = Atotal*ub;
    JCurrent = rhsFCurrent - rhsB;
    %JCurrent(isBdEdge) = tuCurrent(isBdEdge);

    fcurrent = Me*(2*upre - uppre)/deltaT^2 + Ms*uppre/(2*deltaT) - S*uppre/2 + JCurrent;
    fcurrent(isBdEdge) = tuCurrent(isBdEdge);

    option.outsolver = 'cg';
    option.alpha = alpha;
    option.beta = beta;
    option.solver = 'amg';
    edge = mesh.e;
    option.smoother = 'BD';
    option.blklevel = 0;
    option.blkId = NEdof;
    [x,info] = amgMaxwellinterface2(A,fcurrent,mesh.p,edge,option);
    uppre = upre;
    upre = x;

    if mod(NumStep,40) == 0
        UHfit(:,SaveCount+1) = uppre;
        UHfit(:,SaveCount+2) = upre;
        SaveCount = SaveCount+2;
        save('UHfit1','UHfit','-v7.3');
        disp('save data')
    end

end

save('UHfit1','UHfit');

repo: ifem


text:
% STD_READSPECGRAM - returns the stored mean power spectrogram for an ICA component
%                  or a data channel in a specified dataset.  The spectrogram is
%                  assumed to have been saved in a Matlab file,
%                  "[dataset_name].datspecgram", in the same
%                  directory as the dataset file. If this file doesn't exist,
%                  use STD_SPECGRAM to create it.
% Usage:
%  >> [spec, times, freqs] = std_readspecgram(ALLEEG, setindx, component, timerange, freqrange);
%
% Inputs:
%   ALLEEG     - a vector of dataset EEG structures (may also be one dataset).
%                Must contain the dataset of interest (the 'setindx' below).
%   setindx    - [integer] an index of an EEG dataset in the ALLEEG
%                structure for which to read a component spectrum.
%   component  - [integer] index of the component in the selected EEG dataset
%                for which to return the spectrum
%   freqrange  - [min max in Hz] frequency range to return
%
%
% Outputs:
%   spec      - the log-power spectrum of the requested ICA component in the
%               specified dataset (in dB)
%   freqs     - vector of spectral frequencies (in Hz)
%
%  See also  STD_SPEC, POP_PRECLUST, STD_PRECLUST
%
% Authors:  Arnaud Delorme, SCCN, INC, UCSD, February, 2008

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

function [X, t, f] = std_readspecgram(ALLEEG, abset, comp, timerange, freqrange, rmsubjmean);

if nargin < 4
    timerange = [];
end
if nargin < 5
    freqrange = [];
end

X = [];

if iscell(comp)
    % find channel indices list
    % -------------------------
    chanind  = [];
    tmpchanlocs = ALLEEG(abset).chanlocs;
    chanlabs = lower({ tmpchanlocs.labels });
    for index = 1:length(comp)
        tmp = strmatch(lower(comp{index}), chanlabs, 'exact');
        if isempty(tmp)
            error([ 'Channel ''' comp{index} ''' not found in dataset ' int2str(abset)]);
        else
            chanind = [ chanind tmp ];
        end
    end
    filename = fullfile( ALLEEG(abset).filepath,[ ALLEEG(abset).filename(1:end-3) 'datspecgram']);
    prefix = 'chan';
    inds   = chanind;
elseif comp(1) < 0
    filename = fullfile( ALLEEG(abset).filepath,[ ALLEEG(abset).filename(1:end-3) 'datspecgram']);
    prefix = 'chan';
    inds   = -comp;
else
    filename = fullfile( ALLEEG(abset).filepath,[ ALLEEG(abset).filename(1:end-3) 'icaspecgram']);
    prefix = 'comp';
    inds   = comp;
end

for k=1:length(inds)
    try,
        warning backtrace off;
        erpstruct = load( '-mat', filename, [ prefix int2str(inds(k)) ], 'freqs', 'times');
        warning backtrace on;
    catch
        error( [ 'Cannot read file ''' filename '''' ]);
    end

    tmpdat    = getfield(erpstruct, [ prefix int2str(inds(k)) ]);
    if k == 1
        X = zeros([size(tmpdat) length(comp)]);
    end
    X(:,:,k) = tmpdat;
    f        = getfield(erpstruct, 'freqs');
    t        = getfield(erpstruct, 'times');
end

% select frequency range of interest
% ----------------------------------
if ~isempty(freqrange)
    maxind = max(find(f <= freqrange(end)));
    minind = min(find(f >= freqrange(1)));
else
    %if not, use whole spectrum
    maxind = length(f);
    minind = 1;
end
f = f(minind:maxind);
X = X(minind:maxind,:,:);

% select time range of interest
% -----------------------------
if ~isempty(timerange)
    maxind = max(find(t <= timerange(end)));
    minind = min(find(t >= timerange(1)));
else
    %if not, use whole spectrum
    maxind = length(t);
    minind = 1;
end
t = t(minind:maxind);
X = X(:,minind:maxind,:);

return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
scores=mod(score,1000);
d1=mod(scores,10);
d2=(mod(scores,100)-d1)/10;
d3=(mod(scores,1000)-d2*10-d1)/100;
if scores<10
    d2=10;
end
if scores<100
    d3=10;
end
set_digit(his,digs{3},d1);
set_digit(his,digs{2},d2);
set_digit(his,digs{1},d3);


text:
% EX_MAXWELL_SRC_CUBE: solve Maxwell source problem in a thick ring.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_thick_ring.txt';

% Type of boundary conditions
problem_data.nmnn_sides   = [1 2];
problem_data.drchlt_sides = [3 4 5 6];

% Physical parameters
problem_data.c_mass  = @(x, y, z) ones(size(x));
problem_data.c_stiff = @(x, y, z) ones(size(x));

% Source and boundary terms
problem_data.f = @(x, y, z) cat(1, ...
                    reshape (2*sin(y), [1, size(x)]), ...
                    reshape (2*sin(x), [1, size(x)]), ...
                    reshape (x .* y, [1, size(x)]));
problem_data.g = @(x, y, z, ind) test_maxwell_thick_ring_g_nmnn (x, y, z, ind);
problem_data.h = @(x, y, z, ind) cat(1, ...
                    reshape (sin(y), [1, size(x)]), ...
                    reshape (sin(x), [1, size(x)]), ...
                    reshape (x .* y, [1, size(x)]));

% Exact solution (optional)
problem_data.uex     = @(x, y, z) cat(1, ...
                    reshape (sin(y), [1, size(x)]), ...
                    reshape (sin(x), [1, size(x)]), ...
                    reshape (x .* y, [1, size(x)]));
problem_data.curluex = @(x, y, z) cat(1, ...
                    reshape (x, [1, size(x)]), ...
                    reshape (-y, [1, size(x)]), ...
                    reshape (cos(x) - cos(y), [1, size(x)]));

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.degree     = [2 2 2];     % Degree of the bsplines
method_data.regularity = [1 1 1];     % Regularity of the splines
method_data.nsub       = [3 3 3];     % Number of subdivisions
method_data.nquad      = [3 3 3];     % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space, u] = solve_maxwell_src (problem_data, method_data);

% 4) POST-PROCESSING
% 4.1) EXPORT TO PARAVIEW
output_file = 'maxwell_thick_ring_Deg2_Reg1_Sub3';

vtk_pts = {linspace(0, 1, 15), linspace(0, 1, 15), linspace(0, 1, 15)};
fprintf ('The result is saved in the file %s \n \n', output_file);
sp_to_vtk (u, space, geometry, vtk_pts, output_file, 'u')

% 4.2) Comparison with the exact solution
[error_hcurl, error_l2] = ...
    sp_hcurl_error (space, msh, u, problem_data.uex, problem_data.curluex)

%!demo
%! ex_maxwell_src_thick_ring

%!test
%! problem_data.geo_name = 'geo_thick_ring.txt';
%! problem_data.nmnn_sides   = [1 2];
%! problem_data.drchlt_sides = [3 4 5 6];
%! problem_data.c_mass  = @(x, y, z) ones(size(x));
%! problem_data.c_stiff = @(x, y, z) ones(size(x));
%! problem_data.f = @(x, y, z) cat(1, ...
%!                     reshape (2*sin(y), [1, size(x)]), ...
%!                     reshape (2*sin(x), [1, size(x)]), ...
%!                     reshape (x .* y, [1, size(x)]));
%! problem_data.g = @(x, y, z, ind) test_maxwell_thick_ring_g_nmnn (x, y, z, ind);
%! problem_data.h = @(x, y, z, ind) cat(1, ...
%!                     reshape (sin(y), [1, size(x)]), ...
%!                     reshape (sin(x), [1, size(x)]), ...
%!                     reshape (x .* y, [1, size(x)]));
%! problem_data.uex     = @(x, y, z) cat(1, ...
%!                     reshape (sin(y), [1, size(x)]), ...
%!                     reshape (sin(x), [1, size(x)]), ...
%!                     reshape (x .* y, [1, size(x)]));
%! problem_data.curluex = @(x, y, z) cat(1, ...
%!                     reshape (x, [1, size(x)]), ...
%!                     reshape (-y, [1, size(x)]), ...
%!                     reshape (cos(x) - cos(y), [1, size(x)]));
%! method_data.degree     = [2 2 2];     % Degree of the bsplines
%! method_data.regularity = [1 1 1];     % Regularity of the splines
%! method_data.nsub       = [3 3 3];     % Number of subdivisions
%! method_data.nquad      = [3 3 3];     % Points for the Gaussian quadrature rule
%! [geometry, msh, space, u] = solve_maxwell_src (problem_data, method_data);
%! [error_hcurl, error_l2] = sp_hcurl_error (space, msh, u, problem_data.uex, problem_data.curluex);
%! assert (msh.nel, 27)
%! assert (space.ndof, 300)
%! assert (error_l2, 0.0643150154230899, 1e-14)
%! assert (error_hcurl, 0.141137345617213, 1e-14)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function varargout = biharm(varargin)
%BIHARM   Biharmonic operator of a CHEBFUN2.
%   B = BIHARM(F) returns a CHEBFUN2 representing the biharmonic operator
%   applied to F.
%
%   This is shorthand for BIHARMONIC(F).
%
% See also CHEBFUN2/BIHARMONIC.

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

% Call SEPARABLEAPPROX/BIHARM:
[varargout{1:nargout}] = biharm@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [res]=tt_mm(mat1,mat2)
%Matrix-by-matrix product in TT1.0 format
%   [RES]=TT_MM(MAT1,MAT2) Matrix MAT1 by matrix MAT2 product in the
%   TT-format. Please avoid its usage: it will be removed in
%   future releases. Use * operator from the object-oriented version.
%
%
% TT-Toolbox 2.2, 2009-2011
%
%This is TT Toolbox, written by Ivan Oseledets, Olga Lebedeva
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------
d=size(mat1,1);
res=cell(d,1);
k=1;
for s=1:2
n1=size(mat1{k},1);
m1=size(mat1{k},2);
rm=size(mat1{k},3);
n2=size(mat2{k},1);
m2=size(mat2{k},2);
rv=size(mat2{k},3);
%mat1{1} is (n1 x m1 x rm , mat2{1} is n2xm2xrm, result is n1xm2xrv1*rv)
%res{1}=reshape(permute(mat1{1},[1,3,2],[n1*rm1]),m1)*reshape(mat2{1},[n2,m2*rv2]);
res{k}=reshape(permute(mat1{k},[1,3,2]),[n1*rm,m1])*reshape(mat2{k},[n2,m2*rv]);
%res{1} is n1 x rm x m2 x rv
res{k}=reshape(res{k},[n1,rm,m2,rv]);
res{k}=permute(res{k},[1,3,2,4]);
res{k}=reshape(res{k},[n1,m2,rm*rv]);
if ( k == 1 )
 k=d;
end
end
for i=2:d-1
  n1=size(mat1{i},1);
  m1=size(mat1{i},2);
  rm1=size(mat1{i},3);
  rm2=size(mat1{i},4);
  n2=size(mat2{i},1);
  m2=size(mat2{i},2);
  rv1=size(mat2{i},3);
  rv2=size(mat2{i},4);
  %mat1{i} is n1 x m1 x rm1 x rm2
  %mat2{i} is n2 x m2 x rv1 x rv2
  %n1 x rm1 x rm2 x m1 * n2 x m2 x rv1 x rv2
  %n1 x rm1 x rm2 x m2 x rv1 x rv2
  %Need: n1 x m2 x rm1 x rv1 x rm2 x rv2
  %Permut: 1 4 2 5 3 6
  res{i}=reshape(permute(mat1{i},[1,3,4,2]),[n1*rm1*rm2,m1])*reshape(mat2{i},[n2,m2*rv1*rv2]);
  res{i}=reshape(res{i},[n1,rm1,rm2,m2,rv1,rv2]);
  res{i}=permute(res{i},[1,4,2,5,3,6]);
  res{i}=reshape(res{i},[n1,m2,rm1*rv1,rm2*rv2]);
end
return
end

author: oseledets

repo: TT-Toolbox


text:
function tru = getTAFromLineOfNodes(nVect, sma, ecc, inc, raan, arg, gmu)

    func = @(tru) dangBetweenRVectAndLineToAscNode(nVect, sma, ecc, inc, raan, arg, tru, gmu);
    if(ecc < 1)
        tru = fminbnd(func, 0, 2*pi, optimset('TolX',eps));
    else
        maxTA = computeTrueAFromRadiusEcc(Inf, sma, ecc);
        tru = fminbnd(func, -maxTA, maxTA, optimset('TolX',eps));
    end
end

function angle = dangBetweenRVectAndLineToAscNode(nVect, sma, ecc, inc, raan, arg, tru, gmu)
    [rVect,~]=getStatefromKepler(sma, ecc, inc, raan, arg, tru, gmu);
    angle = dang(nVect,rVect);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% CVX: Additional functions added by CVX.
%
%   These functions have been provided to expand the variety of constraints
%   and objectives that can be specified in CVX models. But in fact, they
%   can be used with numeric arguments *outside* of CVX as well. The help
%   text for each of these functions contains general information about the
%   computations it performs, as well as specific information about its
%   proper use in CVX models, as dictated by its convexity/concavity and
%   monotonicity properties.
%
%   Those functions marked with a (*), as well as the exponential and
%   logarithm functions, are supported using a "successive approximation"
%   approach: that is, the solver is called multiple times to refine the
%   solution to the required accuracy. Thus models using these functions
%   should be expected to run more slowly than models of comparable size
%   that do not. See the CVX user guide for details.
%
%   A number of Matlab's built-in functions have been extended to provide
%   CVX support; for example,
%     abs, exp(*), log (*), max, min, norm, prod, sqrt
%   For a full list, type "help cvx/builtins".
%
%   berhu             - Reverse Huber penalty function.
%   det_inv           - Determinant of the inverse of an SPD matrix.
%   det_root2n        - 2nth-root of the determinant of an SPD matrix.
%   det_rootn         - nth-root of the determinant of an SPD matrix.
%   entr              - Scalar entropy. (*)
%   geo_mean          - Geometric mean. (*)
%   huber             - Huber penalty function.
%   huber_circ        - Circularly symmetric version of the Huber penalty.
%   huber_pos         - Monotonic Huber-style function.
%   inv_pos           - Reciprocal of a positive quantity.
%   kl_div            - Scalar Kullback-Leibler distance. (*)
%   lambda_max        - Maximum eigenvalue of a symmetric matrix.
%   lambda_min        - Minimum eigenvalue of a symmetric matrix.
%   log_det           - Logarithm of the determinant. (*)
%   log_normcdf       - Logarithm of the normal CDF. (approximation)
%   log_sum_exp       - log(sum(exp(x))). (*)
%   logsumexp_sdp     - SDP-based approximation of log(sum(exp(x))).
%   matrix_frac       - Matrix fractional function.
%   norm_largest      - Sum of the k largest magnitudes of a vector.
%   norm_nuc          - Nuclear norm of a matrix.
%   norms             - Computation of multiple vector norms.
%   norms_largest     - Computation of multiple norm_largest() norms.
%   poly_env          - Convex or concave envelope of a polynomial.
%   polyval_trig      - Evaluate a trigonometric polynomial.
%   pos               - Positive part.
%   pow_p             - Nonnegative branches of the power function.
%   pow_pos           - Convex/concave branches of the power function.
%   pow_abs           - Absolute value raised to a fixed power.
%   quad_form         - Quadratic form.
%   quad_over_lin     - Sum of squares over linear.
%   quad_pos_over_lin - Sum of squares of positives over linear.
%   rel_entr          - Scalar relative entropy. (*)
%   sigma_max         - Maximum singular value.
%   square            - Square.
%   square_abs        - Square of absolute value.
%   square_pos        - Square of positive part.
%   sum_largest       - Sum of the largest k values of a vector.
%   sum_smallest      - Sum of the smallest k elements of a vector.
%   sum_square        - Sum of squares.
%   sum_square_abs    - sum of squares of absolute values.
%   sum_square_pos    - Sum of squares of positive parts.
%   trace_inv         - Trace of the inverse of a PSD matrix.
%   trace_sqrtm       - Trace of the square root of a PSD matrix.
%   vec               - Vectorize.

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function op = linop_compose( varargin )
%LINOP_COMPOSE Composes two TFOCS linear operators
%OP = LINOP_COMPOSE( OP1, OP2, ..., OPN )
%    Constructs a TFOCS-compatible linear operator from the composition of
%    two or more linear operators and/or matrices. That is,
%       OP(x,1) = OP1(OP2(...(OPN(x,1),...,1),1)
%       OP(x,2) = OPN(...(OP2(OP1(x,2),2)...,2)
%    If matrices are supplied, they must be real; to include complex
%    matrices, convert them first to linear operators using LINOP_MATRIX.

if nargin == 0,
    error( 'Not enough input arguments.' );
end
sz = { [], [] };
for k = 1 : nargin,
    tL = varargin{k};
    if isempty(tL) || ~isa(tL,'function_handle') && ~isnumeric(tL) || ndims(tL) > 2,
        error( 'Arguments must be linear operators, scalars, or matrices.' );
    elseif isnumeric(tL),
        if ~isreal(tL),
            error( 'S or scalar arguments must be real.' );
        elseif numel(tL) == 1,
            tL = linop_scale( tL );
        else
            tL = linop_matrix( tL );
        end
        varargin{k} = tL;
    end
    try
        tsz = tL([],0);
    catch
        error( 'Arguments must be linear operators, scalars, or matrices.' );
    end
    if isempty(tsz)     % i.e. output of linop_identity is []
        tsz = { [], [] };
    end
    if isnumeric(tsz),
        tsz = { [tsz(2),1], [tsz(1),1] };
    end

    % convert [n1;n2] to [n1,n2] if necessary:
    for kk = 1:2
        %if iscolumn( tsz{kk} )
            %tsz{kk} = tsz{kk}.';
        %end
        tsz{kk} = tsz{kk}(:).';
    end

    if ~isempty(sz{1}) && ~isempty(tsz{2}) && ~isequal(tsz{2},sz{1}),
        for kk = 1:min( numel(tsz{2}), numel( sz{1} ) )
            fprintf('Found incompatible sizes: %d ~= %d\n', tsz{2}(kk), sz{1}(kk) );
        end
        error( 'Incompatible dimensions in linear operator composition.' );
    end
    if ~isempty(tsz{1}),
        sz{1} = tsz{1};
    end
    if isempty(sz{2}),
        sz{2} = tsz{2};
    end
end
% Explanation of above code:
% suppose have three inputs, opA, opB, opC; with sizes szA, szB, szC
%   where opA: szA{1} --> szA{2}
%   so we need szC{2} == szB{1} and szB{2} = szA{1}
%
% so  sz{2} = szA{2}
% and sz{1} = szC{1}

if nargin == 1,
    op = varargin{1};
else
    op = @(x,mode)linop_compose_impl( varargin, sz, x, mode );
end

function y = linop_compose_impl( ops, sz, y, mode )
switch mode,
    case 0,
        y = sz;
    case 1,
        for k = numel(ops) : -1 : 1,
            y = ops{k}( y, 1 );
        end
    case 2,
        for k = 1 : numel(ops),
            y = ops{k}( y, 2 );
        end
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
% VL_DEMO_COVDET Demo: VL_COVDET()

% --------------------------------------------------------------------
%                                                               Basics
% --------------------------------------------------------------------

im = vl_impattern('roofs1') ;
im = im(end-128:end,128:320,:) ;

figure(1) ; clf ;
image(im) ; axis image off ;
vl_demo_print('covdet_basic_image') ;

imgs = im2single(rgb2gray(im)) ;
frames = vl_covdet(imgs, 'verbose') ;

hold on ;
vl_plotframe(frames) ;
vl_demo_print('covdet_basic_frames') ;

% --------------------------------------------------------------------
%                                                    Affine adaptation
% --------------------------------------------------------------------

frames = vl_covdet(imgs, 'estimateAffineShape', true, 'verbose') ;

figure(2) ; clf ;
image(im) ; axis image off ; hold on ;
vl_plotframe(frames) ;
vl_demo_print('covdet_affine_frames',.8) ;

% --------------------------------------------------------------------
%                                              Estimating orientations
% --------------------------------------------------------------------

frames = vl_covdet(imgs, 'estimateOrientation', true, 'verbose') ;

figure(3) ; clf ;
image(im) ; axis image off ; hold on ;
vl_plotframe(frames) ;
vl_demo_print('covdet_oriented_frames',.8) ;

% --------------------------------------------------------------------
%                                                   Extracting patches
% --------------------------------------------------------------------

[frames, patches] = vl_covdet(imgs, 'descriptor', 'patch') ;

figure(4) ; clf ;
w = sqrt(size(patches,1)) ;
vl_imarraysc(reshape(patches(:,1:10*10), w,w,[])) ;
axis image off ; hold on ; colormap gray ;
vl_demo_print('covdet_patches') ;

[frames, patches] = vl_covdet(imgs, ...
                              'descriptor', 'patch' ,...
                              'estimateAffineShape', true, ...
                              'estimateOrientation', false) ;

figure(5) ; clf ;
w = sqrt(size(patches,1)) ;
vl_imarraysc(reshape(patches(:,1:10*10), w,w,[])) ;
axis image off ; hold on ; colormap gray ;
vl_demo_print('covdet_affine_patches') ;

% --------------------------------------------------------------------
%                                                  Different detectors
% --------------------------------------------------------------------

names = {'DoG', 'Hessian', ...
         'HarrisLaplace', 'HessianLaplace', ...
         'MultiscaleHarris', 'MultiscaleHessian'} ;
figure(6) ; clf ;
for i = 1:numel(names)
  frames = vl_covdet(imgs, 'method', names{i}) ;

  vl_tightsubplot(3,2,i, 'margintop',0.025, 'marginright', 0.01) ;
  imagesc(im) ; axis image off ;
  hold on ;
  vl_plotframe(frames) ;
  title(names{i}) ;
end

vl_figaspect(3/4) ;
vl_demo_print('covdet_detectors',.9) ;

% --------------------------------------------------------------------
%                                                        Custom frames
% --------------------------------------------------------------------

delta = 15 ;
xr = delta:delta:size(im,2)-delta+1 ;
yr = delta:delta:size(im,1)-delta+1 ;
[x,y] = meshgrid(xr,yr) ;
frames = [x(:)'; y(:)'] ;
frames(end+1,:) = delta/2 ;

[frames, patches] = vl_covdet(imgs, ...
                              'frames', frames, ...
                              'estimateAffineShape', true, ...
                              'estimateOrientation', true) ;

figure(7) ; clf ;
imagesc(im) ;
axis image off ; hold on ; colormap gray ;
vl_plotframe(frames) ;
vl_demo_print('covdet_custom_frames',.8) ;

% --------------------------------------------------------------------
%                                                         Scale spaces
% --------------------------------------------------------------------

[frames, descrs, info] = vl_covdet(imgs) ;

info

figure(8) ; clf ;
vl_plotss(info.gss) ;
colormap gray ;
vl_figaspect(2) ;
vl_demo_print('covdet_gss',.8) ;

figure(9) ; clf ;
vl_plotss(info.css) ;
colormap gray ;
vl_figaspect(2) ;
vl_demo_print('covdet_css',.8) ;

figure(10) ; clf ;
subplot(1,2,1) ;
hist([info.peakScores],10) ;
xlabel('Peak Score') ;
ylabel('Occurences') ;
grid on ;

subplot(1,2,2) ;
hist([info.edgeScores],10) ;
xlabel('Edge Score') ;
ylabel('Occurences') ;
grid on ;

vl_figaspect(2) ;
vl_demo_print('covdet_scores',.9) ;

author: Cloud-CV

repo: object-proposals


text:
function appendYALMIPvariables(lmi_variables,mt,variabletype,hashed_monoms,current_hash);

if nargin == 1
    [mt,variabletype,hashed_monoms,current_hash] = yalmip('monomtable');
end

% Update monomtable and pre-calculated variable type
n_mt = size(mt,1);
m_mt = size(mt,2);
newmt = [];
if min(lmi_variables)>m_mt % New variables
    if size(mt,1)~=size(mt,2)
        mt(size(mt,1),size(mt,1))=0;
    end
    % This was faster before. However in recent versions of matlab, there
    % is a compiled version of blkdiag available
    % fill=spalloc(size(mt,1),length(lmi_variables),0);
    % mt=[mt fill;fill' speye(length(lmi_variables))];
    if isempty(mt)
        mt = speye(length(lmi_variables));
        newmt = mt;
    elseif length(lmi_variables)==1
        % Slightly faster than general case
       newmt = sparse(1);
       mt(end+1,end+1) = newmt;
    else
        newmt = speye(length(lmi_variables));
        mt=blkdiag(mt,newmt);
    end
else
    mt(lmi_variables,lmi_variables) = speye(length(lmi_variables));
end
variabletype(1,size(mt,1)) = 0;
if ~isempty(newmt)
    new_hash = 3*rand_hash(size(mt,2),size(newmt,2),1);
    hashed_monoms = [hashed_monoms;full(newmt*new_hash)];
    current_hash = [current_hash;new_hash];
    yalmip('setmonomtable',mt,variabletype,hashed_monoms,current_hash);
else
    yalmip('setmonomtable',mt,variabletype);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [neighbors nns] = sameclass_nn(data, nn, W)

if (exist('W', 'var'))
  if (~iscell(W))
    B1 = W * data.Xtraining;
  else
    resp1 = compute_NN_output3(W, data.Xtraining);
    B1 = resp1{end};
  end

  B1 = logical(single(B1 > 0));
  B1 = compactbit(B1);
  [nw1 ~] = size(B1);
  nb = nw1*8;
end

nl = zeros(numel(data.labels(:)), 1);
for l = data.labels(:)'
  nl(l+1) = sum(data.Ltraining == l);
end

nns = zeros(data.Ntraining,1);
neighbors = zeros(min(nn, min(nl)-1), data.Ntraining, 'single');

for l = data.labels(:)'
  which_l = data.Ltraining == l;
  id_l = find(which_l);
  nn_l = min(nn, nl(l+1)-1);
  nns(which_l) = nn_l;
  if (~exist('W', 'var'))
    [nnTraining_l ~] = yael_nn(single(data.Xtraining(:,which_l)), ...
			       single(data.Xtraining(:,which_l)), nn_l+1 );
  else
    [btmp nnTraining_l ctmp dtmp] = linscan_sorted_multi_mex(B1(:,which_l), B1(:,which_l), nl(l+1), nb, ceil(nb/2), nn_l+1, 4);
  end

  for i=1:numel(id_l)
    neighbors(1:nn_l,id_l(i)) = id_l(nnTraining_l(2:end, i));
  end
end

author: norouzi

repo: hdml

save_path: github-repos/MATLAB/norouzi-hdml


text:
function lax_wendroff_2d_display ( )

%*****************************************************************************80
%
%% LAX_WENDROFF_2D_DISPLAY illustrates the Lax-Wendroff procedure in 2D.
%
%  Discussion:
%
%    A sequence of ball and stick drawings is displayed in steps 0 through 9.
%
%    Hitting return causes the next item to be drawn.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2003
%
%  Author:
%
%    John Burkardt
%
  for frame = 0 : 9
%
%  Clear the plotting frame.
%
    clf
%
%  Make all the subsequent plotting "cumulative".
%
    hold on

    if ( 0 == frame )
      title ( 'Lax-Wendroff 2D scheme' )
    end
%
%  Central nodes on the blue level.
%
    if ( 1 <= frame )
      plot3 (  0.0,  0.0, 0.0, 'o', 'MarkerFaceColor', 'b', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
    end
    if ( frame == 1 )
      title ( 'A node where we have the solution.' )
    end
%
%  The region associated with the central node, blue level.
%
    if ( 2 <= frame )
      plot3 (  [ -0.5,  0.5 ],  [ -0.5, -0.5 ],  [ 0.0, 0.0 ], 'b', 'LineWidth', 3 )
      plot3 (  [  0.5,  0.5 ],  [ -0.5,  0.5 ],  [ 0.0, 0.0 ], 'b', 'LineWidth', 3 )
      plot3 (  [  0.5, -0.5 ],  [  0.5,  0.5 ],  [ 0.0, 0.0 ], 'b', 'LineWidth', 3 )
      plot3 (  [ -0.5, -0.5 ],  [  0.5, -0.5 ],  [ 0.0, 0.0 ], 'b', 'LineWidth', 3 )
    end
    if ( frame == 2 )
      title ( 'The region associated with the node.' )
    end
%
%  Neighbor nodes on the blue level.
%
    if ( 3 <= frame )
      plot3 ( -1.0,  0.0, 0.0, 'o', 'MarkerFaceColor', 'b', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
      plot3 (  1.0,  0.0, 0.0, 'o', 'MarkerFaceColor', 'b', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
      plot3 (  0.0, -1.0, 0.0, 'o', 'MarkerFaceColor', 'b', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
      plot3 (  0.0,  1.0, 0.0, 'o', 'MarkerFaceColor', 'b', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
    end
    if ( frame == 3 )
      title ( 'The neighboring nodes.' )
    end
%
%  Midside nodes on the blue level.
%
    if ( 4 <= frame )
      plot3 ( -0.5,  0.0, 0.0, 'o', 'MarkerFaceColor', 'b', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
      plot3 (  0.5,  0.0, 0.0, 'o', 'MarkerFaceColor', 'b', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
      plot3 (  0.0, -0.5, 0.0, 'o', 'MarkerFaceColor', 'b', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
      plot3 (  0.0,  0.5, 0.0, 'o', 'MarkerFaceColor', 'b', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
    end
    if ( frame == 4 )
      plot3 (  [ -1.0, +1.0 ],  [  0.0,  0.0 ],  [ 0.0, 0.0 ], 'b', 'LineWidth', 3 )
      plot3 (  [  0.0,  0.0 ],  [ -1.0,  1.0 ],  [ 0.0, 0.0 ], 'b', 'LineWidth', 3 )
    end
    if ( frame == 4 )
      title ( 'Values at midsides by averaging.' )
    end
%
%  Center node on the green level.
%  Line from center node on blue to center node on green.
%
    if ( 5 <= frame )
      plot3 ( -0.5,  0.0, 0.5, 'o', 'MarkerFaceColor', 'g', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
      plot3 (  0.5,  0.0, 0.5, 'o', 'MarkerFaceColor', 'g', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
      plot3 (  0.0, -0.5, 0.5, 'o', 'MarkerFaceColor', 'g', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
      plot3 (  0.0,  0.5, 0.5, 'o', 'MarkerFaceColor', 'g', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
      plot3 (  [ -0.5,  0.5 ],  [ -0.5, -0.5 ],  [ 0.5, 0.5 ], 'g', 'LineWidth', 3 )
      plot3 (  [  0.5,  0.5 ],  [ -0.5,  0.5 ],  [ 0.5, 0.5 ], 'g', 'LineWidth', 3 )
      plot3 (  [  0.5, -0.5 ],  [  0.5,  0.5 ],  [ 0.5, 0.5 ], 'g', 'LineWidth', 3 )
      plot3 (  [ -0.5, -0.5 ],  [  0.5, -0.5 ],  [ 0.5, 0.5 ], 'g', 'LineWidth', 3 )
    end
    if ( frame == 5 )
      plot3 (  [ -0.5, -0.5 ],  [  0.0,  0.0 ],  [ 0.0, 0.5 ], 'k', 'LineWidth', 3 )
      plot3 (  [  0.5,  0.5 ],  [  0.0,  0.0 ],  [ 0.0, 0.5 ], 'k', 'LineWidth', 3 )
      plot3 (  [  0.0,  0.0 ],  [ -0.5, -0.5 ],  [ 0.0, 0.5 ], 'k', 'LineWidth', 3 )
      plot3 (  [  0.0,  0.0 ],  [  0.5,  0.5 ],  [ 0.0, 0.5 ], 'k', 'LineWidth', 3 )
    end
    if ( frame == 5 )
      title ( 'Values at midside, half time step.' )
    end
%
%  Lines from midside nodes to center.
%
    if ( 6 == frame )
      plot3 (  [  0.5,  0.0 ],  [  0.0,  0.0 ],  [ 0.5, 0.5 ], 'g', 'LineWidth', 3 )
      plot3 (  [ -0.5,  0.0 ],  [  0.0,  0.0 ],  [ 0.5, 0.5 ], 'g', 'LineWidth', 3 )
      plot3 (  [  0.0,  0.0 ],  [  0.5,  0.0 ],  [ 0.5, 0.5 ], 'g', 'LineWidth', 3 )
      plot3 (  [  0.0,  0.0 ],  [ -0.5,  0.0 ],  [ 0.5, 0.5 ], 'g', 'LineWidth', 3 )
    end
    if ( frame == 6 )
      title ( 'Estimate fluxes' )
    end
%
%  Carry neighbor nodes to green level, and draw lines.
%
    if ( 7 <= frame )
      plot3 (  0.0,  0.0, 0.5, 'o', 'MarkerFaceColor', 'g', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
    end
    if ( frame == 7 )
      title ( 'Estimate derivative at node, half time step.' )
    end
%
%  Advance solution to full time.
%
    if ( 8 <= frame )
      plot3 (  [  0.0,  0.0 ],  [  0.0,  0.0 ],  [ 0.0, 1.0 ], 'k', 'LineWidth', 3 )
      plot3 (  0.0,  0.0, 1.0, 'o', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
      plot3 (  [ -0.5,  0.5 ],  [ -0.5, -0.5 ],  [ 1.0, 1.0 ], 'r', 'LineWidth', 3 )
      plot3 (  [  0.5,  0.5 ],  [ -0.5,  0.5 ],  [ 1.0, 1.0 ], 'r', 'LineWidth', 3 )
      plot3 (  [  0.5, -0.5 ],  [  0.5,  0.5 ],  [ 1.0, 1.0 ], 'r', 'LineWidth', 3 )
      plot3 (  [ -0.5, -0.5 ],  [  0.5, -0.5 ],  [ 1.0, 1.0 ], 'r', 'LineWidth', 3 )
    end
    if ( frame == 8 )
      title ( 'Take full step to new time.' )
    end
%
%  Advance all data.
%
    if ( 9 <= frame )
      plot3 (  [ -1.0, -1.0 ],  [  0.0,  0.0 ],  [ 0.0, 1.0 ], 'k', 'LineWidth', 3 )
      plot3 (  [  1.0,  1.0 ],  [  0.0,  0.0 ],  [ 0.0, 1.0 ], 'k', 'LineWidth', 3 )
      plot3 (  [  0.0,  0.0 ],  [ -1.0, -1.0 ],  [ 0.0, 1.0 ], 'k', 'LineWidth', 3 )
      plot3 (  [  0.0,  0.0 ],  [  1.0,  1.0 ],  [ 0.0, 1.0 ], 'k', 'LineWidth', 3 )
      plot3 ( -1.0,  0.0, 1.0, 'o', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
      plot3 (  1.0,  0.0, 1.0, 'o', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
      plot3 (  0.0, -1.0, 1.0, 'o', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
      plot3 (  0.0,  1.0, 1.0, 'o', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'k', 'MarkerSize', 12 )
    end
    if ( frame == 9 )
      title ( 'Advance all data to new time.' )
    end

    xlabel ( '- X -' )
    ylabel ( '- Y -' )
    zlabel ( '- Time -' )

    axis ( [ -1, +1, -1, +1, 0, 1] )
    view ( 3 )
%
%  "Release" the plotting screen.
%
    hold off

    pause

  end

  return
end

repo: jburkardt-m


text:
function toolboxSfmGenDoc()
% Generate documentation, must run from dir toolbox.
%
% Requires external/m2html to be in path.
%
% Prior to running, update a few things in the overview.html file in:
%   toolbox/external/m2html/templates/frame-vincent/overview.html
%   1) The version / date
%   2) Link to rar/zip file
% Might be useful to use:
%
% USAGE
%  toolboxGenDoc
%
% INPUTS
%
% OUTPUTS
%
% EXAMPLE
%
% See also
%
% Vincent's Structure From Motion Toolbox      Version 3.0
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

% run m2html
params = {'mfiles',{'@Animation','csfm','linearAlgebra','nrsfm','sfm','visualize'}};
params = {params{:},'htmldir','doc','recursive','on','source','off'};
matlabpathTmp = matlabpath;
tmp = which('m2html','-all');
for i = 1 : length(tmp); rmpath( fileparts( tmp{i} ) ); end
addpath( genpath( fileparts( mfilename('fullpath') ) ) );
params = {params{:},'template','frame-vincent','index','menu','global','on'};
m2html(params{:});
addpath( matlabpathTmp );

% copy custom menu.html
copyfile('external\m2html\templates\menu-for-frame-vincent.html','doc/menu.html');

% copy history file
copyfile('external\history.txt','doc/history.txt');

% remove links to .svn and private in the menu.html files
d = { './doc' };

while ~isempty(d)
  dTmp = dir(d{1});
  for i = 1 : length(dTmp)
    name = dTmp(i).name;
    if strcmp( name,'.') || strcmp( name,'..'); continue; end
    if dTmp(i).isdir; d{end+1} = [ d{1} '/' name ]; continue; end
    if ~strcmp( name,'menu.html'); continue; end
    fid = fopen( [ d{1} '/' name ], 'r' ); c = fread(fid, '*char')'; fclose( fid );
    c = regexprep( c, '
([^<]*[<]?[^<]*)\.svn([^<]*[<]?[^<]*)
', '');
    c = regexprep( c, '
([^<]*[<]?[^<]*)private([^<]*[<]?[^<]*)
', '');
    fid = fopen( [ d{1} '/' name ], 'w' ); fwrite( fid, c ); fclose( fid );
  end
  d(1) = [];
end

% Zip the whole toolbox
s = dir( './' );
goodFolder = { '@Animation' 'csfm' 'doc' 'external' 'linearAlgebra' 'nrsfm' 'sfm' 'visualize' };
for i = length(s) : -1 : 1
  isGood = false;
  for j = 1 : length(goodFolder)
    if strcmp( s(i).name, goodFolder{j} )
      isGood = true;
    end
  end
  if isGood
    s(i).name = [ './' s(i).name ];
  else s(i) = [];
  end
end
files = cell(1,0);

while ~isempty(s)
  if s(1).isdir
    ls = length(s);
    s = [ s; dir( s(1).name ) ];
    for i = length(s) : -1 : ls+1
      if ~strcmp( s(i).name(1), '.' )
        s(i).name = [ s(1).name '/' s(i).name ];
      else s(i) = [];
      end
    end
    s(1) = [];
  else
    if ~strcmp( s(1).name(end), '~' )
      files{end+1} = s(1).name;
    end
    s(1) = [];
  end
end

zip( 'vincentToolbox.zip', files );

author: vrabaud


text:
%  Updates the log of the scaling parameter for mcmc algorithms
%
%  ::
%
%    log_c=update_scaling(log_c,accept_ratio,alpha_range,fixed_scaling,n,xi3)
%    log_c=update_scaling(log_c,accept_ratio,alpha_range,fixed_scaling,n,xi3,c3)
%    log_c=update_scaling(log_c,accept_ratio,alpha_range,fixed_scaling,n,xi3,c3,c_range)
%
%  Args:
%
%     - **log_c** [numeric]: initial scaling
%     - **accept_ratio** [numeric]: acceptance rate
%     - **alpha_range** [interval|{[]}]: target acceptance range
%     - **fixed_scaling** [true|false]: Do not update the scale of the algorithm
%     - **n** [integer]: iteration
%     - **xi3** [numeric]: appears in gam3 = c3*(n+1)^(-xi3). must lie in
%       (0.5,1)
%     - **c3** [numeric|{1}]: appears in gam3 = c3*(n+1)^(-xi3);
%     - **c_range** [interval|{[sqrt(eps),100]}]: range of variation of the
%       scaling parameter (not its log!!!)
%
%  Returns:
%     :
%
%     - **log_c** [numeric]: updated scaling
%
%  Note:
%     Adapts formulae 20 in Blazej Miasojedow, Eric Moulines and Matti
%     Vihola (2012): "Adaptive Parallel Tempering Algorithm"
%
%

author: jmaih


text:
function [ih,xa,ya]=imagehomog(im,h,m,clip)
%IMAGEHOMOG Apply a homography transformation to an image with bilinear interpolation
%Inputs: im(ny,nx,nc)  input image (uint8)
%        h(3,3)        homography
%        m             mode string
%                         i  show image [default if no output arguments]
%                         m  matlab coordinates (1,1) is top left [default]
%                         c  central coordinates (0,0) is the centre = {(1+nx)/2,(1+ny)/2} in 'm'
%                         k  clip to original image dimensions
%                         x  extend zero-fill to the specified clipping rectangle
%        clip(4)       bounding box [xmin xmax ymin ymax]
% Outputs:
%        ih(my,mx,nc)  output image (uint8)
%        xa(mx)        x axis
%        ya(my)        y axis

% Bugs/Suggestions:
% (1) cope with non-uint8
% (2) cope with (a) multiple inputs, (b) multiple transformations
% (3) output a boundary mask as an alpha channel
% (4) do anti-aliasing along the boundary
% (5) check that origin shift is correct for central coordinates

%      Copyright (C) Mike Brookes 2010
%      Version: $Id: imagehomog.m,v 1.3 2010/05/10 15:07:59 dmb Exp $
%
%   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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

maxby=1e7;  % maximum memory to use
[ny,nx,nc]=size(im);
if nargin<4
    clip=[];
    if nargin<3
        m='';
        if nargin<2
            h=eye(3);
        end
    end
end
imr=reshape(im,nx*ny,nc);
t=eye(3);
if any(m=='c')   % convert homography and clipping box to matlab coordinates
    t(7:8)=0.5+[nx ny]/2;  % shift origin to image centre
    h=t*h/t;  % change homography so input and output use MATLAB coordinates
    if numel(clip)
        clip=clip+t([7 7 8 8]);  % make clipping use MATLAB coordinates as well
    end
end
box=h*[1 1 nx nx; 1 ny 1 ny; 1 1 1 1];
box=box(1:2,:)./box([3 3],:);
box=[min(box(1,:)) max(box(1,:)) min(box(2,:)) max(box(2,:))];
if any(m=='k')
    clip=[1 nx 1 ny];
end
if ~numel(clip)
    clip=box;
end
clip=clip(:)';
clip(1:2:3)=floor(clip(1:2:3));
clip(2:2:4)=ceil(clip(2:2:4));
box(1:2:3)=floor(max(clip(1:2:3),box(1:2:3)));  % no point in calculating non-existant points
box(2:2:4)=ceil(min(clip(2:2:4),box(2:2:4)));
g=inv(h);
mx=box(2)-box(1)+1; % number of columns in destination
my=box(4)-box(3)+1; % number of rows in destination

ih=zeros(my*mx,nc,'uint8');
ncol=max(1,min(mx,floor(maxby/(my*nc*8)))); % number of columns to do in a chunk
nloop=ceil(mx/ncol);
cmax=1+rem(mx-1,ncol); % final column of first iteration
jxinc=ncol*my;      % increment target indices each loop
ginc=g(:,1)*ncol; % increment transformed targets each loop
wj=ones(1,jxinc);   % repeat index for ginc
jx=1+cmax*my-jxinc:cmax*my;  % initial target indices (some might be negative)
gk=g*[reshape(repmat(cmax-ncol+box(1):cmax+box(1)-1,my,1),1,jxinc); repmat(box(3):box(4),1,ncol); ones(1,jxinc)];
gn=gk(1:2,:)./gk([3 3],:);   % normalize source coordinates
mn=[zeros(1,jxinc-cmax*my) ones(1,cmax*my)]; % mask for initial iteration
mn=mn & (gn(1,:)>-0.5 & gn(2,:)>-0.5 & gn(1,:)<nx+0.5 & gn(2,:)<ny+0.5); % mask valid pixels
w3=ones(nc,1);
for i=1:nloop
    fn=max(floor(gn(:,mn)),1);
    fn1=min(max(fn(1,:)',1),nx-1);
    fn2=min(max(fn(2,:)',1),ny-1);
    dn=gn(:,mn)-[fn1 fn2]';
    dn1=min(max(dn(1,:)',0),1);
    dn2=min(max(dn(2,:)',0),1);
    dn1c=1-dn1;
    dn2c=1-dn2;
    ih(jx(mn),:)=uint8(dn1c(:,w3).*(dn2c(:,w3).*single(imr(fn2+ny*(fn1-1),:)) ...
        +dn2(:,w3).*single(imr(fn2+1+ny*(fn1-1),:))) ...
        +dn1(:,w3).*(dn2c(:,w3).*single(imr(fn2+ny*(fn1),:)) ...
        +dn2(:,w3).*single(imr(fn2+1+ny*(fn1),:))));
    jx=jx+jxinc;  % target indices
    gk=gk+ginc(:,wj);
    gn=gk(1:2,:)./gk([3 3],:);   % normalize source coordinates
    mn=gn(1,:)>-0.5 & gn(2,:)>-0.5 & gn(1,:)<nx+0.5 & gn(2,:)<ny+0.5; % mask valid pixels
end
ih=reshape(ih,[my,mx,nc]);
if any(m=='x')          % extend blank area to specified clipping rectangle
    ih=[zeros(box(3)-clip(3),clip(2)-clip(1)+1,nc,'uint8'); ...
        zeros(my,box(1)-clip(1),nc,'uint8') ih zeros(my,clip(2)-box(2),nc,'uint8');
        zeros(clip(4)-box(4),clip(2)-clip(1)+1,nc,'uint8')];
    xa=(clip(1):clip(2))-t(7);
    ya=(clip(3):clip(4))-t(8);
else
    xa=(box(1):box(2))-t(7);
    ya=(box(3):box(4))-t(8);
end

if ~nargout || any(m=='i')
    imagesc(xa,ya,ih);
    axis equal
end

author: decouples


text:
%% Example of using KST class for interfacing with KUKA iiwa robots

% An example script, it is used to show how to use the
% non-blocking moiton functions of the KST

% First start the server on the KUKA iiwa controller
% Then run the following script in Matlab

% Note you have 60 seconds to connect to server after starting the
% application (MatlabToolboxServer) from the smart pad on the robot controller.

% This example works with Sunrise application version KST_1.7  and higher.
% Copyright Mohammad SAFEEA, 26th-June-2018

close all;clear;clc;
warning('off')
%% Create the robot object
ip='172.31.1.147'; % The IP of the controller
arg1=KST.LBR7R800; % choose the robot iiwa7R800 or iiwa14R820
arg2=KST.Medien_Flansch_elektrisch; % choose the type of flange
Tef_flange=eye(4); % transofrm matrix of EEF with respect to flange
iiwa=KST(ip,arg1,arg2,Tef_flange); % create the object

%% Start a connection with the server
flag=iiwa.net_establishConnection();
if flag==0
    return;
end
pause(1);
disp('Doing some stuff')

%% Non blocking motion in Cartesian space
for i=1:3
% Go to some initial position
relVel=0.35; % over ride relative joint velocities
jPos={0, -pi / 180 * 10, 0, -pi / 180 * 100, pi / 180 * 90,pi / 180 * 90, 0};   % initial cofiguration
iiwa.movePTPJointSpace(jPos, relVel);

% Get current position of EEF
[ Pos ] = iiwa.getEEFPos(  );
distPos=Pos;
distPos{3}=distPos{3}-100;

% Move non-blocking
vel=25; %25 mm/sec
iiwa.nonBlocking_movePTPLineEEF(distPos, vel);
flag=true;
while ~iiwa.nonBlocking_isGoalReached()
    disp('Cartesian non-blocking motion in progress');
    pause(0.2);
    if(flag)
        iiwa.setBlueOn(  );
        flag=false;
    else
        iiwa.setBlueOff(  );
        flag=true;
    end
end
iiwa.setBlueOff(  );
end
%% Non blocking motion in Cartesian space
for i=1:3
% Go to some initial position
relVel=0.35; % over ride relative joint velocities
jPos={0, -pi / 180 * 10, 0, -pi / 180 * 100, pi / 180 * 90,pi / 180 * 90, 0};   % initial cofiguration
iiwa.movePTPJointSpace(jPos, relVel);

% Distination angular position
distPos=jPos;
distPos{3}=distPos{3}-pi/3;

% Move non-blocking
relVel=0.25;
iiwa.nonBlocking_movePTPJointSpace(distPos, relVel);
flag=true;
while ~iiwa.nonBlocking_isGoalReached()
    disp('Joints non-blocking motion in progress');
    pause(0.2);
    if(flag)
        iiwa.setBlueOn(  );
        flag=false;
    else
        iiwa.setBlueOff(  );
        flag=true;
    end
end
iiwa.setBlueOff(  );
end
%% Turn off server
disp('Distination reached');
disp('Turning off server');
iiwa.net_turnOffServer(  );

repo: KST-Kuka-Sunrise-Toolbox