text:
%This m file is for locating a good enough "M". For a good value of "M",
%"PG" will be high. But it also adds computational load as we increase "M".
%Run this for different frames of "x" and see for yourself. Define the
%starting values of frames by defining a value to "b" at line 25.

%Here,  array of "k"=Reflection C-oefficients
%       array of "a"=LPCs
%       array of "R"=autocorrelation co-efficients
%every notation is according to page 112 [chapter4] of Speech Coding
%Algorithms by W. C. Chu.

%the indexes in this code is "+1" from the indexes in the book in some
%cases because there is no index "0" in matlab

clear all
% close all
clc

for M=3:100,

% INITIALIZATION:
inpfilenm = 's1ofwb.wav';
[x, fs]=wavread(inpfilenm);%"t_16k_2s.wav" is the file I used. Change according
                      %to your case.
% length(x)
b=3841;        %index no. of starting data point of current frame
fsize = 30e-3;    %frame size
frame_length = round(fs .* fsize);  %=number of data points in each framesize
                                    %of "x"
N= frame_length - 1;        %N+1 = frame length = number of data points in
                            %each frame
sk=0;       %initializing summartion term "sk"
a=[zeros(M+1);zeros(M+1)]; %defining a matrix of zeros for "a" for init.

%FRAME SEGMENTATION:
    y1=x(b:b+N);    %"b+N" denotes the end point of current frame.
                    %"y" denotes an array of the data points of the current
                    %frame
    y = filter([1 -.9378], 1, y1);  %pre-emphasis filtering

%MAIN BODY OF THIS PROGRAM STARTS FROM HERE>>>>>>>>>>>>>>
z=xcorr(y);

%finding array of R[l]
R=z( ( (length(z)+1) ./2 ) : length(z)); %R=array of "R[l]", where l=0,1,2,
                                         %...(b+N)-1
                                         %R(1)=R[lag=0], R(2)=R[lag=1],
                                         %R(3)=R[lag=2]... etc

%GETTING OTHER PARAMETERS OF PREDICTOR OF ORDER "0":
s=1;        %s=step no.
J(1)=R(1);          %J=array of "Jl", where l=0,1,2...(b+N)-1
                    %J(1)=J0, J(2)=J1, J(3)=J2 etc

%GETTING OTHER PARAMETERS OF PREDICTOR OF ORDER "(s-1)":
for s=2:M+1,
    sk=0;               %clearing "sk" for each iteration
    for i=2:(s-1),
        sk=sk + a(i,(s-1)).*R(s-i+1);
    end                 %now we know value of "sk", the summation term
                        %of formula of calculating "k(l)"
    k(s)=(R(s) + sk)./J(s-1);
    J(s)=J(s-1).*(1-(k(s)).^2);

    a(s,s)= -k(s);
    a(1,s)=1;
    for i=2:(s-1),
        a(i,s)=a(i,(s-1)) - k(s).*a((s-i+1),(s-1));
    end
end

a_final=a((1:s),s)';
est_y = filter([0 -a_final(2:end)],1,y);    % = s^(n) with a cap on page 92 of the book
e = y - est_y;      %supposed to be a white noise

pg(M) = 10.*log10( (sum(y.^2)) ./ (sum(e.^2)) );       %prediction gain
end
figure;
subplot(2,1,1), plot(x); title(['original speech file, ',inpfilenm]);
subplot(2,1,2),plot(pg);
title('Prediction Gain (PG) vs Prediction Order (M) for frame starting at data point "b"');
xlabel('M');
ylabel('PG');

disp('This m file is simply "func_lev_durb.m" with a little modification. If you run this for many frames, you will see that PG is good enough for M = a value around "10" for most of the frames');


text:
function b = generative_model(A,D,m,modeltype,modelvar,params,epsilon)
% GENERATIVE_MODEL          Run generative model code
%
%   B = GENERATIVE_MODEL(A,D,m,modeltype,modelvar,params)
%
%   Generates synthetic networks using the models described in the study by
%   Betzel et al (2016) in Neuroimage.
%
%   Inputs:
%           A,          binary network of seed connections
%           D,          Euclidean distance/fiber length matrix
%           m,          number of connections that should be present in
%                       final synthetic network
%           modeltype,  specifies the generative rule (see below)
%           modelvar,   specifies whether the generative rules are based on
%                       power-law or exponential relationship
%                       ({'powerlaw'}|{'exponential})
%           params,     either a vector (in the case of the geometric
%                       model) or a matrix (for all other models) of
%                       parameters at which the model should be evaluated.
%           epsilon,    the baseline probability of forming a particular
%                       connection (should be a very small number
%                       {default = 1e-5}).
%
%   Output:
%           B,          m x number of networks matrix of connections
%
%
%   Full list of model types:
%   (each model type realizes a different generative rule)
%
%       1.  'sptl'          spatial model
%       2.  'neighbors'     number of common neighbors
%       3.  'matching'      matching index
%       4.  'clu-avg'       average clustering coeff.
%       5.  'clu-min'       minimum clustering coeff.
%       6.  'clu-max'       maximum clustering coeff.
%       7.  'clu-diff'      difference in clustering coeff.
%       8.  'clu-prod'      product of clustering coeff.
%       9.  'deg-avg'       average degree
%       10. 'deg-min'       minimum degree
%       11. 'deg-max'       maximum degree
%       12. 'deg-diff'      difference in degree
%       13. 'deg-prod'      product of degree
%
%
%   Example usage:
%
%       load demo_generative_models_data
%
%       % get number of bi-directional connections
%       m = nnz(A)/2;
%
%       % get cardinality of network
%       n = length(A);
%
%       % set model type
%       modeltype = 'neighbors';
%
%       % set whether the model is based on powerlaw or exponentials
%       modelvar = [{'powerlaw'},{'powerlaw'}];
%
%       % choose some model parameters
%       params = [-2,0.2; -5,1.2; -1,1.5];
%       nparams = size(params,1);
%
%       % generate synthetic networks
%       B = generative_model(Aseed,D,m,modeltype,modelvar,params);
%
%       % store them in adjacency matrix format
%       Asynth = zeros(n,n,nparams);
%       for i = 1:nparams;
%           a = zeros(n); a(B(:,i)) = 1; a = a + a';
%           Asynth(:,:,i) = a;
%       end
%
%   Reference: Betzel et al (2016) Neuroimage 124:1054-64.
%
%   Richard Betzel, Indiana University/University of Pennsylvania, 2015

if ~exist('epsilon','var')
    epsilon = 1e-5;
end

n = length(D);
nparams = size(params,1);
b = zeros(m,nparams);

switch modeltype

    case 'clu-avg'
        clu = clustering_coef_bu(A);
        Kseed = bsxfun(@plus,clu(:,ones(1,n)),clu')/2;
        for iparam = 1:nparams
            eta = params(iparam,1);
            gam = params(iparam,2);
            b(:,iparam) = fcn_clu_avg(A,Kseed,D,m,eta,gam,modelvar,epsilon);
        end

    case 'clu-diff'
        clu = clustering_coef_bu(A);
        Kseed = abs(bsxfun(@minus,clu(:,ones(1,n)),clu'));
        for iparam = 1:nparams
            eta = params(iparam,1);
            gam = params(iparam,2);
            b(:,iparam) = fcn_clu_diff(A,Kseed,D,m,eta,gam,modelvar,epsilon);
        end

    case 'clu-max'
        clu = clustering_coef_bu(A);
        Kseed = bsxfun(@max,clu(:,ones(1,n)),clu');
        for iparam = 1:nparams
            eta = params(iparam,1);
            gam = params(iparam,2);
            b(:,iparam) = fcn_clu_max(A,Kseed,D,m,eta,gam,modelvar,epsilon);
        end

    case 'clu-min'
        clu = clustering_coef_bu(A);
        Kseed = bsxfun(@min,clu(:,ones(1,n)),clu');
        for iparam = 1:nparams
            eta = params(iparam,1);
            gam = params(iparam,2);
            b(:,iparam) = fcn_clu_min(A,Kseed,D,m,eta,gam,modelvar,epsilon);
        end

    case 'clu-prod'
        clu = clustering_coef_bu(A);
        Kseed = clu*clu';
        for iparam = 1:nparams
            eta = params(iparam,1);
            gam = params(iparam,2);
            b(:,iparam) = fcn_clu_prod(A,Kseed,D,m,eta,gam,modelvar,epsilon);
        end

    case 'deg-avg'
        kseed = sum(A,2);
        Kseed = bsxfun(@plus,kseed(:,ones(1,n)),kseed')/2;
        for iparam = 1:nparams
            eta = params(iparam,1);
            gam = params(iparam,2);
            b(:,iparam) = fcn_deg_avg(A,Kseed,D,m,eta,gam,modelvar,epsilon);
        end

    case 'deg-diff'
        kseed = sum(A,2);
        Kseed = abs(bsxfun(@minus,kseed(:,ones(1,n)),kseed'));
        for iparam = 1:nparams
            eta = params(iparam,1);
            gam = params(iparam,2);
            b(:,iparam) = fcn_deg_diff(A,Kseed,D,m,eta,gam,modelvar,epsilon);
        end

    case 'deg-max'
        kseed = sum(A,2);
        Kseed = bsxfun(@max,kseed(:,ones(1,n)),kseed');
        for iparam = 1:nparams
            eta = params(iparam,1);
            gam = params(iparam,2);
            b(:,iparam) = fcn_deg_max(A,Kseed,D,m,eta,gam,modelvar,epsilon);
        end

    case 'deg-min'
        kseed = sum(A,2);
        Kseed = bsxfun(@min,kseed(:,ones(1,n)),kseed');
        for iparam = 1:nparams
            eta = params(iparam,1);
            gam = params(iparam,2);
            b(:,iparam) = fcn_deg_min(A,Kseed,D,m,eta,gam,modelvar,epsilon);
        end

    case 'deg-prod'
        kseed = sum(A,2);
        Kseed = (kseed*kseed').*~eye(n);
        for iparam = 1:nparams
            eta = params(iparam,1);
            gam = params(iparam,2);
            b(:,iparam) = fcn_deg_prod(A,Kseed,D,m,eta,gam,modelvar,epsilon);
        end

    case 'neighbors'
        Kseed = (A*A).*~eye(n);
        for iparam = 1:nparams
            eta = params(iparam,1);
            gam = params(iparam,2);
            b(:,iparam) = fcn_nghbrs(A,Kseed,D,m,eta,gam,modelvar,epsilon);
        end

    case 'matching'
        Kseed = matching_ind(A);
        Kseed = Kseed + Kseed';
        for iparam = 1:nparams
            eta = params(iparam,1);
            gam = params(iparam,2);
            b(:,iparam) = fcn_matching(A,Kseed,D,m,eta,gam,modelvar,epsilon);
        end

    case 'sptl'
        for iparam = 1:nparams
            eta = params(iparam,1);
            b(:,iparam) = fcn_sptl(A,D,m,eta,modelvar{1});
        end

end

function b = fcn_clu_avg(A,K,D,m,eta,gam,modelvar,epsilon)
K = K + epsilon;
n = length(D);
mseed = nnz(A)/2;
A = A > 0;
mv1 = modelvar{1};
mv2 = modelvar{2};
switch mv1
    case 'powerlaw'
        Fd = D.^eta;
    case 'exponential'
        Fd = exp(eta*D);
end
switch mv2
    case 'powerlaw'
        Fk = K.^gam;
    case 'exponential'
        Fk = exp(gam*K);
end

c = clustering_coef_bu(A);
k = sum(A,2);

Ff = Fd.*Fk.*~A;
[u,v] = find(triu(ones(n),1));
indx = (v - 1)*n + u;
P = Ff(indx);

for i = (mseed + 1):m
    C = [0; cumsum(P)];
    r = sum(rand*C(end) >= C);
    uu = u(r);
    vv = v(r);
    A(uu,vv) = 1;
    A(vv,uu) = 1;
    k([uu,vv]) = k([uu,vv]) + 1;
    bu = A(uu,:);
    su = A(bu,bu);
    bv = A(vv,:);
    sv = A(bv,bv);
    bth = bu & bv;
    c(bth) = c(bth) + 2./(k(bth).^2 - k(bth));
    c(uu) = nnz(su)/(k(uu)*(k(uu) - 1));
    c(vv) = nnz(sv)/(k(vv)*(k(vv) - 1));
    c(k <= 1) = 0;
    bth([uu,vv]) = true;
    K(:,bth) = bsxfun(@plus,c(:,ones(1,sum(bth))),c(bth,:)')/2 + epsilon;
    K(bth,:) = bsxfun(@plus,c(:,ones(1,sum(bth))),c(bth,:)')'/2 + epsilon;

    switch mv2
        case 'powerlaw'
            Ff(bth,:) = Fd(bth,:).*((K(bth,:)).^gam);
            Ff(:,bth) = Fd(:,bth).*((K(:,bth)).^gam);
        case 'exponential'
            Ff(bth,:) = Fd(bth,:).*exp((K(bth,:))*gam);
            Ff(:,bth) = Fd(:,bth).*exp((K(:,bth))*gam);
    end
    Ff = Ff.*~A;
    P = Ff(indx);
end
b = find(triu(A,1));

function b = fcn_clu_diff(A,K,D,m,eta,gam,modelvar,epsilon)
K = K + epsilon;
n = length(D);
mseed = nnz(A)/2;
A = A > 0;
mv1 = modelvar{1};
mv2 = modelvar{2};
switch mv1
    case 'powerlaw'
        Fd = D.^eta;
    case 'exponential'
        Fd = exp(eta*D);
end
switch mv2
    case 'powerlaw'
        Fk = K.^gam;
    case 'exponential'
        Fk = exp(gam*K);
end

c = clustering_coef_bu(A);
k = sum(A,2);

Ff = Fd.*Fk.*~A;
[u,v] = find(triu(ones(n),1));
indx = (v - 1)*n + u;
P = Ff(indx);

for i = (mseed + 1):m
    C = [0; cumsum(P)];
    r = sum(rand*C(end) >= C);
    uu = u(r);
    vv = v(r);
    A(uu,vv) = 1;
    A(vv,uu) = 1;
    k([uu,vv]) = k([uu,vv]) + 1;
    bu = A(uu,:);
    su = A(bu,bu);
    bv = A(vv,:);
    sv = A(bv,bv);
    bth = bu & bv;
    c(bth) = c(bth) + 2./(k(bth).^2 - k(bth));
    c(uu) = nnz(su)/(k(uu)*(k(uu) - 1));
    c(vv) = nnz(sv)/(k(vv)*(k(vv) - 1));
    c(k <= 1) = 0;
    bth([uu,vv]) = true;
    K(:,bth) = abs(bsxfun(@minus,c(:,ones(1,sum(bth))),c(bth,:)')) + epsilon;
    K(bth,:) = abs(bsxfun(@minus,c(:,ones(1,sum(bth))),c(bth,:)'))' + epsilon;

    switch mv2
        case 'powerlaw'
            Ff(bth,:) = Fd(bth,:).*((K(bth,:)).^gam);
            Ff(:,bth) = Fd(:,bth).*((K(:,bth)).^gam);
        case 'exponential'
            Ff(bth,:) = Fd(bth,:).*exp((K(bth,:))*gam);
            Ff(:,bth) = Fd(:,bth).*exp((K(:,bth))*gam);
    end
    Ff = Ff.*~A;
    P = Ff(indx);
end
b = find(triu(A,1));

function b = fcn_clu_max(A,K,D,m,eta,gam,modelvar,epsilon)
K = K + epsilon;
n = length(D);
mseed = nnz(A)/2;
A = A > 0;
mv1 = modelvar{1};
mv2 = modelvar{2};
switch mv1
    case 'powerlaw'
        Fd = D.^eta;
    case 'exponential'
        Fd = exp(eta*D);
end
switch mv2
    case 'powerlaw'
        Fk = K.^gam;
    case 'exponential'
        Fk = exp(gam*K);
end

c = clustering_coef_bu(A);
k = sum(A,2);

Ff = Fd.*Fk.*~A;
[u,v] = find(triu(ones(n),1));
indx = (v - 1)*n + u;
P = Ff(indx);

for i = (mseed + 1):m
    C = [0; cumsum(P)];
    r = sum(rand*C(end) >= C);
    uu = u(r);
    vv = v(r);
    A(uu,vv) = 1;
    A(vv,uu) = 1;
    k([uu,vv]) = k([uu,vv]) + 1;
    bu = A(uu,:);
    su = A(bu,bu);
    bv = A(vv,:);
    sv = A(bv,bv);
    bth = bu & bv;
    c(bth) = c(bth) + 2./(k(bth).^2 - k(bth));
    c(uu) = nnz(su)/(k(uu)*(k(uu) - 1));
    c(vv) = nnz(sv)/(k(vv)*(k(vv) - 1));
    c(k <= 1) = 0;
    bth([uu,vv]) = true;
    K(:,bth) = bsxfun(@max,c(:,ones(1,sum(bth))),c(bth,:)') + epsilon;
    K(bth,:) = bsxfun(@max,c(:,ones(1,sum(bth))),c(bth,:)')' + epsilon;

    switch mv2
        case 'powerlaw'
            Ff(bth,:) = Fd(bth,:).*((K(bth,:)).^gam);
            Ff(:,bth) = Fd(:,bth).*((K(:,bth)).^gam);
        case 'exponential'
            Ff(bth,:) = Fd(bth,:).*exp((K(bth,:))*gam);
            Ff(:,bth) = Fd(:,bth).*exp((K(:,bth))*gam);
    end
    Ff = Ff.*~A;
    P = Ff(indx);
end
b = find(triu(A,1));

function b = fcn_clu_min(A,K,D,m,eta,gam,modelvar,epsilon)
K = K + epsilon;
n = length(D);
mseed = nnz(A)/2;
A = A > 0;
mv1 = modelvar{1};
mv2 = modelvar{2};
switch mv1
    case 'powerlaw'
        Fd = D.^eta;
    case 'exponential'
        Fd = exp(eta*D);
end
switch mv2
    case 'powerlaw'
        Fk = K.^gam;
    case 'exponential'
        Fk = exp(gam*K);
end

c = clustering_coef_bu(A);
k = sum(A,2);

Ff = Fd.*Fk.*~A;
[u,v] = find(triu(ones(n),1));
indx = (v - 1)*n + u;
P = Ff(indx);

for i = (mseed + 1):m
    C = [0; cumsum(P)];
    r = sum(rand*C(end) >= C);
    uu = u(r);
    vv = v(r);
    A(uu,vv) = 1;
    A(vv,uu) = 1;
    k([uu,vv]) = k([uu,vv]) + 1;
    bu = A(uu,:);
    su = A(bu,bu);
    bv = A(vv,:);
    sv = A(bv,bv);
    bth = bu & bv;
    c(bth) = c(bth) + 2./(k(bth).^2 - k(bth));
    c(uu) = nnz(su)/(k(uu)*(k(uu) - 1));
    c(vv) = nnz(sv)/(k(vv)*(k(vv) - 1));
    c(k <= 1) = 0;
    bth([uu,vv]) = true;
    K(:,bth) = bsxfun(@min,c(:,ones(1,sum(bth))),c(bth,:)') + epsilon;
    K(bth,:) = bsxfun(@min,c(:,ones(1,sum(bth))),c(bth,:)')' + epsilon;

    switch mv2
        case 'powerlaw'
            Ff(bth,:) = Fd(bth,:).*((K(bth,:)).^gam);
            Ff(:,bth) = Fd(:,bth).*((K(:,bth)).^gam);
        case 'exponential'
            Ff(bth,:) = Fd(bth,:).*exp((K(bth,:))*gam);
            Ff(:,bth) = Fd(:,bth).*exp((K(:,bth))*gam);
    end
    Ff = Ff.*~A;
    P = Ff(indx);
end
b = find(triu(A,1));

function b = fcn_clu_prod(A,K,D,m,eta,gam,modelvar,epsilon)
K = K + epsilon;
n = length(D);
mseed = nnz(A)/2;
A = A > 0;
mv1 = modelvar{1};
mv2 = modelvar{2};
switch mv1
    case 'powerlaw'
        Fd = D.^eta;
    case 'exponential'
        Fd = exp(eta*D);
end
switch mv2
    case 'powerlaw'
        Fk = K.^gam;
    case 'exponential'
        Fk = exp(gam*K);
end

c = clustering_coef_bu(A);
k = sum(A,2);

Ff = Fd.*Fk.*~A;
[u,v] = find(triu(ones(n),1));
indx = (v - 1)*n + u;
P = Ff(indx);

for i = (mseed + 1):m
    C = [0; cumsum(P)];
    r = sum(rand*C(end) >= C);
    uu = u(r);
    vv = v(r);
    A(uu,vv) = 1;
    A(vv,uu) = 1;
    k([uu,vv]) = k([uu,vv]) + 1;
    bu = A(uu,:);
    su = A(bu,bu);
    bv = A(vv,:);
    sv = A(bv,bv);
    bth = bu & bv;
    c(bth) = c(bth) + 2./(k(bth).^2 - k(bth));
    c(uu) = nnz(su)/(k(uu)*(k(uu) - 1));
    c(vv) = nnz(sv)/(k(vv)*(k(vv) - 1));
    c(k <= 1) = 0;
    bth([uu,vv]) = true;
    K(bth,:) = (c(bth,:)*c') + epsilon;
    K(:,bth) = (c*c(bth,:)') + epsilon;

    switch mv2
        case 'powerlaw'
            Ff(bth,:) = Fd(bth,:).*((K(bth,:)).^gam);
            Ff(:,bth) = Fd(:,bth).*((K(:,bth)).^gam);
        case 'exponential'
            Ff(bth,:) = Fd(bth,:).*exp((K(bth,:))*gam);
            Ff(:,bth) = Fd(:,bth).*exp((K(:,bth))*gam);
    end
    Ff = Ff.*~A;
    P = Ff(indx);
end
b = find(triu(A,1));

function b = fcn_deg_avg(A,K,D,m,eta,gam,modelvar,epsilon)
n = length(D);
mseed = nnz(A)/2;
k = sum(A,2);
[u,v] = find(triu(ones(n),1));
indx = (v - 1)*n + u;
D = D(indx);
mv1 = modelvar{1};
mv2 = modelvar{2};
switch mv1
    case 'powerlaw'
        Fd = D.^eta;
    case 'exponential'
        Fd = exp(eta*D);
end
K = K + epsilon;
switch mv2
    case 'powerlaw'
        Fk = K.^gam;
    case 'exponential'
        Fk = exp(gam*K);
end
P = Fd.*Fk(indx).*~A(indx);
b = zeros(m,1);
b(1:mseed) = find(A(indx));
for i = (mseed + 1):m
    C = [0; cumsum(P)];
    r = sum(rand*C(end) >= C);
    w = [u(r),v(r)];
    k(w) = k(w) + 1;
    switch mv2
        case 'powerlaw'
            Fk(:,w) = [((k + k(w(1)))/2) + epsilon, ((k + k(w(2)))/2) + epsilon].^gam;
            Fk(w,:) = ([((k + k(w(1)))/2) + epsilon, ((k + k(w(2)))/2) + epsilon].^gam)';
        case 'exponential'
            Fk(:,w) = exp([((k + k(w(1)))/2) + epsilon, ((k + k(w(2)))/2) + epsilon]*gam);
            Fk(w,:) = exp([((k + k(w(1)))/2) + epsilon, ((k + k(w(2)))/2) + epsilon]*gam)';
    end
    P = Fd.*Fk(indx);
    b(i) = r;
    P(b(1:i)) = 0;
end
b = indx(b);

function b = fcn_deg_diff(A,K,D,m,eta,gam,modelvar,epsilon)
n = length(D);
mseed = nnz(A)/2;
k = sum(A,2);
[u,v] = find(triu(ones(n),1));
indx = (v - 1)*n + u;
D = D(indx);
mv1 = modelvar{1};
mv2 = modelvar{2};
switch mv1
    case 'powerlaw'
        Fd = D.^eta;
    case 'exponential'
        Fd = exp(eta*D);
end
K = K + epsilon;
switch mv2
    case 'powerlaw'
        Fk = K.^gam;
    case 'exponential'
        Fk = exp(gam*K);
end
P = Fd.*Fk(indx).*~A(indx);
b = zeros(m,1);
b(1:mseed) = find(A(indx));
for i = (mseed + 1):m
    C = [0; cumsum(P)];
    r = sum(rand*C(end) >= C);

    w = [u(r),v(r)];
    k(w) = k(w) + 1;
    switch mv2
        case 'powerlaw'
            Fk(:,w) = (abs([k - k(w(1)), k - k(w(2))]) + epsilon).^gam;
            Fk(w,:) = ((abs([k - k(w(1)), k - k(w(2))]) + epsilon).^gam)';
        case 'exponential'
            Fk(:,w) = exp((abs([k - k(w(1)), k - k(w(2))]) + epsilon)*gam);
            Fk(w,:) = exp((abs([k - k(w(1)), k - k(w(2))]) + epsilon)*gam)';
    end
    P = Fd.*Fk(indx);
    b(i) = r;
    P(b(1:i)) = 0;
end
b = indx(b);

function b = fcn_deg_min(A,K,D,m,eta,gam,modelvar,epsilon)
n = length(D);
mseed = nnz(A)/2;
k = sum(A,2);
[u,v] = find(triu(ones(n),1));
indx = (v - 1)*n + u;
D = D(indx);
mv1 = modelvar{1};
mv2 = modelvar{2};
switch mv1
    case 'powerlaw'
        Fd = D.^eta;
    case 'exponential'
        Fd = exp(eta*D);
end
K = K + epsilon;
switch mv2
    case 'powerlaw'
        Fk = K.^gam;
    case 'exponential'
        Fk = exp(gam*K);
end
P = Fd.*Fk(indx).*~A(indx);
b = zeros(m,1);
b(1:mseed) = find(A(indx));
for i = (mseed + 1):m
    C = [0; cumsum(P)];
    r = sum(rand*C(end) >= C);
    w = [u(r),v(r)];
    k(w) = k(w) + 1;
    switch mv2
        case 'powerlaw'
            Fk(:,w) = [min(k,k(w(1))) + epsilon, min(k,k(w(2))) + epsilon].^gam;
            Fk(w,:) = ([min(k,k(w(1))) + epsilon, min(k,k(w(2))) + epsilon].^gam)';
        case 'exponential'
            Fk(:,w) = exp([min(k,k(w(1))) + epsilon, min(k,k(w(2))) + epsilon]*gam);
            Fk(w,:) = exp([min(k,k(w(1))) + epsilon, min(k,k(w(2))) + epsilon]*gam)';
    end
    P = Fd.*Fk(indx);
    b(i) = r;
    P(b(1:i)) = 0;
end
b = indx(b);

function b = fcn_deg_max(A,K,D,m,eta,gam,modelvar,epsilon)
n = length(D);
mseed = nnz(A)/2;
k = sum(A,2);
[u,v] = find(triu(ones(n),1));
indx = (v - 1)*n + u;
D = D(indx);
mv1 = modelvar{1};
mv2 = modelvar{2};
switch mv1
    case 'powerlaw'
        Fd = D.^eta;
    case 'exponential'
        Fd = exp(eta*D);
end
K = K + epsilon;
switch mv2
    case 'powerlaw'
        Fk = K.^gam;
    case 'exponential'
        Fk = exp(gam*K);
end
P = Fd.*Fk(indx).*~A(indx);
b = zeros(m,1);
b(1:mseed) = find(A(indx));
for i = (mseed + 1):m
    C = [0; cumsum(P)];
    r = sum(rand*C(end) >= C);
    w = [u(r),v(r)];
    k(w) = k(w) + 1;
    switch mv2
        case 'powerlaw'
            Fk(:,w) = [max(k,k(w(1))) + epsilon, max(k,k(w(2))) + epsilon].^gam;
            Fk(w,:) = ([max(k,k(w(1))) + epsilon, max(k,k(w(2))) + epsilon].^gam)';
        case 'exponential'
            Fk(:,w) = exp([max(k,k(w(1))) + epsilon, max(k,k(w(2))) + epsilon]*gam);
            Fk(w,:) = exp([max(k,k(w(1))) + epsilon, max(k,k(w(2))) + epsilon]*gam)';
    end
    P = Fd.*Fk(indx);
    b(i) = r;
    P(b(1:i)) = 0;
end
b = indx(b);

function b = fcn_deg_prod(A,K,D,m,eta,gam,modelvar,epsilon)
n = length(D);
mseed = nnz(A)/2;
k = sum(A,2);
[u,v] = find(triu(ones(n),1));
indx = (v - 1)*n + u;
D = D(indx);
mv1 = modelvar{1};
mv2 = modelvar{2};
switch mv1
    case 'powerlaw'
        Fd = D.^eta;
    case 'exponential'
        Fd = exp(eta*D);
end
K = K + epsilon;
switch mv2
    case 'powerlaw'
        Fk = K.^gam;
    case 'exponential'
        Fk = exp(gam*K);
end
P = Fd.*Fk(indx).*~A(indx);
b = zeros(m,1);
b(1:mseed) = find(A(indx));
for i = (mseed + 1):m
    C = [0; cumsum(P)];
    r = sum(rand*C(end) >= C);
    w = [u(r),v(r)];
    k(w) = k(w) + 1;
    switch mv2
        case 'powerlaw'
            Fk(:,w) = ([k*k(w(1)) + epsilon, k*k(w(2)) + epsilon].^gam);
            Fk(w,:) = (([k*k(w(1)) + epsilon, k*k(w(2)) + epsilon].^gam)');
        case 'exponential'
            Fk(:,w) = exp([k*k(w(1)) + epsilon, k*k(w(2)) + epsilon]*gam);
            Fk(w,:) = exp([k*k(w(1)) + epsilon, k*k(w(2)) + epsilon]*gam)';
    end
    P = Fd.*Fk(indx);
    b(i) = r;
    P(b(1:i)) = 0;
end
b = indx(b);

function b = fcn_nghbrs(A,K,D,m,eta,gam,modelvar,epsilon)
K = K + epsilon;
n = length(D);
mseed = nnz(A)/2;
A = A > 0;
mv1 = modelvar{1};
mv2 = modelvar{2};
switch mv1
    case 'powerlaw'
        Fd = D.^eta;
    case 'exponential'
        Fd = exp(eta*D);
end
switch mv2
    case 'powerlaw'
%         gam = abs(gam);
        Fk = K.^gam;
    case 'exponential'
        Fk = exp(gam*K);
end
Ff = Fd.*Fk.*~A;
[u,v] = find(triu(ones(n),1));
indx = (v - 1)*n + u;
P = Ff(indx);
for i = (mseed + 1):m
    C = [0; cumsum(P)];
    r = sum(rand*C(end) >= C);
    uu = u(r);
    vv = v(r);
    x = A(uu,:);
    y = A(:,vv);
    A(uu,vv) = 1;
    A(vv,uu) = 1;
    K(uu,y) = K(uu,y) + 1;
    K(y,uu) = K(y,uu) + 1;
    K(vv,x) = K(vv,x) + 1;
    K(x,vv) = K(x,vv) + 1;
    switch mv2
        case 'powerlaw'
            Ff(uu,y) = Fd(uu,y).*(K(uu,y).^gam);
            Ff(y,uu) = Ff(uu,y)';
            Ff(vv,x) = Fd(vv,x).*(K(vv,x).^gam);
            Ff(x,vv) = Ff(vv,x)';
        case 'exponential'
            Ff(uu,y) = Fd(uu,y).*exp(K(uu,y)*gam);
            Ff(y,uu) = Ff(uu,y)';
            Ff(vv,x) = Fd(vv,x).*exp(K(vv,x)*gam);
            Ff(x,vv) = Ff(vv,x)';
    end
    Ff(A) = 0;
    P = Ff(indx);
end
b = find(triu(A,1));

function b = fcn_matching(A,K,D,m,eta,gam,modelvar,epsilon)
K = K + epsilon;
n = length(D);
mseed = nnz(A)/2;
mv1 = modelvar{1};
mv2 = modelvar{2};
switch mv1
    case 'powerlaw'
        Fd = D.^eta;
    case 'exponential'
        Fd = exp(eta*D);
end
switch mv2
    case 'powerlaw'
        Fk = K.^gam;
    case 'exponential'
        Fk = exp(gam*K);
end
Ff = Fd.*Fk.*~A;
[u,v] = find(triu(ones(n),1));
indx = (v - 1)*n + u;
P = Ff(indx);
for ii = (mseed + 1):m

    C = [0; cumsum(P)];
    r = sum(rand*C(end) >= C);
    uu = u(r);
    vv = v(r);

    A(uu,vv) = 1;
    A(vv,uu) = 1;

    updateuu = find(A*A(:,uu));
    updateuu(updateuu == uu) = [];
    updateuu(updateuu == vv) = [];

    updatevv = find(A*A(:,vv));
    updatevv(updatevv == uu) = [];
    updatevv(updatevv == vv) = [];

    c1 = [A(:,uu)', A(uu,:)];
    for i = 1:length(updateuu)
        j = updateuu(i);
        c2 = [A(:,j)' A(j,:)];
        use = ~(~c1&~c2);
        use(uu) = 0;  use(uu+n) = 0;
        use(j) = 0;  use(j+n) = 0;
        ncon = sum(c1(use))+sum(c2(use));
        if (ncon==0)
            K(uu,j) = epsilon;
            K(j,uu) = epsilon;
        else
            K(uu,j) = (2*(sum(c1(use)&c2(use))/ncon)) + epsilon;
            K(j,uu) = K(uu,j);
        end

    end

    c1 = [A(:,vv)', A(vv,:)];
    for i = 1:length(updatevv)
        j = updatevv(i);
        c2 = [A(:,j)' A(j,:)];
        use = ~(~c1&~c2);
        use(vv) = 0;  use(vv+n) = 0;
        use(j) = 0;  use(j+n) = 0;
        ncon = sum(c1(use))+sum(c2(use));
        if (ncon==0)
            K(vv,j) = epsilon;
            K(j,vv) = epsilon;
        else
            K(vv,j) = (2*(sum(c1(use)&c2(use))/ncon)) + epsilon;
            K(j,vv) = K(vv,j);
        end
    end
    switch mv2
        case 'powerlaw'
            Fk = K.^gam;
        case 'exponential'
            Fk = exp(gam*K);
    end
    Ff = Fd.*Fk.*~A;
    P = Ff(indx);
end
b = find(triu(A,1));

function b = fcn_sptl(A,D,m,eta,modelvar)
n = length(D);
mseed = nnz(A)/2;
switch modelvar
    case 'powerlaw'
        Fd = D.^eta;
    case 'exponential'
        Fd = exp(eta*D);
end
[u,v] = find(triu(ones(n),1));
indx = (v - 1)*n + u;
P = Fd(indx).*~A(indx);
b = zeros(m,1);
b(1:mseed) = find(A(indx));
for i = (mseed + 1):m
    C = [0; cumsum(P)];
    r = sum(rand*C(end) >= C);
    b(i) = r;
    P = Fd(indx);
    P(b(1:i)) = 0;
end
b = indx(b);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [varargout]=subQuad(varargin)

% function [Fs,Vs,C,CV]=subQuad(F,V,n,splitMethod)
% ------------------------------------------------------------------------
% Sub-devides the quadrilateral faces defined by the patch format data F
% (faces) and V (vertices). Each face is split n times using the specified
% split method (splitMethod). The user may request the following outputs:
% The new faces: Fs
% The new coordinates: Vs
% Face color labels: C
% Nodal labels: CV
%
% Four split methods are defined:
% 1: General linear resampling
% 2: Linear resampling only in the first direction
% 3: Linear resampling only in the second direction
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2010/06/01 Created
% 2017/11/29 Fixed single face input related bug
% 2018/11/06 Added splitMethod for splitting in a certain direction
% 2018/11/06 Added varargin based input handling
% 2018/11/06 Added color data handling
% 2018/11/06 Added variable output handling
% 2021/07/11 Fixed handling of CV
% ------------------------------------------------------------------------

%% parse input

switch nargin
    case 2
        F=varargin{1};
        V=varargin{2};
        n=1;
        splitMethod=1;
    case 3
        F=varargin{1};
        V=varargin{2};
        n=varargin{3};
        splitMethod=1;
    case 4
        F=varargin{1};
        V=varargin{2};
        n=varargin{3};
        splitMethod=varargin{4};
end

C=(1:1:size(F,1))'; %Face colors or indices

%%

if n>0
    for qIter=1:1:n
        switch splitMethod
            case 1

                edgeMat=[F(:,[1 2]); F(:,[2 3]);  F(:,[3 4]); F(:,[4 1])]; %Edges matrix
                Es=sort(edgeMat,2); %Sorted edges matrix
                [~,ind1,~]=unique(Es,'rows');
                edgeMat=edgeMat(ind1,:);

                numPoints = size(V,1);
                numEdges = size(edgeMat,1);

                % Get indices of the three edges associated with each face
                A = sparse(edgeMat(:,1),edgeMat(:,2),(1:numEdges)+numPoints,numPoints,numPoints,numEdges);
                A = max(A,A'); %Copy symmetric

                %Indices for A matrix
                indA_12=F(:,1)+(F(:,2)-1)*numPoints;
                indA_23=F(:,2)+(F(:,3)-1)*numPoints;
                indA_34=F(:,3)+(F(:,4)-1)*numPoints;
                indA_41=F(:,4)+(F(:,1)-1)*numPoints;

                %Get indices for vertex array
                indV_12=full(A(indA_12));
                indV_23=full(A(indA_23));
                indV_34=full(A(indA_34));
                indV_41=full(A(indA_41));

                indV_mid=(1:1:size(F,1))'+numPoints+size(edgeMat,1);

                %Create faces array
                Fs=[[F(:,1)  indV_12 indV_mid indV_41];...
                    [F(:,2)  indV_23 indV_mid indV_12];...
                    [F(:,3)  indV_34 indV_mid indV_23];...
                    [F(:,4)  indV_41 indV_mid indV_34]];

                %Create vertex array
                Vn=0.5*(V(edgeMat(:,1),:)+V(edgeMat(:,2),:)); %new mid-edge points

                Vm=patchCentre(F,V);

                Vs = [V; Vn; Vm]; %Join point sets

                if qIter>1
                    CV=[CV; 1*ones(size(Vn,1),1); 2*ones(size(Vm,1),1);];
                else
                    CV=[0*ones(size(V,1),1); 1*ones(size(Vn,1),1); 2*ones(size(Vm,1),1);];
                end

            case 2
                edgeMat=[F(:,[1 2]); F(:,[3 4]);]; %Edges matrix
                Es=sort(edgeMat,2); %Sorted edges matrix
                [~,ind1,~]=unique(Es,'rows');
                edgeMat=edgeMat(ind1,:);

                numPoints = size(V,1);
                numEdges = size(edgeMat,1);

                % Get indices of the three edges associated with each face
                A = sparse(edgeMat(:,1),edgeMat(:,2),(1:numEdges)+numPoints,numPoints,numPoints,numEdges);
                A = max(A,A'); %Copy symmetric

                %Indices for A matrix
                indA_12=F(:,1)+(F(:,2)-1)*numPoints;
                indA_34=F(:,3)+(F(:,4)-1)*numPoints;

                %Get indices for vertex array
                indV_12=full(A(indA_12));
                indV_34=full(A(indA_34));

                %Create faces array
                Fs=[[F(:,1)  indV_12 indV_34 F(:,4)];...
                    [indV_12 F(:,2)  F(:,3)  indV_34]];

                %Create vertex array
                Vn=0.5*(V(edgeMat(:,1),:)+V(edgeMat(:,2),:)); %new mid-edge points

                Vs = [V; Vn; ]; %Join point sets

                if qIter>1
                    CV=[CV; 1*ones(size(Vn,1),1);];
                else
                    CV=[0*ones(size(V,1),1); 1*ones(size(Vn,1),1); ];
                end
            case 3
                edgeMat=[F(:,[2 3]); F(:,[4 1]);]; %Edges matrix
                Es=sort(edgeMat,2); %Sorted edges matrix
                [~,ind1,~]=unique(Es,'rows');
                edgeMat=edgeMat(ind1,:);

                numPoints = size(V,1);
                numEdges = size(edgeMat,1);

                % Get indices of the three edges associated with each face
                A = sparse(edgeMat(:,1),edgeMat(:,2),(1:numEdges)+numPoints,numPoints,numPoints,numEdges);
                A = max(A,A'); %Copy symmetric

                %Indices for A matrix
                indA_23=F(:,2)+(F(:,3)-1)*numPoints;
                indA_41=F(:,4)+(F(:,1)-1)*numPoints;

                %Get indices for vertex array
                indV_23=full(A(indA_23));
                indV_41=full(A(indA_41));

                %Create faces array
                Fs=[[F(:,1)  F(:,2) indV_23 indV_41];...
                    [indV_41 indV_23 F(:,3)  F(:,4)]];

                %Create vertex array
                Vn=0.5*(V(edgeMat(:,1),:)+V(edgeMat(:,2),:)); %new mid-edge points

                Vs = [V; Vn; ]; %Join point sets

                if qIter>1
                    CV=[CV; 1*ones(size(Vn,1),1);];
                else
                    CV=[0*ones(size(V,1),1); 1*ones(size(Vn,1),1); ];
                end

        end

        %Override input for looping
        C=repmat(C,[size(Fs,1)/size(F,1),1]);
        F=Fs;
        V=Vs;

    end
else
    CV=zeros(size(V,1),1);
end

varargout{1}=F;
varargout{2}=V;
varargout{3}=C;
varargout{4}=CV;

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function res = SymPoly2Sym(F,Dexp,MaxCtr)
%SYMPOLY2SYM   Similar to poly2sym but can handle symbolic coefficients and
%indeterminate variables specified as equations
%
%   SymPoly2Sym('a*z^2 + b*z + c','z') = [a b c]
%   SymPoly2Sym('a*sin(x)^2 + b*sin(x) + c','sin(x)') = [a b c]

%%This function iteratively takes the derivative then evaluates at 0 to
%%determine the coefficients

if nargin<2,
    Dexp = 't'; %Default dependent variable
end
F = sym(subs(F,Dexp,'XX')); %Substitute independent variable 'XX' for dependent expression Dexp

if nargin<3,
    MaxCtr = 20; %Default maximum order
end

res = sym([]); %Initialize res as empty sym
ctr = 0; %Initialize counter
while F~=0,
    res(end+1) = subs(F,'XX',0); %Evaluate at XX=0
    F = diff(F,'XX',1); %Take derivative

    %%%%Check ctr as a defense against recursive loops
    ctr = ctr+1;
    if ctr>=MaxCtr,
        disp(['Error: Exceeded maximum number of terms (' num2str(MaxCtr) ')'])
        disp(['Inputs may be in invalid form (not a polynomial)'])
        disp('To increase maximum number, use SymPoly2Sym(F,DV,Max#)')
        res = [];
        return
    end
    %%%%%%%%

end
Fact = [1 cumprod(1:length(res)-1)]; %Create factorial vector
res = res./Fact; %Divide coefficients by the amount contributed by derivative operations
res = res(end:-1:1); %Reverse order to match poly2sym convention

author: Sable

repo: mcbench-benchmarks


text:
function [yi, ypi] = trigint(x, y, xi, c)

% hRIGINT 1-D piecewise trigonometric interpolation
%    TRIGINT(X,Y,XI,C) interpolates to find YI, the values of the
%    underlying function Y at the points in the array XI, using
%    piecewise trigonometric interpolation.  X and Y must be vectors
%    of length N.
%
%    C specifies the number of data points to use in the
%    interpolation.  The default is to use all points.
%
%    [YI,YPI] = TRIGINT() also returns the interpolated derivative
%    of the underlying function Y at points XI.

% Joe Henning - Fall 2012

% On the Interpolation Trigonometric Polynomial with an Arbitrary Even Number of Nodes
% Ernest Scheiber
% 2011 13th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing
% 978-0-7695-4630-8/11 2011 IEEE
% DOI 10.1109/SYNASC.2011.13

if (nargin < 4)
   c = 0;
end

n = length(x);

if (n < c)
   fprintf('??? Bad c input to trigint ==> c <= length(x)\n');
   yi = [];
   ypi = [];
   return
end

for i = 1:length(xi)
   % Find the right place in the table by means of a bisection.
   klo = 1;
   khi = n;
   while (khi-klo > 1)
      k = fix((khi+klo)/2.0);
      if (x(k) > xi(i))
         khi = k;
      else
         klo = k;
      end
   end

   h = x(khi) - x(klo);
   if (h == 0.0)
      fprintf('??? Bad x input to trigint ==> x values must be distinct\n');
      yi(i) = NaN;
      ypi(i) = NaN;
      continue;
   end

   % Evaluate lagrange polynomial
   yi(i) = 0;
   ypi(i) = 0;
   if (c == 0)
      if (mod(n,2) == 0)   % even
         for k = 1:n
            sumx = 0;
            for m = 1:n
               sumx = sumx + x(m);
            end
            term = y(k)*(cos(0.5*(xi(i)-x(k))) + sin(0.5*(xi(i)-x(k)))*cot(0.5*sumx));
            termp = 0;
            termp2 = y(k)*0.5*(-sin(0.5*(xi(i)-x(k))) + cos(0.5*(xi(i)-x(k)))*cot(0.5*sumx));
            for m = 1:n
               if (k ~= m)
                  term = term*sin(0.5*(xi(i)-x(m)))/sin(0.5*(x(k)-x(m)));
                  prod = 0.5*cos(0.5*(xi(i)-x(m)));
                  for j = 1:n
                     if ((k ~= j) && (m ~= j))
                        prod = prod*sin(0.5*(xi(i)-x(j)))/sin(0.5*(x(k)-x(j)));
                     end
                  end
                  termp = termp + y(k)*(cos(0.5*(xi(i)-x(k))) + sin(0.5*(xi(i)-x(k)))*cot(0.5*sumx))*prod/sin(0.5*(x(k)-x(m)));
                  termp2 = termp2*sin(0.5*(xi(i)-x(m)))/sin(0.5*(x(k)-x(m)));
               end
            end
            yi(i) = yi(i) + term;
            ypi(i) = ypi(i) + termp + termp2;
         end
      else   % odd
         for k = 1:n
            term = y(k);
            termp = 0;
            for m = 1:n
               if (k ~= m)
                  term = term*sin(0.5*(xi(i)-x(m)))/sin(0.5*(x(k)-x(m)));
                  prod = 0.5*cos(0.5*(xi(i)-x(m)));
                  for j = 1:n
                     if ((k ~= j) && (m ~= j))
                        prod = prod*sin(0.5*(xi(i)-x(j)))/sin(0.5*(x(k)-x(j)));
                     end
                  end
                  termp = termp + y(k)*prod/sin(0.5*(x(k)-x(m)));
               end
            end
            yi(i) = yi(i) + term;
            ypi(i) = ypi(i) + termp;
         end
      end
   else
      if (mod(c,2) == 0)   % even
         c2 = c/2;
         if (klo < c2)
            klo = c2;
         end
         if (klo > n-c2)
            klo = n-c2;
         end
         khi = klo + 1;
         for k = klo-(c2-1):klo+c2
            sumx = 0;
            for m = klo-(c2-1):klo+c2
               sumx = sumx + x(m);
            end
            term = y(k)*(cos(0.5*(xi(i)-x(k))) + sin(0.5*(xi(i)-x(k)))*cot(0.5*sumx));
            termp = 0;
            termp2 = y(k)*0.5*(-sin(0.5*(xi(i)-x(k))) + cos(0.5*(xi(i)-x(k)))*cot(0.5*sumx));
            for m = klo-(c2-1):klo+c2
               if (k ~= m)
                  term = term*sin(0.5*(xi(i)-x(m)))/sin(0.5*(x(k)-x(m)));
                  prod = 0.5*cos(0.5*(xi(i)-x(m)));
                  for j = 1:n
                     if ((k ~= j) && (m ~= j))
                        prod = prod*sin(0.5*(xi(i)-x(j)))/sin(0.5*(x(k)-x(j)));
                     end
                  end
                  termp = termp + y(k)*(cos(0.5*(xi(i)-x(k))) + sin(0.5*(xi(i)-x(k)))*cot(0.5*sumx))*prod/sin(0.5*(x(k)-x(m)));
                  termp2 = termp2*sin(0.5*(xi(i)-x(m)))/sin(0.5*(x(k)-x(m)));
               end
            end
            yi(i) = yi(i) + term;
            ypi(i) = ypi(i) + termp + termp2;
         end
      else   % odd
         c2 = floor(c/2);
         if (klo < c2+1)
            klo = c2+1;
         end
         if (klo > n-c2)
            klo = n-c2;
         end
         khi = klo + 1;
         for k = klo-c2:klo+c2
            term = y(k);
            termp = 0;
            for m = klo-c2:klo+c2
               if (k ~= m)
                  term = term*sin(0.5*(xi(i)-x(m)))/sin(0.5*(x(k)-x(m)));
                  prod = 0.5*cos(0.5*(xi(i)-x(m)));
                  for j = 1:n
                     if ((k ~= j) && (m ~= j))
                        prod = prod*sin(0.5*(xi(i)-x(j)))/sin(0.5*(x(k)-x(j)));
                     end
                  end
                  termp = termp + y(k)*prod/sin(0.5*(x(k)-x(m)));
               end
            end
            yi(i) = yi(i) + term;
            ypi(i) = ypi(i) + termp;
         end
      end
   end
end


text:
function [maxVal, maxPos] = max(f)
%MAX   Global maximum of a DELTAFUN.
%   MAXVAL = MAX(F) returns the global maximum of the DELTAFUN F.
%
%   [MAXVAL, MAXPOS] = MAX(F) returns also a value such that MAXVAL = F(MAXPOS).
%
% See also MIN, MINANDMAX.

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

% To avoid code duplication, we simply call MINANDMAX():
[maxVal, maxPos] = minandmax(f);
maxVal = maxVal(2,:);
maxPos = maxPos(2,:);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function fx = p00_fun ( problem, option, nvar, x )

%*****************************************************************************80
%
%% P00_FUN evaluates the function for any problem.
%
%  Discussion:
%
%    These problems were collected by Professor Werner Rheinboldt, of the
%    University of Pittsburgh, and were used in the development of the
%    PITCON program.
%
%  Index:
%
%     1  The Freudenstein-Roth function
%     2  The Boggs function
%     3  The Powell function
%     4  The Broyden function
%     5  The Wacker function
%     6  The Aircraft stability function
%     7  The Cell kinetic function
%     8  The Riks mechanical problem
%     9  The Oden mechanical problem
%    10  Torsion of a square rod, finite difference solution
%    11  Torsion of a square rod, finite element solution
%    12  The materially nonlinear problem
%    13  Simpson's mildly nonlinear boundary value problem
%    14  Keller's boundary value problem
%    15  The Trigger Circuit
%    16  The Moore-Spence Chemical Reaction Integral Equation
%    17  The Bremermann Propane Combustion System
%    18  The semiconductor problem
%    19  The nitric acid absorption flash
%    20  The Buckling Spring
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Rami Melhem, Werner Rheinboldt,
%    A Comparison of Methods for Determining Turning Points of Nonlinear Equations,
%    Computing,
%    Volume 29, Number 3, September 1982, pages 201-226.
%
%    Werner Rheinboldt,
%    Numerical Analysis of Parameterized Nonlinear Equations,
%    Wiley, 1986,
%    ISBN: 0-471-88814-1,
%    LC: QA372.R54.
%
%    Werner Rheinboldt,
%    Sample Problems for Continuation Processes,
%    Technical Report ICMA-80-?,
%    Institute for Computational Mathematics and Applications,
%    Department of Mathematics,
%    University of Pittsburgh, November 1980.
%
%    Werner Rheinboldt, John Burkardt,
%    A Locally Parameterized Continuation Process,
%    ACM Transactions on Mathematical Software,
%    Volume 9, Number 2, June 1983, pages 215-235.
%
%    Werner Rheinboldt, John Burkardt,
%    Algorithm 596:
%    A Program for a Locally Parameterized
%    Continuation Process,
%    ACM Transactions on Mathematical Software,
%    Volume 9, Number 2, June 1983, pages 236-241.
%
%    Werner Rheinboldt,
%    Computation of Critical Boundaries on Equilibrium Manifolds,
%    SIAM Journal on Numerical Analysis,
%    Volume 19, Number 3, June 1982, pages 653-669.
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem index.
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, real X(NVAR), the argument of the function.
%
%    Output, real FX(NVAR-1), the value of the function at X.
%
  if ( problem == 1 )
    fx = p01_fun ( option, nvar, x );
  elseif ( problem == 2 )
    fx = p02_fun ( option, nvar, x );
  elseif ( problem == 3 )
    fx = p03_fun ( option, nvar, x );
  elseif ( problem == 4 )
    fx = p04_fun ( option, nvar, x );
  elseif ( problem == 5 )
    fx = p05_fun ( option, nvar, x );
  elseif ( problem == 6 )
    fx = p06_fun ( option, nvar, x );
  elseif ( problem == 7 )
    fx = p07_fun ( option, nvar, x );
  elseif ( problem == 8 )
    fx = p08_fun ( option, nvar, x );
  elseif ( problem == 9 )
    fx = p09_fun ( option, nvar, x );
  elseif ( problem == 10 )
    fx = p10_fun ( option, nvar, x );
  elseif ( problem == 11 )
    fx = p11_fun ( option, nvar, x );
  elseif ( problem == 12 )
    fx = p12_fun ( option, nvar, x );
  elseif ( problem == 13 )
    fx = p13_fun ( option, nvar, x );
  elseif ( problem == 14 )
    fx = p14_fun ( option, nvar, x );
  elseif ( problem == 15 )
    fx = p15_fun ( option, nvar, x );
  elseif ( problem == 16 )
    fx = p16_fun ( option, nvar, x );
  elseif ( problem == 17 )
    fx = p17_fun ( option, nvar, x );
  elseif ( problem == 18 )
    fx = p18_fun ( option, nvar, x );
  elseif ( problem == 19 )
    fx = p19_fun ( option, nvar, x );
  elseif ( problem == 20 )
    fx = p20_fun ( option, nvar, x );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_FUN - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized problem number = %d\n', problem );
    error ( 'P00_FUN - Fatal error!' );
  end
%
%  Ensure the result is a COLUMN vector.
%
  fx = fx(:);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = sinh(varargin)
%SINH   Hyperbolic sine of a DISKFUN.
%
% See also DISKFUN/SIN and DISKFUN/COSH.

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

[varargout{1:nargout}] = sinh@separableApprox(varargin{:});

end


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

%%
clear; close all; clc;

%% Syntax
% |[M,linearInd]=subImage(M,nSub);|

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


text:
function X = cholsolve(A,B)
% Solves (A'*A)*X=B with respect to X.
% A should be an upper triangular matrix (Cholesky decomposition).

opts.UT = true;
opts.TRANSA = true;
Z = linsolve(A, B, opts);
opts.TRANSA = false;
X = linsolve(A, Z, opts);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function hydro = combineBEM(hydro)
% Combines multiple BEM outputs into one hydrodynamic ‘system.’ This function
% requires that all BEM outputs have the same water depth, wave frequencies,
% and wave headings. This function would be implemented following multiple
% readWAMIT, readNEMOH, readCAPYTAINE, or readAQWA and before radiationIRF,
% radiationIRFSS, excitationIRF, writeBEMIOH5, or plotBEMIO function calls.
%
% See ``WEC-Sim\examples\BEMIO\NEMOH`` for examples of usage.
%
% Parameters
% ----------
%     hydro : [1 x n] struct
%         Structures of hydro data that will be combined into a single
%         structure
%
% Returns
% -------
%     hydro : [1 x 1] struct
%         Combined structure.

p = waitbar(0,'Combining multiple BEM results...');  % Progress bar

[m,n] = size(hydro);
tol = 1e-2;
if n>1
    for i = 2:n
        if max(abs(hydro(i).h-hydro(i-1).h)) > tol
            error('Error: Inconsistent water depth');
        elseif (size(hydro(i).theta) ~= size(hydro(i-1).theta)) | ...
                (max(abs(sort(hydro(i).theta)-sort(hydro(i-1).theta))) > tol)
            error('Error: Inconsistent wave headings');
        elseif (size(hydro(i).T) ~= size(hydro(i-1).T)) | ...
                (max(abs(sort(hydro(i).T)-sort(hydro(i-1).T))) > tol)
            error('Error: Inconsistent wave frequencies');
        end
    end

    n = n+1;
    hydro(n).theta = hydro(1).theta;
    hydro(n).body = [];
    hydro(n).Khs = [];
    hydro(n).cb = [];
    hydro(n).cg = [];
    hydro(n).code = [];
    hydro(n).dof = [];
    hydro(n).ex_im = [];
    hydro(n).ex_ma = [];
    hydro(n).ex_ph = [];
    hydro(n).ex_re = [];
    hydro(n).fk_im = [];
    hydro(n).fk_ma = [];
    hydro(n).fk_ph = [];
    hydro(n).fk_re = [];
    hydro(n).sc_im = [];
    hydro(n).sc_ma = [];
    hydro(n).sc_ph = [];
    hydro(n).sc_re = [];
    hydro(n).file = [];
    hydro(n).g = 9.807;
    hydro(n).h = hydro(1).h;
    hydro(n).Nb = 0;
    hydro(n).Nf = hydro(1).Nf;
    hydro(n).Nh = hydro(1).Nh;
    hydro(n).rho = 1000;
    hydro(n).T = hydro(1).T;
    hydro(n).Vo = [];
    hydro(n).w = hydro(1).w;
    for i = 1:(n-1)
        hydro(n).body = [hydro(n).body, hydro(i).body];
        hydro(n).Khs = cat(3, hydro(n).Khs, hydro(i).Khs);
        hydro(n).cb = [hydro(n).cb, hydro(i).cb];
        hydro(n).cg = [hydro(n).cg, hydro(i).cg];
        hydro(n).code = [hydro(n).code, hydro(i).code];
        hydro(n).dof = [hydro(n).dof, hydro(i).dof];
        hydro(n).ex_im = cat(1, hydro(n).ex_im, hydro(i).ex_im);
        hydro(n).ex_ma = cat(1, hydro(n).ex_ma, hydro(i).ex_ma);
        hydro(n).ex_ph = cat(1, hydro(n).ex_ph, hydro(i).ex_ph);
        hydro(n).ex_re = cat(1, hydro(n).ex_re, hydro(i).ex_re);
        hydro(n).fk_im = cat(1, hydro(n).fk_im, hydro(i).fk_im);
        hydro(n).fk_ma = cat(1, hydro(n).fk_ma, hydro(i).fk_ma);
        hydro(n).fk_ph = cat(1, hydro(n).fk_ph, hydro(i).fk_ph);
        hydro(n).fk_re = cat(1, hydro(n).fk_re, hydro(i).fk_re);
        hydro(n).sc_im = cat(1, hydro(n).sc_im, hydro(i).sc_im);
        hydro(n).sc_ma = cat(1, hydro(n).sc_ma, hydro(i).sc_ma);
        hydro(n).sc_ph = cat(1, hydro(n).sc_ph, hydro(i).sc_ph);
        hydro(n).sc_re = cat(1, hydro(n).sc_re, hydro(i).sc_re);
        hydro(n).file = [hydro(n).file, hydro(i).file];
        hydro(n).Nb = hydro(n).Nb + hydro(i).Nb;
        hydro(n).Vo = [hydro(n).Vo, hydro(i).Vo];
    end
    hydro(n).Ainf = zeros(sum(hydro(n).dof),sum(hydro(n).dof));
    hydro(n).A = zeros(sum(hydro(n).dof),sum(hydro(n).dof),hydro(n).Nf);
    hydro(n).B = zeros(sum(hydro(n).dof),sum(hydro(n).dof),hydro(n).Nf);
    if isfield(hydro,'gbm')==1
        hydro(n).gbm = zeros(sum(hydro(n).dof),sum(hydro(n).dof),4);
    end
    k = 0;
    for i = 1:(n-1)
        j = sum(hydro(i).dof);
        hydro(n).Ainf((k+1):(k+j),(k+1):(k+j)) = hydro(i).Ainf;
        hydro(n).A((k+1):(k+j),(k+1):(k+j),:) = hydro(i).A;
        hydro(n).B((k+1):(k+j),(k+1):(k+j),:) = hydro(i).B;
        if isfield(hydro,'gbm')==1
            hydro(n).gbm((k+1):(k+j),(k+1):(k+j),:) = hydro(i).gbm;
        end
        k = k + j;
    end

    hydro(1) = hydro(n);
    hydro(2:n) = [];

end
close(p);
end

author: WEC-Sim

repo: WEC-Sim

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


text:
function rs = medianfilt(s, len)

%tstoolbox/@signal/medianfilt
%   Syntax:
%     * rs = medianfilt(s, len)
%
%   Moving median filter of width len samples for a scalar time series
%   (len should be odd).
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(2,2);

c = medianfilt(s.core, len);
rs = signal(c, s);				% special constructor calling syntax for working routines

a = getaxis(s,1);
a = setfirst(a, ((len+1)/2) * delta(a));

rs = setaxis(rs, 1, a);
rs = addhistory(rs,  ['Moving average with window of ' num2str(len) ' samples'] );
rs = addcommandlines(rs, 's = medianfilt(s', len);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%% Example Scene
%
% This script demonstrates using a combination of functions to create
% terrain, add a sky, set some lighting options, and provide the ability to
% rotate around the figure. It uses several of MATLAB's 3D graphics
% commands, including patches, lighting, camera positioning, and
% projection. If any of these are unfamiliar, the documentation covers them
% well. MATLAB can produce much more than simple plots.
%
% This script is likely to take between 3 and 30 seconds to complete on
% most computers.
%
% Tucker McClure
% Copyright 2012, The MathWorks, Inc.

%%
% Create or clear a figure for the rendering.
figure(1);
clf();

%% Land and Sea

% Create the terrain height map.
[~, ~, ~, hm, xm, ym] = generate_terrain(7, 513, 0, 0.1);

% Generate appropriate colors.
cm = generate_terrain_colors(hm);

% Flatten the oceans.
hmp = max(hm, 0);

% Draw the land and calculate its normal vectors.
h_land = patch(surf2patch(xm, ym, hmp, cm));
[nx, ny, nz] = surfnorm(hmp);
land_normals = [nx(:), ny(:), nz(:)];

% Set it's material properties for interpolated colors and appropriate
% lighting.
set(h_land, 'VertexNormals',    land_normals, ...
            'DiffuseStrength',  0.8, ...      % Reacts to light direction
            'SpecularStrength', 0, ...        % Not shiny
            'AmbientStrength',  0.3, ...      % Reacts to ambient light
            'BackFaceLighting', 'unlit');     % Don't illuminate reverse

% Add a black backdrop beneath the land. Sometimes a vertex seems to miss
% its position slightly, and sky colors creep through, and that's weird.
patch('Faces',           [1 2 3 4], ...
      'Vertices',        [-1 -1 -0.01; ...
                           1 -1 -0.01; ...
                           1  1 -0.01; ...
                          -1  1 -0.01], ...
      'FaceVertexCData', 0.25*ones(4, 3));

%% Sky

% Draw a random time of day between sunrise and sunset. Get the
% corresponding light color. The sky is always the color of the sun at noon
% (and is then affected by the sun color at the current time).
current_time = 0.24 + 0.52 * rand();
sun_color    = sun_tones(current_time);
sky_color    = sun_tones(0.5);

% Create the sky patch by scaling a sphere. Note that the nothing in the
% scene can be rendered outside of [-100, 100] on any axis due to a
% documented rendering bug. Therefore, the sphere will be scaleld up to
% 99, which is plenty far away since the terrain is limited to [-1, 1].
[xs, ys, zs] = sphere(16);
sky_scale = 99;
sky_patch = surf2patch(sky_scale*xs, sky_scale*ys, sky_scale*zs, ...
                       repmat(reshape(sky_color, [1 1 3]), [17 17 1]));
h_sky = patch(sky_patch);

% Set appropriate lighting options for the sky.
set(h_sky, 'DiffuseStrength',  0.3, ...
           'SpecularStrength', 0, ...
           'AmbientStrength',  1, ...
           'BackFaceLighting', 'unlit');

% Vertex normals are the opposite of what I expect. Reverse them.
set(h_sky,  'VertexNormals', -get(h_sky, 'VertexNormals'));

%% Sun

% Create a light for the sun.
h_light = lightangle(90, 360*(current_time - 0.25));

% Set the light's color.
set(h_light, 'Color', sun_color);

% Set the ambient color used in the scene.
set(gca(), 'AmbientLightColor', sun_color);

% Use decent lighting.
lighting gouraud;

%% Axes and Global Settings

% Set axes options.
camera_target   = [0 0 mean(hmp(:)) + 0.5*std(hmp(:))];
camera_position = camera_target + [1.15 0 0.5];
set(gca, 'DataAspectRatio', [1 1 1], ...
         'Visible',         'off', ...
         'Projection',      'Perspective', ...
         'Position',        [0 0 1 1], ...
         'CameraTarget',    camera_target, ...
         'CameraViewAngle', 45, ...
         'CameraUpVector',  [0 0 1], ...
         'CameraPosition',  camera_position, ...
         'XLim',            [-100 100], ...
         'YLim',            [-100 100], ...
         'ZLim',            [-100 100]);
camorbit(360*rand(), 0);

% Patches should have no edges and interpolated face colors.
shading interp;

% Force the drawing.
drawnow();

% Add FigureRotator to allow the user to rotate around the middle of the
% scene and zoom in and out using a mouse. This is available in the extras
% zip file or from File Exchange. Use |unzip terrain_generation_extras.zip|
% in this directory to extract FigureRotator along with an example file.
if exist('FigureRotator', 'class');
    f = FigureRotator(gca());
end


text:
function H = comp_nyquistfilt(wintype,fs,chan_max,freqtoscale,scaletofreq,bwmul,bins,Ls)
%COMP_NYQUISTFILT high-pass filter for warped filter banks

    kk = chan_max;
    while scaletofreq(kk-bwmul) < fs/2;
      kk = kk+1/bins;
    end
    Maxfilt = kk;

    Minpos = ceil(Ls/fs*scaletofreq(chan_max+1/bins-bwmul));
    samples = freqtoscale((Minpos-1:floor(Ls/2))*fs/Ls);

    FILTS = zeros(round(bins*(Maxfilt-chan_max)),numel(samples));
    for kk = 1:size(FILTS,1)
       FILTS(kk,:) = firwin(wintype,(samples-(chan_max+kk/bins))/(2*bwmul));
    end
    H = zeros(2*numel(samples)-1,1);
    H(1:numel(samples)) = sqrt(sum(abs(FILTS.^2),1));
    H(numel(samples)+1:end) = H(numel(samples)-1:-1:1);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%--- help for dsge/simulate_nonlinear ---
%
%  INTERNAL FUNCTION
%
%

author: jmaih


text:
function rhs = globRHSIFE3DFace(fun, funB, fem, femIF, d1)

%% USAGE: generate global matrix on a 3D mesh Face
%
% INPUTS:
% fun --- coefficient function
% mesh --- a struct data contains very rich mesh information.
% fem1 --- global DoF for test function space
% fem2 --- global DoF for trial function space
% d1 --- derivative info for test function
% d2 --- derivative info for trial function
%            d = [0,0,0]: function value
%            d = [1,0,0]: Dx value
%            d = [0,1,0]: Dy value
%            d = [0,0,1]: Dz value
%
% OUTPUTS:
% matrix --- global (mass, stiffness, ...) matrix.

% Last Modified: 08/02/2020 by Xu Zhang

%% 1. RHS on Boundary Interface Faces
dof = fem.ldof;
if strcmp(fem.type,'P1')||strcmp(fem.type,'DGP1')||strcmp(fem.type,'CR')
    feEvalBas = @evalP1Bas3D;
elseif strcmp(fem.type,'P2')||strcmp(fem.type,'DGP2')
    feEvalBas = @evalP2Bas3D;
end
rhs = zeros(length(fem.p),1);

if size(femIF.tB,1) ~= 0

    nmB = femIF.normalB;
    AB = femIF.areaB; gxB = femIF.gxB; gyB = femIF.gyB; gzB = femIF.gzB; gw = femIF.gw;
    nt = size(femIF.tB,1); % not number of interface element, but quadrature element
    ntot = nt*dof;
    IB = zeros(ntot,1); XB = zeros(ntot,1);

    coef = feval(fun,gxB,gyB,gzB);
    bval = feval(funB,gxB,gyB,gzB);
    IbasB = cell(dof,1);
    if d1 == 0
        for i = 1:dof
            IbasB{i} = feEvalBas(femIF.basB(:,:,i), gxB, gyB, gzB, [0,0,0]);
        end
    elseif d1 == 1
        for i = 1:dof
            IbasBx = feEvalBas(femIF.basB(:,:,i), gxB, gyB, gzB, [1,0,0]);
            IbasBy = feEvalBas(femIF.basB(:,:,i), gxB, gyB, gzB, [0,1,0]);
            IbasBz = feEvalBas(femIF.basB(:,:,i), gxB, gyB, gzB, [0,0,1]);
            IbasB{i} = IbasBx.*nmB(:,1) + IbasBy.*nmB(:,2) + IbasBz.*nmB(:,3);
        end
    end

    ind = 0;
    for i = 1:dof
        IB(ind+1:ind+nt) = femIF.tB(:,i);
        XB(ind+1:ind+nt) = AB.*sum((((IbasB{i}.*coef).*bval).*gw'),2);
        ind = ind + nt;
    end
    ID = find(XB~=0);
    rhs = sparse(IB(ID),1,XB(ID),length(fem.p),1);
end

repo: ifem


text:
function stats = applyStatsfun(problem, x, storedb, key, options, stats)
% Apply the statsfun function to a stats structure (for solvers).
%
% function stats = applyStatsfun(problem, x, storedb, key, options, stats)
%
% Applies the options.statsfun user supplied function (if it was provided)
% to the stats structure, and returns the (possibly) modified stats
% structure.
%
% storedb is a StoreDB object, key is the StoreDB key to point x.
%
% Note: if statsfun accepts a store structure as input, this structure can
% be read but not modified (modifications will be lost) ; the store
% structure will contain the store.shared field.
%
% See also: applyHook

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, April 3, 2013.
% Contributors:
% Change log:
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.
%
%   July 19, 2020 (NB):
%       Creates a field called 'hooked' in stats structures, for applyHook.

    if isfield(options, 'statsfun')

        switch nargin(options.statsfun)
            case 3
                stats = options.statsfun(problem, x, stats);
            case 4
                % Obtain, pass along, and save the store for x.
                % get/setWithShared must come in pairs.
                store = storedb.getWithShared(key);
                stats = options.statsfun(problem, x, stats, store);
                storedb.setWithShared(store, key);
            otherwise
                warning('manopt:statsfun', ...
                        'statsfun unused: wrong number of inputs');
        end
    end

    % Always create a field called 'hooked' for the benefit of applyHook.
    % This is to be set to either true or false by applyHook.
    stats.hooked = NaN;

end

author: NicolasBoumal

repo: manopt


text:
function r8utp_print_some ( n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R8UTP_PRINT_SOME prints some of an R8UTP matrix.
%
%  Discussion:
%
%    The R8UTP storage format is appropriate for an upper triangular
%    matrix.  Only the upper triangle of the matrix is stored,
%    by successive partial columns, in an array of length (N*(N+1))/2,
%    which contains (A11,A12,A22,A13,A23,A33,A14,...,ANN)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, real A((N*(N+1))/2), the matrix.
%
%    Input, integer ILO, JLO, IHI, JHI, the first row and
%    column, and the last row and column to be printed.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  incx = 5;
%
%  Print the columns of the matrix, in strips of 5.
%
  for j2lo = jlo : incx : jhi

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    inc = j2hi + 1 - j2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );

    for j = j2lo : j2hi
      j2 = j + 1 - j2lo;
      fprintf ( 1, '%7d       ', j );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );
    fprintf ( 1, '  ---\n' );
%
%  Determine the range of the rows in this strip.
%
    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, n );

    for i = i2lo : i2hi

      fprintf ( 1, '%4d  ', i );
%
%  Print out (up to) 5 entries in row I, that lie in the current strip.
%
      for j2 = 1 : inc

        j = j2lo - 1 + j2;

        if ( i <= j )
          aij = a(i+(j*(j-1))/2);
        else
          aij = 0.0;
        end

        fprintf ( 1, '%12g  ', aij );

      end

      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
function varargout = gallery(name)
%CHEB.GALLERY   Chebfun example functions.
%   F = CHEB.GALLERY(NAME) returns a chebfun or a quasimatrix corresponding to
%   NAME.  See the listing below for available names.
%
%   For example,  plot(cheb.gallery('zigzag'))  plots a degree 10000 polynomial
%   that doesn't look like a polynomial, and  plot(cheb.gallery('gamma'))  shows
%   a chebfun with poles. For details of how each function is constructed, try
%   type +cheb/gallery  or  edit cheb.gallery.
%
%   [F,FA] = CHEB.GALLERY(NAME) also returns the anonymous function FA used to
%   define the function. Some gallery functions are generated by operations
%   beyond the usual Chebfun constructor (e.g. by solving ODEs), so FA in those
%   cases simply evaluates the chebfun.
%
%   CHEB.GALLERY with no input argument returns a random function from the
%   gallery.
%
%   CHEB.GALLERY with no output argument creates a plot of the selected
%   function.
%
%   airy         Airy Ai function on [-40,40]
%   bessel       Bessel function J_0 on [-100,100]
%   bump         C-infinity function with compact support
%   blasius      Blasius function on [0,10]
%   chirp        Sine with exponentially increasing frequency
%   daubechies   Approximation to Daubechies phi_2 wavelet scaling function
%   erf          Error function on [-10,10]
%   fishfillet   Wild oscillations from Extreme Extrema example
%   gamma        Gamma function on [-4,4]
%   gaussian     Gaussian function on [-Inf,Inf]
%   jitter       A piecewise constant function generated by ROUND
%   kahaner      Challenging integrand with four spikes
%   motto        Chebfun motto (Gilbert Strang)
%   random       Polynomial interpolant through random data in Chebyshev points
%   rose         A complex-valued sinusoid
%   runge        Runge function
%   seismograph  Tanh plus growing oscillation
%   Si           Sine integral on [-50,50]
%   sinefun1     As smooth as it looks
%   sinefun2     Not as smooth as it looks
%   spikycomb    25 peaks, each sharper than the last
%   stegosaurus  max(wiggly, x/10)
%   vandercheb   Chebyshev-Vandermonde quasimatrix
%   vandermonde  Vandermonde quasimatrix
%   wiggly       One of the Chebfun team's favorites
%   wild         An iteratively defined function on [-1 1]
%   zigzag       Degree 10000 polynomial that looks piecewise linear
%
%   Gallery functions are subject to change in future releases of Chebfun.
%
% See also CHEB.GALLERYTRIG, CHEB.GALLERY2, CHEB.GALLERY3, CHEB.GALLERYDISK, CHEB.GALLERYSPHERE.

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

% If the user did not supply an input, return a random function from the
% gallery.
if ( nargin == 0 )
    % NOTE [AB, 2014/12/01]: This is not a particularly scalable way of
    % implementing cheb.gallery. A better approach would be that each gallery
    % function lived in a separate m-file, and could return it's name,
    % description, and the required outputs, F and FA. That way, we would get
    % rid of the switch statement below. Furthermore, this would make the
    % testing of gallery easier, as we wouldn't have to update the list of
    % input options there.
    names = {'airy', 'bessel', 'blasius', 'bump', ...
        'chirp', 'daubechies', 'erf', 'fishfillet', 'gamma', 'gaussian', ...
        'jitter', 'kahaner', 'motto', 'random', 'rose', 'runge', ...
        'seismograph', 'Si', 'sinefun1', 'sinefun2', 'spikycomb', ...
        'stegosaurus', 'vandercheb', 'vandermonde', 'wiggly', 'wild', 'zigzag'};
    name = names{randi(length(names))};
end

% If nargout == 0, then the function is plotted. Each function in the gallery
% has its own plotting preferences.
ylims = [];
axispref = {};

% The main switch statement.
switch lower(name)

    % Airy Ai function on [-40,40]:
    case 'airy'
        fa = @airy;
        f = chebfun(fa, [-40 40]);
        ylims = [-.75 .75];

    % Bessel function with parameter 0 on [-100,100]:
    case 'bessel'
        fa = @(x) besselj(0, x);
        f = chebfun(fa, [-100 100]);
        ylims = [-.5 1.1];

    % Blasius function on [0,10]:
    case 'blasius'
        op = @(u) 2*diff(u,3) + u.*diff(u,2);
        bc  = @(x,u) [u(0); feval(diff(u),0); feval(diff(u),10)-1];
        N = chebop(op, [0 10], bc);
        N.init = chebfun([3.55542; 4.25907; 0.43669; -0.21367;
                0.06382; -0.00118; -0.00865; 0.00306], [0 10], 'coeffs');
        f = N\0;
        fa = @(x) f(x);
        ylims = [-.5 9];

    % Bump function:
    case 'bump'
        fa = @(x) (abs(x) < 1).*exp(-1./(1-x.^2));
        f = chebfun(fa, [-2 2]);

    % Sine with exponentially increasing frequency:
    case 'chirp'
        fa = @(x) sin(x.*exp(x));
        f = chebfun(fa, [0 5]);
        ylims = [-2 2];

    % Approx to Daubechies phi_2 wavelet scaling function:
    case 'daubechies'
        f = daubechies(10);
        fa = @(x) f(x);
        ylims = [-.5 1.5];

    % Error function on [-10,10]:
    case 'erf'
        fa = @erf;
        f = chebfun(fa, [-10 10]);

    % Wild oscillations from Extreme Extrema example:
    case 'fishfillet'
        fa = @(x) cos(x).*sin(exp(x));
        f = chebfun(fa, [0 6]);
        ylims = [-1.2 1.2];

    % Gamma function on [-4,4]:
    case 'gamma'
        fa = @gamma;
        f = chebfun(fa, [-4 4], 'blowup', 'on', 'splitting', 'on');

    % Gaussian function on [-Inf,Inf]:
    case 'gaussian'
        fa = @(x) exp(-x.^2/2)/sqrt(2*pi);
        f = chebfun(fa, [-Inf Inf]);
        ylims = [-.05 .45];

    % A piecewise constant function generated by ROUND:
    case 'jitter'
        fa = @(x) round(exp(x)*2.*sin(8*x));
        f = chebfun(fa, 'splitting', 'on');
        ylims = [-5 6];

    % Challenging integrand with four spikes:
    case 'kahaner'
        fa = @(x) sech(10*(x-0.2)).^2 + sech(100*(x-0.4)).^4 + ...
         sech(1000*(x-0.6)).^6 + sech(1000*(x-0.8)).^8;
        f = chebfun(fa, [0 1]);
        ylims = [-.1 1.2];

    % (Scribbled) Chebfun motto by Gilbert Strang:
    case 'motto'
        f = exp(3i*scribble('there is no fun like chebfun'));
        fa = @(x) f(x);
        axispref = {'equal', [-1 1 -1 1]*1.2, 'off'};

    % Polynomial interpolant through random data in Chebyshev points:
    case 'random'
        f = chebfun(rand(100,1));
        fa = @(x) f(x);

    % Rose curve:
    case 'rose'
        m = 5;
        n = 4;
        fa = @(t) cos(m/n*t).*cos(t) + 1i*cos(m/n*t).*sin(t);
        f = chebfun(fa, [0, 8*pi], 'trig');
        axispref = {'equal', [-1 1 -1 1]*1.2};

    % Runge function:
    case 'runge'
        fa = @(x) 1./(1 + x.^2);
        f = chebfun(fa, [-5 5]);
        ylims = [0 1.2];

    % Tanh plus growing oscillation from ATAP, Chapter 5:
    case 'seismograph'
        fa = @(x) tanh(20*sin(12*x)) + .02*exp(3*x).*sin(300*x);
        f = chebfun(fa);

    % Sine integral:
    case 'si'
        f = cumsum(chebfun(@(x) sin(x)./(x), [-50, 50]));
        fa = @(x) f(x);

    % As smooth as it looks:
    case 'sinefun1'
        fa = @(x) (1.75 + sin(50*x));
        f = chebfun(fa);

    % Not as smooth as it looks:
    case 'sinefun2'
        fa = @(x) (1.75 + sin(50*x)).^1.0001;
        f = chebfun(fa);

    % 25 peaks, each sharper than the last from ATAP, Chapter 18:
    case 'spikycomb'
        fa = @(x) exp(x).*sech(4*sin(40*x)).^exp(x);
        f = chebfun(fa);

    % max(wiggly, x/10):
    case 'stegosaurus'
        fa = @(x) max(sin(x)+sin(x.^2), x/10);
        f = chebfun(fa, [0 10], 'splitting', 'on');
        ylims = [-.2 2.2];

    % Chebyshev-Vandermonde quasimatrix:
    case 'vandercheb'
        f = chebpoly(0:5);
        f = simplify(cheb2quasi(f));
        fa = @(x) f(x);
        ylims = [-1.2 1.2];

    % Vandermonde quasimatrix:
    case 'vandermonde'
        fa = @(x) x.^(0:5);
        f = chebfun(fa, 'vectorize');
        f = simplify(cheb2quasi(f));
        ylims = [-1.2 1.2];

    % One of the Chebfun team's favorites:
    case 'wiggly'
        fa = @(x) sin(x) + sin(x.^2);
        f = chebfun(fa, [0 10]);
        ylims = [-2.4 2.4];

    % An example from one of the first Chebfun papers (Trefethen, 2007):
    case 'wild'
        fa = @(x) wild(x); % Defined below.
        f = chebfun(fa, [-1 1]);
        ylims = [5.5 9.5];

    % Degree 10000 polynomial that looks piecewise linear from ATAP appendix:
    case 'zigzag'
        f = cumsum(chebfun(@(t) sign(sin(100*t./(2-t))), 10000));
        fa = @(x) f(x);

    % Raise an error if the input is unknown.
    otherwise
        error('CHEB:GALLERY:unknown:unknownFunction', ...
            'Unknown function.')
end

% Only return something if there is an output argument.
if ( nargout > 0 )
    varargout = {f, fa};
else
    % Otherwise, plot the function.
    plot(f)
    title([name ', length = ' num2str(length(f))])
    if ( ~isempty(ylims) )
        ylim(ylims)
    end
    if ( ~isempty(axispref) )
        axis(axispref{:})
    end
    shg
end

end

function f = daubechies(n)
% nth order polynomial approximation to Daubechies scaling function phi_2
[x, phi] = daub(n);
f = chebfun([phi'; 0], [0 3], 'equi');

    function [x, phi] = daub(n)
        s = sqrt(3);
        if ( n == 0 )
            % Base case of recursion.
            x = 0:2;
            phi = [0, (1+s)/2, (1-s)/2];
        else
            c = [1+s, 3+s, 3-s, 1-s]/4;
            % Recursively call the daub() method.
            [x2, phi2] = daub(n-1);
            pp = [phi2, 0*phi2];
            N = 2^(n-1);
            ii = 1:6*N;
            phi = c(1)*pp(ii);
            for k = 2:4
                ii = ii([(5*N + 1):6*N, 1:5*N]);
                phi = phi + c(k)*pp(ii);
            end
            x = [x2, x2+3]/2;
        end
    end

end

function s = wild(x)
% The 'wild' function from Computing Numerically with Functions, Trefethen 2007.
f = sin(pi*x);
s = f;
for j = 1:15
    f = (3/4)*(1 - 2*f.^4);
    s = s + f;
end
end


text:
function [am,em]=lpcrr2am(rr);
%LPCRR2AM Convert autocorrelation coefs to ar coef matrix [AM,EM]=(RR)
%AM is a 3-dimensional matrix of size (p+1,p+1,nf) where p is the lpc order
%and nf the number of frames.
%The matrix AM(:,:,*) is upper triangular with 1's on the main diagonal
%and contains the lpc coefficients for all orders from p down to 0.
%
%For lpc order p+1-r, AM(r,r:p+1,*), AM(p+1:-1:r,p+1,*) and EM(*,r) contain
%the lpc coefficients, reflection coefficients and the residual energy respectively.
%
%If A=am(:,:,*), R=toeplitz(rr(*,:)) and E=diag(em(*,:)), then
% A*R*A'=E; inv(R)=A'*(1/E)*A; A*R is lower triangular with the same diagonal as E
%
% This routine is equivalent to: c=chol(inv(toeplitz(rr))); d=diag(c).^-1; em=d.^2; am=diag(d)*c

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: lpcrr2am.m,v 1.4 2007/05/04 07:01:39 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p1]=size(rr);
p=p1-1;
p2=p1+1;
am=zeros(nf,p1,p1);
em=zeros(nf,p1);
am(:,p1,p1)=1;
em(:,p1)=rr(:,1);
ar=ones(nf,p1);
ar(:,2) = -rr(:,2)./rr(:,1);
e = rr(:,1).*(ar(:,2).^2-1);
for n = 2:p
   q=p2-n;
   em(:,q)=-e;
   am(:,q:p1,q)=ar(:,1:n);
   k = (rr(:,n+1)+sum(rr(:,n:-1:2).*ar(:,2:n),2)) ./ e;
   ar(:,2:n) = ar(:,2:n)+k(:,ones(1,n-1)).*ar(:,n:-1:2);
   ar(:,n+1) = k;
   e = e.*(1-k.^2);
end
em(:,1)=-e;
am(:,:,1)=ar;
am=permute(am,[3 2 1]);

author: decouples


text:
function linpack_c_test04 ( )

%*****************************************************************************80
%
%% TEST04 tests CGBCO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  ml = 1;
  mu = 1;
  n = 3;
  lda = 2*ml+mu+1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04\n' );
  fprintf ( 1, '  For a complex general band storage matrix:\n' );
  fprintf ( 1, '  CGBCO factors the matrix and estimates the\n' );
  fprintf ( 1, '  reciprocal condition number.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix order is N = %d\n', n );
  fprintf ( 1, '  The lower band is ML =  %d\n', ml );
  fprintf ( 1, '  The upper band is MU =  %d\n', mu );
%
%  Set the values of the matrix A.
%
  a_save(1:n,1:n) = 0.0;

  m = ml + mu + 1;

  seed = 123456789;

  for j = 1 : n
    i1 = max ( 1, j - mu );
    i2 = min ( n, j + ml );
    for i = i1 : i2
      k = i - j + m;
      [ a(k,j), seed ] = c4_uniform_01 ( seed );
      a_save(i,j) = a(k,j);
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix A:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  (%8f  %8f)', real ( a_save(i,j) ), imag ( a_save(i,j) ) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Factor the matrix A.
%
  [ a, ipvt, rcond ] = cgbco ( a, lda, n, ml, mu );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimated reciprocal condition RCOND = %f\n', rcond );

  return
end


text:
function dtiResetMouseImage(handles)
% Reset callbacks on the image planes in the main window
%
%    dtiResetMouseImage(handles)
%
% When the image is redrawn, matlab resets the callbacks.  So, we reset
% them again
%
% See also:  dtiRefreshFigure
%
% Brian (c) Stanford VISTASOFT Team, 2012

set(handles.x_cut_img,'HitTest','off');
set(handles.x_cut,'ButtonDownFcn','dtiFiberUI(''x_cut_click_Callback'',gcbo,[],guidata(gcbo))');

set(handles.y_cut_img,'HitTest','off');
set(handles.y_cut,'ButtonDownFcn','dtiFiberUI(''y_cut_click_Callback'',gcbo,[],guidata(gcbo))');

set(handles.z_cut_img,'HitTest','off');
set(handles.z_cut,'ButtonDownFcn','dtiFiberUI(''z_cut_click_Callback'',gcbo,[],guidata(gcbo))');

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ r, pc ] = tetrahedron_insphere_3d ( tetra )

%*****************************************************************************80
%
%% TETRAHEDRON_INSPHERE_3D finds the insphere of a tetrahedron in 3D.
%
%  Discussion:
%
%    The insphere of a tetrahedron is the inscribed sphere, which touches
%    each face of the tetrahedron at a single point.
%
%    The points of contact are the centroids of the triangular faces
%    of the tetrahedron.  Therefore, the point of contact for a face
%    can be computed as the average of the vertices of that face.
%
%    The sphere can then be determined as the unique sphere through
%    the four given centroids.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Schneider, David Eberly,
%    Geometric Tools for Computer Graphics,
%    Elsevier, 2002,
%    ISBN: 1558605940,
%    LC: T385.G6974.
%
%  Parameters:
%
%    Input, real TETRA(3,4), the vertices of the tetrahedron.
%
%    Output, real R, PC(3,1), the radius and the center
%    of the sphere.
%
  dim_num = 3;

  v21(1:dim_num,1) = tetra(1:dim_num,2) - tetra(1:dim_num,1);
  v31(1:dim_num,1) = tetra(1:dim_num,3) - tetra(1:dim_num,1);
  v41(1:dim_num,1) = tetra(1:dim_num,4) - tetra(1:dim_num,1);
  v32(1:dim_num,1) = tetra(1:dim_num,3) - tetra(1:dim_num,2);
  v42(1:dim_num,1) = tetra(1:dim_num,4) - tetra(1:dim_num,2);
  v43(1:dim_num,1) = tetra(1:dim_num,4) - tetra(1:dim_num,3);

  n123 = r8vec_cross_product_3d ( v21, v31 );
  n124 = r8vec_cross_product_3d ( v41, v21 );
  n134 = r8vec_cross_product_3d ( v31, v41 );
  n234 = r8vec_cross_product_3d ( v42, v32 );

  l123 = norm ( n123 );
  l124 = norm ( n124 );
  l134 = norm ( n134 );
  l234 = norm ( n234 );

  pc(1:dim_num,1) = ( l234 * tetra(1:dim_num,1)   ...
                    + l134 * tetra(1:dim_num,2)   ...
                    + l124 * tetra(1:dim_num,3)   ...
                    + l123 * tetra(1:dim_num,4) ) ...
                  / ( l234 + l134 + l124 + l123 );

  b(1:dim_num,1:4) = tetra(1:dim_num,1:4);
  b(4,1:4) = 1.0;

  gamma = abs ( r8mat_det_4d ( b ) );

  r = gamma / ( l234 + l134 + l124 + l123 );

  return
end


text:
function [V,E,P,BE,CE,PE] = readTGF(filename)
  % READTGF
  %
  % [V,E,P,BE,CE,PE] = readTGF(filename)
  %
  % Read a graph from a .tgf file
  %
  % Input:
  %  filename  .tgf file name
  % Ouput:
  %  V  # vertices by 3 list of vertex positions
  %  E  # edges by 2 list of edge indices
  %  P  # point-handles list of point handle indices
  %  BE # bone-edges by 2 list of bone-edge indices
  %  CE # cage-edges by 2 list of cage-edge indices
  %  PE # pseudo-edges by 2 list of pseudo-edge indices
  %
  % Assumes that graph vertices are 3 dimensional

  V = [];
  E = [];
  BE = [];
  CE = [];
  PE = [];
  fp = fopen(filename,'r');
  % read next whole line
  line = fscanf(fp,' %[^\n]s');
  % read vertices until seeing separator
  while(line(1,1) ~= '#')
    [v,count] = sscanf(line,'%d %g %g %g',4);
    % first number is always index
    v = v(2:min(count,4));
    V = [V;v'];
    line = fscanf(fp,' %[^\n]s');
  end
  % read edges until file is through
  line = fscanf(fp,' %[^\n]s');
  while(sum(size(line)) > 0 && line(1,1) ~= '#')
    [e,count] = sscanf(line,'%d %d %d %d %d',5);
    assert(count >= 2);
    % add edge
    % .tgf format is 1-indexed
    E = [E;e(1:2)'];
    if count >= 3 && e(3)
      BE = [BE;e(1:2)'];
    end
    % bone edges trump pseudo-edges
    if count >= 4 && e(4) && ~e(3)
      PE = [PE;e(1:2)'];
    end
    % bone edges trump cage edges
    if count >= 5 && e(5) && ~e(3)
      CE = [CE;e(1:2)'];
    end
    line = fscanf(fp,' %[^\n]s');
  end

  % point handles are points not attached to a bone
  P = 1:size(V,1);
  P = P(~ismember(P,BE(:)));
  fclose(fp);
end

author: alecjacobson

repo: gptoolbox


text:
function timeSliderStateChanged(src,evt, lvdData, handles, app)
    s = dbstack;
    matches = strfind({s.name},'timeSliderStateChanged');
    matches = cell2mat(matches(2:end));
    tf = any(matches == 1);

    if(not(tf))
        try
            time = evt.Value;
        catch
            time = src.Value;
        end
        hAx = handles.dispAxes;

        markerTrajData = lvdData.viewSettings.selViewProfile.markerTrajData;
        markerBodyData = lvdData.viewSettings.selViewProfile.markerBodyData;
        markerBodyAxesData = lvdData.viewSettings.selViewProfile.markerTrajAxesData;
        markerGrdObjData = lvdData.viewSettings.selViewProfile.markerGrdObjData;
        centralBodyData = lvdData.viewSettings.selViewProfile.centralBodyData;
        pointData = lvdData.viewSettings.selViewProfile.pointData;
        vectorData = lvdData.viewSettings.selViewProfile.vectorData;
        refFrameData = lvdData.viewSettings.selViewProfile.refFrameData;
        angleData = lvdData.viewSettings.selViewProfile.angleData;
        planeData = lvdData.viewSettings.selViewProfile.planeData;
        sensorData = lvdData.viewSettings.selViewProfile.sensorData;
        sensorTgtData = lvdData.viewSettings.selViewProfile.sensorTgtData;

        notify(app, 'GenericStatusLabelUpdate', GenericStatusLabelUpdate('Drawing Scene...'));

%         notify(app, 'GenericStatusLabelUpdate', GenericStatusLabelUpdate('Drawing Vehicle Trajectory...'));
        markerTrajData.plotBodyMarkerAtTime(time, hAx);

%         notify(app, 'GenericStatusLabelUpdate', GenericStatusLabelUpdate('Drawing Celestial Bodies...'));
        for(i=1:length(markerBodyData)) %#ok<*NO4LP>
            markerBodyData(i).plotBodyMarkerAtTime(time, hAx);
        end

%         notify(app, 'GenericStatusLabelUpdate', GenericStatusLabelUpdate('Drawing Vehicle Body Axes...'));
        markerBodyAxesData.plotBodyAxesAtTime(time, hAx);

%         notify(app, 'GenericStatusLabelUpdate', GenericStatusLabelUpdate('Drawing Ground Objects...'));
        for(i=1:length(markerGrdObjData))
            markerGrdObjData(i).plotBodyMarkerAtTime(time, hAx);
        end

%         notify(app, 'GenericStatusLabelUpdate', GenericStatusLabelUpdate('Drawing Geometric Reference Frames...'));
        for(i=1:length(refFrameData))
            refFrameData(i).plotRefFrameAtTime(time, hAx);
        end

%         notify(app, 'GenericStatusLabelUpdate', GenericStatusLabelUpdate('Drawing Geometric Vectors...'));
        for(i=1:length(vectorData))
            vectorData(i).plotVectorAtTime(time, hAx);
        end

%         notify(app, 'GenericStatusLabelUpdate', GenericStatusLabelUpdate('Drawing Geometric Angles...'));
        for(i=1:length(angleData))
            angleData(i).plotAngleAtTime(time, hAx);
        end

%         notify(app, 'GenericStatusLabelUpdate', GenericStatusLabelUpdate('Drawing Geometric Planes...'));
        for(i=1:length(planeData))
            planeData(i).plotPlaneAtTime(time, hAx);
        end

%         notify(app, 'GenericStatusLabelUpdate', GenericStatusLabelUpdate('Drawing Geometric Points...'));
        for(i=1:length(pointData))
            pointData(i).plotPointAtTime(time, hAx);
        end

%         notify(app, 'GenericStatusLabelUpdate', GenericStatusLabelUpdate('Drawing Sensors...'));
        sensorTgtResults = SensorTargetResults.empty(1,0);
        for(i=1:length(sensorData))
            subResults = sensorData(i).plotSensorAtTime(time, hAx);
            sensorTgtResults = horzcat(sensorTgtResults, subResults(:)'); %#ok<AGROW>
        end

        if(numel(sensorTgtResults) > 0)
%             notify(app, 'GenericStatusLabelUpdate', GenericStatusLabelUpdate('Drawing Sensor Targets...'));

            mergedResults = SensorTargetResults.mergeResults(sensorTgtResults);
            for(i=1:length(sensorTgtData))
                sensorTgtData(i).plotTargetResults(mergedResults, hAx);
            end
        end

        centralBodyData.setCentralBodyRotation(time);

%         notify(app, 'GenericStatusLabelUpdate', GenericStatusLabelUpdate('Creating Light Source...'));
        lvdData.viewSettings.selViewProfile.updateLightPosition(time);

        [year, day, hour, minute, sec] = convertSec2YearDayHrMnSec(time);
        epochStr = formDateStr(year, day, hour, minute, sec);

        [~, timeEvtsListboxStrs] = lvdData.script.getAllEvtsThatOccurAtTime(time);
        for(i=1:length(timeEvtsListboxStrs))
            timeEvtsListboxStrs(i) = sprintf('Event: %s',timeEvtsListboxStrs(i));
        end
        evtsStr = strjoin(timeEvtsListboxStrs,'\n');

        if(not(isempty(timeEvtsListboxStrs)))
            dashes = repmat('-',1,max(strlength(timeEvtsListboxStrs)));
        else
            dashes = '';
        end

        tooltipStr = sprintf('UT = %0.3f sec\n%s\n%s', time, dashes, evtsStr);

        handles.timeSliderValueLabel.String = epochStr;
        handles.timeSliderValueLabel.TooltipString = tooltipStr;

        setappdata(app.DispAxesTimeSlider,'lastTime',time);
        drawnow limitrate;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function EEG = smooth_prob(EEG, smoothlen)
if nargin<2
    smoothlen = 2;
end

llt_smooth = zeros(size(EEG.etc.amica.LLt));
smoothwnd = hann(smoothlen*EEG.srate); smoothwnd = smoothwnd/sum(smoothwnd);
for i = 1:size(llt_smooth,1)
    llt_smooth(i,:) = filtfilt(smoothwnd,1,EEG.etc.amica.LLt(i,:)')';
end
for i = 1:size(llt_smooth,1)
    EEG.etc.amica.LLt_smooth(i,:) = llt_smooth(i,:);
end

v_smooth = zeros(size(llt_smooth));
for m = 1:size(v_smooth,1)-1
    v_smooth(m,:) = 1./sum(exp(bsxfun(@minus,llt_smooth,llt_smooth(m,:))),1);
end
v_smooth(m+1,:) = ones(1,size(v_smooth,2)) - sum(v_smooth(1:m,:));
for i = 1:size(v_smooth,1)
    EEG.etc.amica.v_smooth(i,:) = v_smooth(i,:);
end

if ~isempty(EEG.etc.amica.v_cropped)
    EEG.etc.amica.v_smooth_cropped = cropvector(EEG.etc.amica.v_smooth,EEG.etc.amica.crop_start,EEG.etc.amica.crop_end);
end

author: PatternRecognition

repo: OpenBMI


text:
function CPD = set_params(CPD, varargin)
% SET_PARAMS Set the parameters (fields) for a tabular_decision_node object
% CPD = set_params(CPD, name/value pairs)
%
% The following optional arguments can be specified in the form of name/value pairs:
%
% policy - the table containing the policy
%
% e.g., CPD = set_params(CPD, 'policy', T)

args = varargin;
nargs = length(args);
for i=1:2:nargs
  switch args{i},
   case 'policy',   CPD.CPT = args{i+1};
   otherwise,
    error(['invalid argument name ' args{i}]);
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% IMAGE
% See also
%
% Display:
%   clusterMontage    - Used for visualization of clusters of images and videos.
%   im                - Function for displaying grayscale images.
%   filmStrip         - Used to display R stacks of T images as a "filmstrip".
%   makeGif           - Writes a matlab movie to an animated GIF.
%   montage2          - Used to display collections of images and videos.
%   movieToImages     - Creates a stack of images from a matlab movie M.
%   playMovie         - Shows/makes an/several movie(s) from an image sequence.
%
% Histograms:
%   assignToBins      - Quantizes A according to values in edges.
%   histc2            - Multidimensional histogram count with weighted values.
%   histcImWin        - Calculates local histograms at every point in an image I.
%   histcImLoc        - Creates a series of locally position dependent histograms.
%   histMontage       - Used to display multiple 1D histograms.
%
% Generalized correlation:
%   normxcorrn        - Normalized n-dimensional cross-correlation.
%   xcorrn            - n-dimensional cross-correlation.  Generalized version of xcorr2.
%   xeucn             - n-dimensional euclidean distance between each window in A and template T.
%
% Image deformation:
%   imagesAlign       - Fast and robust estimation of homography relating two images.
%   imNormalize       - Various ways to normalize a (multidimensional) image.
%   imShrink          - Used to shrink a multidimensional array I by integer amount.
%   imtransform2      - Applies a linear or nonlinear transformation to an image I.
%   jitterImage       - Creates multiple, slightly jittered versions of an image.
%   opticalFlow       - Coarse-to-fine optical flow using Lucas&Kanade or Horn&Schunck.
%   textureMap        - Maps texture in I according to rsDst and csDst.
%
% Generalized nonmaximal suppression:
%   nonMaxSupr        - Applies nonmaximal suppression on an image of arbitrary dimension.
%   nonMaxSuprList    - Applies nonmaximal suppression to a list.
%   nonMaxSuprWin     - Nonmaximal suppression of values outside of a given window.
%
% Seq files:
%   seqIo             - Utilities for reading and writing seq files.
%   seqReaderPlugin   - Plugin for seqIo and videoIO to allow reading of seq files.
%   seqWriterPlugin   - Plugin for seqIo and videoIO to allow writing of seq files.
%   seqPlayer         - Simple GUI to play seq files.
%
% Behavior annotation for seq files:
%   behaviorAnnotator - Caltech Behavior Annotator.
%   behaviorData      - Retrieve and manipulate behavior annotation of a video.
%
% Binary mask creation:
%   maskCircle        - Creates an image of a 'pie slice' of a circle.
%   maskEllipse       - Creates a binary image of an ellipse.
%   maskGaussians     - Divides a volume into softly overlapping gaussian windows.
%   maskSphere        - Creates an 'image' of a n-dimensional hypersphere.
%
% Linear filtering:
%   convnFast         - Fast convolution, replacement for both conv2 and convn.
%   gaussSmooth       - Applies Gaussian smoothing to a (multidimensional) image.
%   localSum          - Fast routine for box filtering.
%
% Miscellaneous:
%   imMlGauss         - Calculates max likelihood params of Gaussian that gave rise to image G.
%   imrectLite        - A 'lite' version of imrect [OBSOLETE: use imrectRot].
%   imRectRot         - Create a draggable, resizable, rotatable rectangle or ellipse.
%   imwrite2          - Similar to imwrite, except follows a strict naming convention.
%   kernelTracker     - Kernel Tracker from Comaniciu, Ramesh and Meer PAMI 2003.

author: zouchuhang

repo: LayoutNet


text:
function value = s_eqi ( s1, s2 )

%*****************************************************************************80
%
%% S_EQI is a case insensitive comparison of two strings for equality.
%
%  Example:
%
%    S_EQI ( 'Anjana', 'ANJANA' ) is TRUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, S2, the strings to compare.
%
%    Output, logical VALUE, is TRUE if the strings are equal.
%
  FALSE = 0;
  TRUE = 1;

  len1 = length ( s1 );
  len2 = length ( s2 );
  lenc = min ( len1, len2 );

  value = FALSE;

  for i = 1 : lenc

    c1 = ch_cap ( s1(i) );
    c2 = ch_cap ( s2(i) );

    if ( c1 ~= c2 )
      value = FALSE;
      return
    end

  end

  for i = lenc + 1 : len1
    if ( s1(i) ~= ' ' )
      value = FALSE;
      return
    end
  end

  for i = lenc + 1 : len2
    if ( s2(i) ~= ' ' )
      value = FALSE;
      return
    end
  end

  value = TRUE;

  return
end


text:
function tetrahedron_unit_monomial_test ( degree_max )

%*****************************************************************************80
%
%% TETRAHEDRON_UNIT_MONOMIAL_TEST tests TETRAHEDRON_UNIT_MONOMIAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DEGREE_MAX, the maximum total degree of the
%    monomials to check.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TETRAHEDRON_UNIT_MONOMIAL_TEST\n' );
  fprintf ( 1, '  For the unit tetrahedron,\n' );
  fprintf ( 1, '  TETRAHEDRON_UNIT_MONOMIAL returns the exact value of the\n' );
  fprintf ( 1, '  integral of X^ALPHA Y^BETA Z^GAMMA\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Volume = %f\n', tetrahedron_unit_volume ( ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     ALPHA      BETA     GAMMA      INTEGRAL\n' );
  fprintf ( 1, '\n' );

  for alpha = 0 : degree_max
    expon(1) = alpha;
    for beta = 0 : degree_max - alpha
      expon(2) = beta;
      for gamma = 0 : degree_max - alpha - beta
        expon(3) = gamma;
        value = tetrahedron_unit_monomial ( expon );
        fprintf ( 1, '  %8d  %8d  %8d  %14.6e\n', expon(1:3), value );
      end
    end
  end

  return
end


text:
function clusterComp = gacZetaEntropy (subIminuszW)
%% Compute structural complexity from the subpart of the weighted adjacency matrix
% Input:
%   - subIminuszW: the subpart of (I - z*P)
% Output:
%	- clusterComp - strucutral complexity of a cluster.
% by Wei Zhang (wzhang009 at gmail.com), June, 8, 2011

clusterComp = sum(log(real(diag(inv(subIminuszW))))) / size(subIminuszW,1); % x_c in Thm. 2

end

author: jwyang


text:
function toms446_test02 ( )

%*****************************************************************************80
%
%% TOMS446_TEST02 tests MULTPLY, which multiplies two Chebyshev series.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2011
%
%  Author:
%
%    John Burkardt
%
  nf = 5;
  npl = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TOMS446_TEST02\n' );
  fprintf ( 1, '  Test MLTPLY, which computes the\n' );
  fprintf ( 1, '  product of two Chebyshev series.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Multiply series for SIN(X) and COS(X)\n' );
  fprintf ( 1, '  and compare with series for 1/2*SIN(2X).\n' );

  x = cheby ( nf, npl, @functn );

  for i = 1 : npl
    x1(i) = x(i,1);
    x2(i) = x(i,2);
    x(i,3) = 0.5 * x(i,3);
  end

  x3 = mltply ( x1, x2, npl );

  fprintf ( 1, '\n' );
  fprintf ( 1, '        Sin(x)      Cos(x)   1/2*Sin(2x)     RESULT\n' );
  fprintf ( 1, '\n' );

  for i = 1 : npl
    fprintf ( 1, '  %10.4f  %10.4f  %10.4f  %10.4f\n', x(i,1:3), x3(i) );
  end

  return
end


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

%*****************************************************************************80
%
%% JED_TO_YMDF_ISLAMIC_B converts a JED to an Islamic B YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 July 2012
%
%  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 + 7664;

  y_prime = floor ( ( 30 * j_prime + 15 ) / 10631 );
  t_prime = floor ( mod ( 30 * j_prime + 15, 10631 ) / 30 );
  m_prime = floor ( ( 100 * t_prime + 10 ) / 2951 );
  d_prime = floor ( mod ( 100 * t_prime + 10, 2951 ) / 100 );
%
%  Convert the computational date to a calendar date.
%
  d = d_prime + 1;
  m = mod ( m_prime, 12 ) + 1;
  y = y_prime - 5519 + floor ( ( 12 - m ) / 12 );

  return
end