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);


text:
function exp = exp_releasehold(exp)
% Peel off a layer of hold expressions from some expression.
% [Out-Expression] = exp_releasehold(In-Expression)
%
% In:
%   In-Expression      : some expression
%
% Out:
%   Out-Expression     : In-Expression with the first level of Hold's peeled off
%
% Notes:
%   Descends into cell arrays.
%
% See also:
%   exp_hold
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-05-19

if ~exp_beginfun('symbolic') return; end

exp = utl_releasehold(exp);

exp_endfun;

repo: REST


text:
function hP = plugin_board()
% Zynq Platform PCore
% Use Plugin API to create board plugin object

%   Copyright 2015 The MathWorks, Inc.

% Call the common board definition function
hP = AnalogDevices.fmcomms5.common.plugin_board('ZC706');

author: analogdevicesinc


text:
function ellipsoid_monte_carlo_test01 ( )

%*****************************************************************************80
%
%% ELLIPSOID_MONTE_CARLO_TEST01 uses ELLIPSOID_SAMPLE on a 2D ellipse centered at (0,0).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2014
%
%  Author:
%
%    John Burkardt
%
  m = 2;

  a = [ ...
    9.0, 1.0;
    1.0, 4.0 ]';
  e_test = [ ...
    0, 0;
    1, 0;
    0, 1;
    2, 0;
    1, 1;
    0, 2;
    3, 0 ]';
  r = 2.0;
  v = [ 0.0, 0.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ELLIPSOID_MONTE_CARLO_TEST01\n' );
  fprintf ( 1, '  Use ELLIPSOID_SAMPLE to estimate integrals\n' );
  fprintf ( 1, '  in a 2D ellipse x'' * A * x <= r^2.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Ellipsoid radius R = %g\n', r );
  r8vec_print ( m, v, '  Ellipsoid center V:' );
  r8mat_print ( m, m, a, '  Ellipsoid matrix A:' );

  volume = ellipsoid_volume ( m, a, v, r );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Ellipsoid volume = %g\n', volume );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '         N        1              X               Y  ' );
  fprintf ( 1, '             X^2               XY             Y^2             X^3\n' );
  fprintf ( 1, '\n' );

  n = 1;

  while ( n <= 65536 )

    [ x, seed ] = ellipsoid_sample ( m, n, a, v, r, seed );

    fprintf ( 1, '  %8d', n );

    for j = 1 : 7

      e(1:m) = e_test(1:m,j);

      value = monomial_value ( m, n, e, x );

      result = volume * sum ( value(1:n) ) / n;

      fprintf ( 1, '  %14.6g', result );

    end

    fprintf ( 1, '\n' );

    n = 2 * n;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%Implements an example MAP decoder, hard-input soft-output. Uses a
%brute-force algorithm, which is useful for demonstration/learning
%purposes.
%
% Copyright Colin O'Flynn, 2011. All rights reserved.
% http://www.newae.com
%
% Redistribution and use in source and binary forms, with or without modification, are
% permitted provided that the following conditions are met:
%
%    1. Redistributions of source code must retain the above copyright notice, this list of
%       conditions and the following disclaimer.
%
%    2. Redistributions in binary form must reproduce the above copyright notice, this list
%       of conditions and the following disclaimer in the documentation and/or other materials
%       provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY COLIN O'FLYNN ''AS IS'' AND ANY EXPRESS OR IMPLIED
% WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL COLIN O'FLYNN OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
% ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

function [llrs, codeword] = brute_force_map(feedback, feedforward, input, pberr, nbits)

    if nbits > 16
        error('nbits too high')
    end

    %Generate every possible information sequence
    D = [0:2^nbits - 1];
    B = dec2bin(D);

    allValidInputs = zeros(2^nbits, nbits);

    %Convert from string to matrix
    for i=1:nbits
        allValidInputs(:,i) = str2num(B(:,i));
    end

    len = (nbits+3)*2;
    all_codewords = zeros(2^nbits, len);

    %Generate every possible codeword
    for i=1:2^nbits
        codeword = rsc_encode([feedback; feedforward], allValidInputs(i,:), 1);
        all_codewords(i,:) =  reshape(codeword, 1, len);
    end

    pcodeword = zeros(2^nbits, 1);

    %For every possible codeword & ours: find out Pcodeword
    for i=1:2^nbits
        bitsInDiff = sum(abs(input - all_codewords(i,:)));
        bitsOK = len - bitsInDiff;
        %Find APP Pr{X | Y}
        % X = Codeword that was transmitted
        % Y = Codeword that was receieved
        pcodeword(i) = pberr^bitsInDiff * (1-pberr)^bitsOK;
    end

    %Limited valid Tx codewords, so normalize probability to add up to 1.0
    pcodeword = pcodeword ./ (sum(pcodeword));

    %Find resulting maximum likilhood codeword
    [Pmax, Imax] = max(pcodeword);

    %The suggested codeword
    codeword = all_codewords(Imax, :);

    %Reshape to make division between systematic & parity obvious
    %Now looks like: [systematic ; parity]
    codeword = reshape(codeword, 2, len/2);

    %Calculate individual probability of error
    p0Systematic = zeros(1, nbits);
    for bitindex=1:nbits
        psum = 0;

        %Sum over all codewords
        for i=1:2^nbits
            %Reshape codeword to easily get systematic part out
            codewords_reshaped = reshape(all_codewords(i,:), 2, len/2);

            %If codeword has bit n as zero, count it in summation
            if codewords_reshaped(1, bitindex) == 0
                psum = psum + pcodeword(i);
            end
        end

        %Save total probability of ALL codewords with bit 'n' is zero
        p0Systematic(bitindex) = psum;
    end

    %Find probability any given bit is ONE
    p1Systematic = 1.0 - p0Systematic;

    %From P1 & P0 calculate LLR
    llrs = log(p1Systematic ./ p0Systematic);

author: Sable

repo: mcbench-benchmarks


text:
function SO3VF = sqrt(SO3VF, varargin)
% square root of a SO3VectorField
%
% Syntax
%   SO3VF = sqrt(SO3VF)
%
% Input
%  SO3VF - @SO3VectorField
%
% Output
%  SO3VF - @SO3VectorField
%

SO3VF = SO3VF.^(1/2);

end

author: mtex-toolbox

repo: mtex

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


text:
function p = addMultipliedInequalityCuts(p)

newRows = [];
for i =  p.linears
    if p.lb(i) >= -1000 | p.ub(i) <= 1000 % Avoid numerically tricky stuff
        for j = 1+p.K.f:p.K.f+p.K.l
            row = p.F_struc(j,:);
            used = find(row(2:end));
            if ~any(p.variabletype(used))
                % Linear constraint
                ok = 1;
                monom = zeros(1,length(used));
                for k = 1:length(used)
                    find1 = findrows(p.bilinears(:,2:end),[i used(k)]);
                    find2 = findrows(p.bilinears(:,2:end),[used(k) i]);
                    if (isempty(find1) & isempty(find2))
                        ok = 0;
                        break
                    else
                        monom(k) = unique([find1 find2]);
                    end
                end
                if ok
                    monom = p.bilinears(monom,1);
                    row0 = row;
                    if p.lb(i) >= -1000
                        row(1+monom) = row(1+used);
                        row(1+used) = 0;
                        row(1+i) = row(1);
                        row(1)=0;
                        newRows = [newRows;row-p.lb(i)*row0];
                    end
                    if p.ub(i) <= 1000
                        row = row0;
                        row(1+monom) = row(1+used);
                        row(1+used) = 0;
                        row(1+i) = row(1);
                        row(1)=0;
                        newRows = [newRows;-row+p.ub(i)*row0];
                    end
                end
            end
        end
    end
end
if ~isempty(newRows)
    p.F_struc = [p.F_struc;newRows];
    p.K.l = p.K.l + size(newRows,1);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [grid3d, ge, osc] = read_grid(filenamebase)

chkarg(istypesizeof(filenamebase, 'char', [1 0]), '"filenamebase" should be string.');

inputfile = ['./', filenamebase, '.h5'];

lambda = h5read(inputfile, '/lambda');
L0 = h5read(inputfile, '/L0');
Npml = h5read(inputfile, '/Npml').';
ge = h5read(inputfile, '/ge');
ge = GT.elems(ge+1);
bc = BC.elems(h5read(inputfile, '/bc').' + 1);
bc = bc(:, Sign.n).';

lprim = cell(1, Axis.count);
for w = Axis.elems
	lprim{w} = h5read(inputfile, ['/', char(w), '_prim']).';
end

unit = PhysUnit(L0);
osc = Oscillation(lambda, unit);
grid3d = Grid3d(unit, lprim, Npml, bc);

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
%% BANG_WHO_TEST tests the BANG facility for issuing operating system commands.
%
%  Discussion:
%
%    WHO is a legal command on UNIX systems, and returns a list of the
%    users currently logged in.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 July 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BANG_WHO_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Demonstrate a simple use of the ! command,\n' );
  fprintf ( 1, '  which allows MATLAB to invoke operating system commands.\n' );

  ! who

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


text:
function p = createMultiframeTestphan()

% load the k-space of a predifined multiecho testphantom from file
basePath = fileparts(mfilename('fullpath'));
dat = load(fullfile(basePath, 'testPhan64.mat'));
p = dat.p;

% create a movie with different (simulated) shifts of the phantom
NF = 30; % number of frames
dim = size(p,1);
P = zeros([size(p),NF]);

x = linspace(0,1,dim);
[X,Y] = meshgrid(x,x');
shiftMax = dim * 2*pi;
shiftInc = shiftMax / NF;

shift = 0;
for i = 1 : NF
    shiftPhX = exp(1i * -shift * X);
    shiftPhY = exp(1i * shift * Y);
    shift = shift + shiftInc;
    currP = ftimes(p,shiftPhX);
    currP = ftimes(currP,shiftPhY);
    P(:,:,:,i) = currP;
end

% do a Fourier reconstruction of the data
P = asDataClass.mrIfft(P);

% simulate 4 different coil profiles
coils = cat(5,X,Y,1-X, 1-Y);
P = ftimes(P,coils);

% simulate phase perturbations
phase = exp(1i *2 * pi * (1.5 - coils));
P = ftimes(P,phase);

% create pseudo k-space data from the fft of the images
p = asDataClass.mrFft(P);
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function proj2d = compute_one_projection_method_2(xs,ys,zs,data3d,psrc,pcdet,su,sv,nu,nv)
%
%	proj2d = compute_one_projection_method_2(xs,ys,zs,data3d,psrc,pcdet,su,sv,nu,nv)
%
%	Method 2: compute the project based on the stack of planes perpendicular to the src-detector line
%	This method should be faster than the line-integral method
%
%   Deshan Yang, PhD
%	Department of radiation oncology
%	Washington University in Saint Louis
%   01/16/2011, Saint Louis, MO, USA
%

corner_points = [...
	xs(1) ys(1) zs(1);...
	xs(1) ys(1) zs(end);...
	xs(1) ys(end) zs(1);...
	xs(1) ys(end) zs(end);...
	xs(end) ys(1) zs(1);...
	xs(end) ys(1) zs(end);...
	xs(end) ys(end) zs(1);...
	xs(end) ys(end) zs(end)...
	];

corner_points_t = zeros(8,1);
for k =1:8
	p = corner_points(k,:);
	[q,corner_points_t(k)] = project_1_point_to_a_line(psrc,pcdet,p);
end

t_min = min(corner_points_t);
t_max = max(corner_points_t);

dist = norm(psrc-pcdet);

N = ceil(dist*(t_max-t_min));	% N planes
dt = (t_max-t_min)/N;
ts = t_min:dt:t_max;

proj2d = zeros(nv,nu,'single');
vecu = [0 0 1];

vecn = psrc-pcdet;
vecv = cross(vecn,vecu);
vecv = vecv/norm(vecv);

us = ((-nu/2+0.5):1:(nu/2-0.5))*su/nu;
vs = ((-nv/2+0.5):1:(nv/2-0.5))*sv/nv;
[uu,vv] = meshgrid(us,vs);

for T = 1:N
	vecu_t = vecu*ts(T);
	vecv_t = vecv*ts(T);
	pcdet_t = psrc+ts(T)*(pcdet-psrc);

	vdet_t_xs = pcdet_t(1) + vecu_t(1)*uu + vecv_t(1)*vv;
	vdet_t_ys = pcdet_t(2) + vecu_t(2)*uu + vecv_t(2)*vv;
	vdet_t_zs = pcdet_t(3) + vecu_t(3)*uu + vecv_t(3)*vv;

	proj2d = proj2d + interp3(xs,ys,zs,data3d,vdet_t_xs,vdet_t_ys,vdet_t_zs,'linear',0)*dt*dist;
end

% adjust the value according to the projection path length
vdet_t_xs = pcdet(1) + vecu(1)*uu + vecv(1)*vv;
vdet_t_ys = pcdet(2) + vecu(2)*uu + vecv(2)*vv;
vdet_t_zs = pcdet(3) + vecu(3)*uu + vecv(3)*vv;

dists = sqrt((vdet_t_xs-psrc(1)).^2+(vdet_t_ys-psrc(2)).^2+(vdet_t_zs-psrc(3)).^2);
dists = dists / min(dists(:));
proj2d = proj2d .* dists;


text:
function [res, isCancel] = java_dialog( msgType, msg, msgTitle, jParent, varargin )
% JAVA_DIALOG: Display a Java-based dialog.
%
% USAGE:  [res, isCancel] = java_dialog( msgType, msg )
%         [res, isCancel] = java_dialog( msgType, msg, msgTitle )
%         [res, isCancel] = java_dialog( msgType, msg, msgTitle, jParent )
%         [res, isCancel] = java_dialog( msgType, msg, msgTitle, jParent, OPTIONS )
%
% INPUT:
%    - msgType  : String that defines the icon and controls displayed in the dialog box. Possible values:
%                 {'error', 'warning', 'msgbox', 'confirm', 'question', 'input', 'checkbox', 'radio', 'combo', 'hotkey', 'color'}
%    - msg      : Message in the dialog box
%                 (can be a cell array of strings for 'input' message type)
%    - msgTitle : Title of the dialog box
%    - jParent  : Handle to the parent JFrame.
%                 If no parent, or unknown, set to [].
%    - OPTIONS  : Depends on the dialog type
%                  - 'question' : OPTIONS = buttonList
%                                 OPTIONS = buttonList, defaultButton
%                  - 'combo'    : OPTIONS = buttonList
%                                 OPTIONS = buttonList, defaultButton
%                  - 'input'    : OPTIONS = defaultVals
%                  - 'checkbox' : OPTIONS = buttonList
%                                 OPTIONS = buttonList, defaultVals
%                  - 'radio'    : OPTIONS = buttonList
%                                 OPTIONS = buttonList, defaultInd
%                  - 'color'    : No options
% OUTPUT:
%    - res      : Depends on the dialog type
%                  - 'error', 'warning', 'msgbox', 'confirm':
%                        1 if user validated the dialog, 0 else.
%                  - 'question':
%                        'Yes' if user clicked yes button, 'No' for no button,
%                        [] if user closed the dialog.
%                  - 'input':
%                        str or cell array of str corresponding to selected input
%                        [] if dialog canceled
%                  - 'radio':
%                        [] if user closed dialog or cancelled, else index of choice.
%                  - 'combo':
%                        value of choice (empty if user cancelled or closed dialog)
%                  - 'hotkey':
%                        [] if user closed dialog or cancelled, else chosen hotkey
%    - isCancel : 1 if the user pressed the cancel button or closed the
%                 window (0 else). Always 0 for dialog types 'error', 'warning' and 'msgbox'

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2008-2019

global GlobalData;

% Java imports
import org.brainstorm.dialogs.*;

% Get brainstorm frame
jBstFrame = bst_get('BstFrame');

% Parse inputs
if (nargin < 1)
    error('Invalid call to java_dialog()');
end
if (nargin < 2)
    msg = '';
end
if (nargin < 3)
    msgTitle = '';
end
if (nargin < 4) || isempty(jParent)
    jParent = jBstFrame;
end
isCancel = 1;

% Hide progress bar
isProgress = bst_progress('isvisible');
if isProgress
    bst_progress('hide');
end

% Headless dialog displayed in command line
if bst_get('GuiLevel') < 0
    switch(lower(msgType))
        case 'error'
            disp(['Error: ' msg]);
            isCancel = 0;
            res = 1;
        case 'warning'
            disp(['Warning: ' msg]);
            isCancel = 0;
            res = 1;
        case 'msgbox'
            disp(msg);
            isCancel = 0;
            res = 1;
        case 'confirm'
            response = lower(strtrim(input([msg ' (y/n) '], 's')));
            if ismember(response, {'y', 'yes'})
                isCancel = 0;
                res = 1;
            elseif ismember(response, {'n', 'no'})
                isCancel = 0;
                res = 0;
            else
                isCancel = 1;
                res = 0;
            end
        case 'question'
            if length(varargin) < 1 || isempty(varargin{1})
                answers = {'Yes', 'No'};
            else
                answers = varargin{1};
            end
            nAnswers = length(answers);

            % Display question and possible answers
            disp(msg);
            for iAnswer = 1:nAnswers
                disp([num2str(iAnswer) '. ' answers{iAnswer}]);
            end
            disp(' ');

            % Validate user selection
            response = str2num(input('Select the number of your answer: ', 's'));
            if ~isempty(response) && response >= 1 && response <= nAnswers
                res = answers{response};
                isCancel = 0;
            else
                res = [];
                isCancel = 1;
            end
        otherwise
            error('Unsupported call of java_dialog in headless mode.');
    end
    return;
end

% Get all the modal JDialogs
if ~isempty(jBstFrame)
    jDialogModal = [];
    jDialogAlwaysOnTop = [];
    for i=1:length(GlobalData.Program.GUI.panels)
        panelContainer = get(GlobalData.Program.GUI.panels(i), 'container');
        panelContainer = panelContainer.handle{1};
        if isa(panelContainer, 'javax.swing.JDialog') && panelContainer.isModal()
            % A modal JDialog is found => Set it non non-modal
            jDialogModal = panelContainer;
            jDialogModal.setModal(0);
        end
        if (isa(panelContainer, 'javax.swing.JDialog') || isa(panelContainer, 'javax.swing.JFrame')) && panelContainer.isAlwaysOnTop()
            % An AlwaysOnTop frame is found => Remove always on top attribute
            jDialogAlwaysOnTop = panelContainer;
            jDialogAlwaysOnTop.setAlwaysOnTop(0);
        end
    end
end

% Check that parent frame is visible
if ~isempty(jParent) && ~jParent.isVisible()
    jParent = [];
end

% Show java dialog
switch(lower(msgType))
    case 'error'
        if isempty(msgTitle)
            msgTitle = 'Error';
        end
        java_call('javax.swing.JOptionPane', 'showMessageDialog', 'Ljava.awt.Component;Ljava.lang.Object;Ljava.lang.String;I', jParent, msg, msgTitle, javax.swing.JOptionPane.ERROR_MESSAGE);
        isCancel = 0;
        res = 1;
    case 'errorhelp'
        if isempty(msgTitle)
            msgTitle = 'Error';
        end
        if ~java_call('org.brainstorm.dialogs.MsgServer', 'dlgErrorHelp', 'Ljava.awt.Component;Ljava.lang.String;Ljava.lang.String;[Ljava.lang.String;Ljava.lang.String;', jParent, msg, msgTitle);
            web('https://neuroimage.usc.edu/forums/', '-browser');
        end
        isCancel = 0;
        res = 1;
    case 'warning'
        if isempty(msgTitle)
            msgTitle = 'Warning';
        end
        java_call('javax.swing.JOptionPane', 'showMessageDialog', 'Ljava.awt.Component;Ljava.lang.Object;Ljava.lang.String;I', jParent, msg, msgTitle, javax.swing.JOptionPane.WARNING_MESSAGE);
        isCancel = 0;
        res = 1;
    case 'msgbox'
        if isempty(msgTitle)
            msgTitle = 'Information';
        end
        java_call('javax.swing.JOptionPane', 'showMessageDialog', 'Ljava.awt.Component;Ljava.lang.Object;Ljava.lang.String;I', jParent, msg, msgTitle, javax.swing.JOptionPane.INFORMATION_MESSAGE);
        isCancel = 0;
        res = 1;
    case 'confirm'
        if isempty(msgTitle)
            msgTitle = 'Confirmation';
        end
        reponse = java_call('javax.swing.JOptionPane', 'showConfirmDialog', 'Ljava.awt.Component;Ljava.lang.Object;Ljava.lang.String;I', jParent, msg, msgTitle, javax.swing.JOptionPane.YES_NO_OPTION);
        if reponse == javax.swing.JOptionPane.CLOSED_OPTION
            res = 0;
            isCancel = 1;
        else
            res = (reponse == javax.swing.JOptionPane.OK_OPTION);
            isCancel = 0;
        end
    % OPTIONS: buttonList; buttonDefault
    case 'question'
        % Button list
        if (length(varargin) < 1)
            buttonList = {'Yes','No'};
        elseif iscell(varargin{1})
            buttonList = varargin{1};
        else
            error('Invalid call.');
        end
        % Default button
        if (length(varargin) < 2)
            buttonDefault = buttonList{1};
        elseif ischar(varargin{2})
            buttonDefault = varargin{2};
        else
            error('Invalid call.');
        end
        % Show dialog
        java_res = java_call('org.brainstorm.dialogs.MsgServer', 'dlgQuest', 'Ljava.awt.Component;Ljava.lang.String;Ljava.lang.String;[Ljava.lang.String;Ljava.lang.String;', jParent, msg, msgTitle, buttonList, buttonDefault);
        if isempty(char(java_res))
            res = [];
            isCancel = 1;
        else
            res = char(java_res);
            isCancel = 0;
        end

    % OPTIONS: defaultVals
    case 'input'
        % Default values
        if (length(varargin) < 1)
            if ischar(msg)
                defaultVals = '';
            elseif iscell(msg)
                defaultVals = repmat({''}, 1, length(msg));
            end
        elseif ~((ischar(msg) && ischar(varargin{1})) || (iscell(msg) && iscell(varargin{1}) && (length(varargin{1}) == length(msg))))
            error('First and third argument must have the same type and the same size.');
        else
            defaultVals = varargin{1};
        end
        % === CALL INPUT DIALOG ===
        if iscell(msg)
            java_res = java_call('org.brainstorm.dialogs.MsgServer', 'dlgInput', 'Ljava/awt/Component;[Ljava/lang/String;Ljava/lang/String;[Ljava/lang/String;', jParent, msg, msgTitle, defaultVals);
            if ~isempty(java_res)
                isCancel = 0;
            end
        else
            java_res = java_call('org.brainstorm.dialogs.MsgServer', 'dlgInput', 'Ljava/awt/Component;Ljava/lang/String;Ljava/lang/String;Ljava/lang/String;', jParent, msg, msgTitle, defaultVals);
            if (isa(java_res, 'java.lang.String') || ischar(java_res)) && ~isequal(char(java_res), '#cancel#')
                isCancel = 0;
            end
        end
        % === GET RESULTS ===
        if isCancel
            res = [];
        elseif isempty(java_res)
            res = '';
        elseif ischar(msg)
            res = char(java_res);
            isCancel = 0;
        elseif iscell(msg)
            nbInputs = numel(java_res);
            res = cell(1, nbInputs);
            for i=1:nbInputs
                res{i} = char(java_res(i));
            end
            isCancel = 0;
        end

    % Options: buttonList, defaultVal
    case 'checkbox'
        % List of options
        if iscell(varargin{1})
            buttonList = varargin{1};
        else
            error('Invalid call.');
        end
        % Default values
        if (length(varargin) >= 2)
            defaultVal = varargin{2};
        else
            defaultVal = [];
        end

        % Create a dialog message
        if ~isempty(buttonList)
            jPanel = gui_river([2,2], [3,3,3,3]);
            jCheck = javaArray('javax.swing.JCheckBox', length(buttonList));
            isLongText = any(cellfun(@length, buttonList) > 20);
            for i = 1:length(buttonList)
                % Create check box
                jCheck(i) = javax.swing.JCheckBox(buttonList{i});
                % Default value
                if ~isempty(defaultVal)
                    jCheck(i).setSelected(defaultVal(i));
                end
                % Add: right after of on the next row ?
                if isLongText
                    jPanel.add('br', jCheck(i));
                else
                    if (mod(i-1, 5) == 0)
                        jPanel.add('tab br', jCheck(i));
                    else
                        jPanel.add('tab', jCheck(i));
                    end
                end
            end
        end
        message = javaArray('java.lang.Object',2);
        message(1) = java.lang.String(msg);
        message(2) = jPanel;
        % Show question
        answer = java_call('javax.swing.JOptionPane', 'showConfirmDialog', 'Ljava.awt.Component;Ljava.lang.Object;Ljava.lang.String;I', jParent, message, msgTitle, javax.swing.JOptionPane.OK_CANCEL_OPTION);
        if (answer ~= javax.swing.JOptionPane.OK_OPTION)
            res = [];
        else
            isCancel = 0;
            % Else look for the selected options
            if isempty(defaultVal)
                res = {};
                for i = 1:length(buttonList)
                    if (jCheck(i).isSelected())
                        res{end+1} = buttonList{i};
                    end
                end
            else
                res = [];
                for i = 1:length(buttonList)
                    res(end+1) = double(jCheck(i).isSelected());
                end
            end
        end

    % Options: buttonList, defaultInd
    case 'radio'
        % List of options
        if iscell(varargin{1})
            buttonList = varargin{1};
        else
            error('Invalid call.');
        end
        % Default values
        if (length(varargin) >= 2)
            defaultInd = varargin{2};
        else
            defaultInd = [];
        end

        % Create a dialog message
        if ~isempty(buttonList)
            jPanel = gui_river([2,2], [3,3,3,3]);
            jRadio = javaArray('javax.swing.JRadioButton', length(buttonList));
            jGroup = javax.swing.ButtonGroup;
            for i = 1:length(buttonList)
                % Create check box
                jRadio(i) = javax.swing.JRadioButton(buttonList{i});
                % Add it to button group
                jGroup.add(jRadio(i));
                % Box size
                isLongText = (length(buttonList{i}) > 10);
                if ~isLongText
                    jRadio(i).setPreferredSize(java_scaled('dimension', 80, 20));
                end
                % Default value
                if ~isempty(defaultInd) && (i == defaultInd)
                    jRadio(i).setSelected(1);
                end
                % Add: right after of on the next row ?
                if (mod(i-1, 5) == 0) || isLongText
                    jPanel.add('br', jRadio(i));
                else
                    jPanel.add(jRadio(i));
                end
            end
        end
        message = javaArray('java.lang.Object',2);
        message(1) = java.lang.String(msg);
        message(2) = jPanel;
        % Show question
        answer = java_call('javax.swing.JOptionPane', 'showConfirmDialog', 'Ljava.awt.Component;Ljava.lang.Object;Ljava.lang.String;I', jParent, message, msgTitle, javax.swing.JOptionPane.OK_CANCEL_OPTION);
        if (answer ~= javax.swing.JOptionPane.OK_OPTION)
            res = [];
        else
            isCancel = 0;
            % Else look for the selected option
            for i = 1:length(jRadio)
                if (jRadio(i).isSelected())
                    res = i;
                    break;
                end
            end
        end

    % OPTIONS: buttonList; buttonDefault
    case 'combo'
        % List of options
        if iscell(varargin{1})
            buttonList = varargin{1};
        else
            error('Invalid call.');
        end
        % Default values
        if (length(varargin) >= 2)
            iSelect = find(strcmpi(varargin{2}, buttonList));
        else
            iSelect = [];
        end

        % Create a dialog message
        jCombo = gui_component('ComboBox', [], [], [], {buttonList}, [], [], []);
        jCombo.setMaximumRowCount(30);
        % Default selection
        if ~isempty(iSelect)
            jCombo.setSelectedIndex(iSelect - 1);
        end
        % Message
        message = javaArray('java.lang.Object',2);
        message(1) = java.lang.String(msg);
        message(2) = jCombo;
        % Show question
        answer = java_call('javax.swing.JOptionPane', 'showConfirmDialog', 'Ljava.awt.Component;Ljava.lang.Object;Ljava.lang.String;I', jParent, message, msgTitle, javax.swing.JOptionPane.OK_CANCEL_OPTION);
        if (answer ~= javax.swing.JOptionPane.OK_OPTION)
            res = [];
        else
            isCancel = 0;
            iSelect = jCombo.getSelectedIndex() + 1;
            res = buttonList{iSelect};
        end

    % Dialog for hotkey input from user
    case 'hotkey'
        % Open up hotkey dialog
        fontSize = round(11 * bst_get('InterfaceScaling') / 100);
        dialog = HotkeyDialog(fontSize);
        % If a valid key was selected, return it. Otherwise, cancel
        if dialog.hasKey()
            res = dialog.getKey();
            isCancel = 0;
        else
            res = [];
            isCancel = 1;
        end

    % Color picker
    case 'color'
        % Create color chooser
        if ~isfield(GlobalData.Program, 'ColorChooser') || isempty(GlobalData.Program.ColorChooser)
            jColorChooser = javax.swing.JColorChooser();
            GlobalData.Program.ColorChooser = javax.swing.JColorChooser();
        else
            jColorChooser = GlobalData.Program.ColorChooser;
        end
        % Show question
        answer = java_call('javax.swing.JOptionPane', 'showConfirmDialog', 'Ljava.awt.Component;Ljava.lang.Object;Ljava.lang.String;I', jParent, jColorChooser, msgTitle, javax.swing.JOptionPane.OK_CANCEL_OPTION, javax.swing.JOptionPane.PLAIN_MESSAGE);
        if (answer ~= javax.swing.JOptionPane.OK_OPTION)
            res = [];
        else
            isCancel = 0;
            color = jColorChooser.getColor();
            res = double([color.getRed(), color.getGreen(), color.getBlue()]) ./ 255;
        end
end

% Restore modal panels
if ~isempty(jBstFrame)
    if ~isempty(jDialogModal)
        jDialogModal.setModal(1);
    end
    if ~isempty(jDialogAlwaysOnTop)
        jDialogAlwaysOnTop.setAlwaysOnTop(1);
    end
end

% Restore progress bar
if isProgress
    bst_progress('show');
end

author: brainstorm-tools


text:
% Kernel Normalized Least-Mean-Square algorithm with Coherence Criterion
%
% C. Richard, J.C.M. Bermudez, and P. Honeine, "Online Prediction of Time
% Series Data With Kernels," IEEE Transactions on Signal Processing,
% vol. 57, no. 3, pp. 1058-1067, March 2009,
% http://dx.doi.org/10.1109/TSP.2008.2009895
%
% Remark: memories are initialized empty in this implementation
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

classdef knlms < kernel_adaptive_filter

    properties (GetAccess = 'public', SetAccess = 'private') % parameters
        eta = .5; % step size
        mu0 = .95; % coherence criterion threshold
        eps = 1E-2; % regularization
        kerneltype = 'gauss'; % kernel type
        kernelpar = 1; % kernel parameter
    end

    properties (GetAccess = 'public', SetAccess = 'private') % variables
        dict = []; % dictionary
        modict = []; % modulus of the dictionary elements
        alpha = []; % expansion coefficients
    end

    methods
        function kaf = knlms(parameters) % constructor
            if (nargin > 0) % copy valid parameters
                for fn = fieldnames(parameters)'
                    if ismember(fn,fieldnames(kaf))
                        kaf.(fn{1}) = parameters.(fn{1});
                    end
                end
            end
        end

        function y_est = evaluate(kaf,x) % evaluate the algorithm
            if size(kaf.dict,1)>0
                k = kernel(kaf.dict,x,kaf.kerneltype,kaf.kernelpar);
                y_est = k'*kaf.alpha;
            else
                y_est = zeros(size(x,1),1);
            end
        end

        function train(kaf,x,y) % train the algorithm
            if size(kaf.dict,2)==0 % initialize
                k = kernel(x,x,kaf.kerneltype,kaf.kernelpar);
                kaf.dict = x;
                kaf.modict = sqrt(k);
                kaf.alpha = 0;
            else
                k = kernel(kaf.dict,x,kaf.kerneltype,kaf.kernelpar);
                kx = kernel(x,x,kaf.kerneltype,kaf.kernelpar);
                C = k./(sqrt(kx)*kaf.modict); % coherence
                if (max(C) <= kaf.mu0) % coherence criterion
                    kaf.dict = [kaf.dict; x]; % order increase
                    kaf.modict = [kaf.modict; sqrt(kx)];
                    kaf.alpha = [kaf.alpha; 0]; % reserve spot
                end
            end

            k = kernel(kaf.dict,x,kaf.kerneltype,kaf.kernelpar);
            kaf.alpha = kaf.alpha + ... % update coefficients
                kaf.eta / (kaf.eps + k'*k) * (y - k'*kaf.alpha) * k;
        end
    end
end

repo: kafbox


text:
function dz = dynamics(z,u,param)
% dz = dynamics(z,u,param)
%
% Computes the first-order form of the dynamics for the combined chain
% integrator and pendulum system
%
% INPUTS:
%   z = [x;v1;v2];
%   u = [u1;u2];
%
% OUTPUTS:
%   dz = dz/dt
%

x = z(1,:);
v1 = z(2,:);
% v2 = z(3,:);   %Unused
u1 = u(1,:);
u2 = u(2,:);

% Pendulum physics
dv1 = pendulum(x,v1,u1,param);

% Integrator chain physics:
dx = v1;
dv2 = u2;

% Combine:
dz = [dx;dv1;dv2];

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function [ point_coord, face_order, face_point ] = soccer_shape_3d ( ...
  point_num, face_num, face_order_max )

%*****************************************************************************80
%
%% SOCCER_SHAPE_3D describes a truncated icosahedron in 3D.
%
%  Discussion:
%
%    The shape is a truncated icosahedron, which is the design used
%    on a soccer ball.  There are 12 pentagons and 20 hexagons.
%
%    Call SOCCER_SIZE_3D to get the values of POINT_NUM, FACE_NUM, and
%    FACE_ORDER_MAX, so you can allocate space for the arrays.
%
%    For each face, the face list must be of length FACE_ORDER_MAX.
%    In cases where a face is of lower than maximum order (the
%    12 pentagons, in this case), the extra entries are listed as
%    "-1".
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    http://mathworld.wolfram.com/TruncatedIcosahedron.html
%
%  Parameters:
%
%    Input, integer POINT_NUM, the number of points in the shape (60).
%
%    Input, integer FACE_NUM, the number of faces in the shape (32).
%
%    Input, integer FACE_ORDER_MAX, the maximum order of any face (6).
%
%    Output, real POINT_COORD(3,POINT_NUM), the vertices.
%
%    Output, integer FACE_ORDER(FACE_NUM), the number of vertices per face.
%
%    Output, integer FACE_POINT(FACE_ORDER_MAX,FACE_NUM); FACE_POINT(I,J)
%    contains the index of the I-th point in the J-th face.  The
%    points are listed in the counter-clockwise direction defined
%    by the outward normal at the face.
%
  dim_num = 3;
%
%  Set the point coordinates.
%
  point_coord(1:dim_num,1:point_num) = [ ...
       -1.00714,    0.153552,   0.067258; ...
       -0.960284,   0.0848813, -0.33629; ...
       -0.95172,   -0.153552,   0.33629; ...
       -0.860021,   0.529326,   0.150394; ...
       -0.858,     -0.290893,  -0.470806; ...
       -0.849436,  -0.529326,   0.201774; ...
       -0.802576,  -0.597996,  -0.201774; ...
       -0.7842,     0.418215,  -0.502561; ...
       -0.749174,  -0.0848813,  0.688458; ...
       -0.722234,   0.692896,  -0.201774; ...
       -0.657475,   0.597996,   0.502561; ...
       -0.602051,   0.290893,   0.771593; ...
       -0.583675,  -0.692896,   0.470806; ...
       -0.579632,  -0.333333,  -0.771593; ...
       -0.52171,   -0.418215,   0.771593; ...
       -0.505832,   0.375774,  -0.803348; ...
       -0.489955,  -0.830237,  -0.33629; ...
       -0.403548,   0.,        -0.937864; ...
       -0.381901,   0.925138,  -0.201774; ...
       -0.352168,  -0.666667,  -0.688458; ...
       -0.317142,   0.830237,   0.502561; ...
       -0.271054,  -0.925138,   0.33629; ...
       -0.227464,   0.333333,   0.937864; ...
       -0.224193,  -0.993808,  -0.067258; ...
       -0.179355,   0.993808,   0.150394; ...
       -0.165499,   0.608015,  -0.803348; ...
       -0.147123,  -0.375774,   0.937864; ...
       -0.103533,   0.882697,  -0.502561; ...
       -0.0513806,  0.666667,   0.771593; ...
        0.0000000,  0.,         1.021; ...
        0.0000000,  0.,        -1.021; ...
        0.0513806, -0.666667,  -0.771593; ...
        0.103533,  -0.882697,   0.502561; ...
        0.147123,   0.375774,  -0.937864; ...
        0.165499,  -0.608015,   0.803348; ...
        0.179355,  -0.993808,  -0.150394; ...
        0.224193,   0.993808,   0.067258; ...
        0.227464,  -0.333333,  -0.937864; ...
        0.271054,   0.925138,  -0.33629; ...
        0.317142,  -0.830237,  -0.502561; ...
        0.352168,   0.666667,   0.688458; ...
        0.381901,  -0.925138,   0.201774; ...
        0.403548,   0.,         0.937864; ...
        0.489955,   0.830237,   0.33629; ...
        0.505832,  -0.375774,   0.803348; ...
        0.521710,   0.418215,  -0.771593; ...
        0.579632,   0.333333,   0.771593; ...
        0.583675,   0.692896,  -0.470806; ...
        0.602051,  -0.290893,  -0.771593; ...
        0.657475,  -0.597996,  -0.502561; ...
        0.722234,  -0.692896,   0.201774; ...
        0.749174,   0.0848813, -0.688458; ...
        0.784200,  -0.418215,   0.502561; ...
        0.802576,   0.597996,   0.201774; ...
        0.849436,   0.529326,  -0.201774; ...
        0.858000,   0.290893,   0.470806; ...
        0.860021,  -0.529326,  -0.150394; ...
        0.951720,   0.153552,  -0.33629; ...
        0.960284,  -0.0848813,  0.33629; ...
        1.007140,  -0.153552,  -0.067258 ]';
%
%  Set the face orders.
%
  face_order(1:face_num) = [ ...
    6, 6, 5, 6, 5, 6, 5, 6, 6, 6, ...
    5, 6, 5, 6, 5, 6, 6, 6, 5, 6, ...
    5, 5, 6, 6, 6, 5, 6, 5, 6, 6, ...
    5, 6 ];
%
%  Set faces.
%
  face_point(1:face_order_max,1:face_num) = [ ...
       30, 43, 47, 41, 29, 23; ...
       30, 23, 12,  9, 15, 27; ...
       30, 27, 35, 45, 43, -1; ...
       43, 45, 53, 59, 56, 47; ...
       23, 29, 21, 11, 12, -1; ...
       27, 15, 13, 22, 33, 35; ...
       47, 56, 54, 44, 41, -1; ...
       45, 35, 33, 42, 51, 53; ...
       12, 11,  4,  1,  3,  9; ...
       29, 41, 44, 37, 25, 21; ...
       15,  9,  3,  6, 13, -1; ...
       56, 59, 60, 58, 55, 54; ...
       53, 51, 57, 60, 59, -1; ...
       11, 21, 25, 19, 10,  4; ...
       33, 22, 24, 36, 42, -1; ...
       13,  6,  7, 17, 24, 22; ...
       54, 55, 48, 39, 37, 44; ...
       51, 42, 36, 40, 50, 57; ...
        4, 10,  8,  2,  1, -1; ...
        3,  1,  2,  5,  7,  6; ...
       25, 37, 39, 28, 19, -1; ...
       55, 58, 52, 46, 48, -1; ...
       60, 57, 50, 49, 52, 58; ...
       10, 19, 28, 26, 16,  8; ...
       36, 24, 17, 20, 32, 40; ...
        7,  5, 14, 20, 17, -1; ...
       48, 46, 34, 26, 28, 39; ...
       50, 40, 32, 38, 49, -1; ...
        8, 16, 18, 14,  5,  2; ...
       46, 52, 49, 38, 31, 34; ...
       16, 26, 34, 31, 18, -1; ...
       32, 20, 14, 18, 31, 38 ]';

  return
end


text:
%TIME_INTOVERSAMP
%
%   The purpose of this test is to determine whether special code for
%   handling the integer oversampling case in the factorization routines
%   pays off in any way.

Lr=[1048576*2,2400];
ar=[      512,   2];
Mr=[     1024, 800];

for ii=1:length(Lr)

  L=Lr(ii);

  M=Mr(ii);
  a=ar(ii);

  [L, a, M]

  N=L/a;
  c=gcd(a,M);
  p=a/c;
  q=M/c;
  d=N/q;

  f=rand(L,1);
  gf=rand(p*q,c*d);
  c1=mex_dgt_fac_1(f,gf,a,M);

  f=rand(L,1);
  gf=rand(p*q,c*d);
  c2=mex_dgt_fac_2(f,gf,a,M);

end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function path = grFindMaximalLengthPath(nodes, edges, edgeWeights)
%GRFINDMAXIMALLENGTHPATH Find a path that maximizes sum of edge weights
%
%   PATH = grFindMaximalLengthPath(NODES, EDGES, EDGE_WEIGHTS);
%   Finds a greatest geodesic path in the graph. A path between two nodes
%   is a succession of adjacent edges that link the first and last nodes.
%   the length of the path is the sum of weights of edges that constitute
%   the path.
%   A geodesic path is a path that minimizes the length of the path among
%   the set of paths between the nodes.
%   A maximal length path maximizes the length of the geodesic path between
%   couples of nodes in the graph
%
%   The result PATH is the list of edge indices that constitutes the path.
%
%   PATH = grFindMaximalLengthPath(NODES, EDGES);
%   Assumes each edge has a weight equal to 1.
%
%   See Also
%   grFindGeodesicPath
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-05-22,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% ensure weights are defined
if ~exist('edgeWeights', 'var')
    edgeWeights = ones(size(edges, 1), 1);
end

% find an extremity vertex
inds = graphPeripheralVertices(nodes, edges, edgeWeights);
ind0 = inds(end);

% find a vertex opposite to the first extremity
dists = grPropagateDistance(nodes, edges, ind0, edgeWeights);
ind1 = find(dists == max(dists), 1, 'first');

% iterate on neighbors of current node: choose next neighbor with smallest
% cumulated weight, until we are back on source node
path = [];
while true
    % find neighbor with lowest cumulated distance
    neighs = grAdjacentNodes(edges, ind1);
    neighDists = dists(neighs);
    indN = find(neighDists == min(neighDists), 1);
    ind2 = neighs(indN);

    % add edge index to the path
    indE = find(sum(ismember(edges, [ind1 ind2]), 2) == 2, 1);
    path = [path indE]; %#ok<AGROW>

    % test if path is finished or not
    if ind2 == ind0
        break;
    end
    ind1 = ind2;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function ok=checkScans(vw,scanList)
%
% Check that all scans in scanList have the same slices, numFrames, cropSizes
for iscan = 2:length(scanList)
    if find(sliceList(vw,scanList(1)) ~= sliceList(vw,scanList(iscan)))
        mrErrorDlg('Can not average these scans; they have different slices.');
    end
    if (numFrames(vw,scanList(1)) ~= numFrames(vw,scanList(iscan)))
        mrErrorDlg('Can not average these scans; they have different numFrames.');
    end
    if find(viewGet(vw,'slice dims', scanList(1)) ~= viewGet(vw,'slice dims', scanList(iscan)))
        mrErrorDlg('Can not average these scans; they have different cropSizes.');
    end
end

ok = true;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [outputScore] = getSophisticatedEdgeScore(edgeImg, labels, labelIndices, spSet)

[~,labelCount] = size(spSet);
k = 0.4;
totalSize = 0;
edgeCount = 0;
for i=1:labelCount
    tmpLabel = spSet(1,i);
    locations = labelIndices{1,tmpLabel};
    [pixelSize,~] = size(locations);
    totalSize = totalSize + pixelSize;

    for a=1:pixelSize
        if edgeImg(locations(a,1),locations(a,2)) > 0
            edgeCount = edgeCount + 1;
        end
    end
end

outputScore = edgeCount / ((totalSize).^k);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm

path: github-repos/MATLAB/zhangqianqianQQ-MachineVisionAlgorithm/MachineVisionAlgorithm-683338f6c3b1aab9fa2b80026915fe936aebf0ee/识别算法/Object Recognition based on super pixel/getSophisticatedEdgeScore.m


text:
% DEMO_fitSpline.m
%
% Fit a spline to data, which can then be evaluated by ppval

tSpan = [0,10];
nData = 100;
nKnot = 10;

tData = linspace(tSpan(1),tSpan(2),nData);
xData = sin(tData);
yData = cos(tData);

tKnot = linspace(tSpan(1),tSpan(2),nKnot);

pp = fitSpline(tData, [xData;yData], tKnot);

zSpline = ppval(pp,tData);

figure(15); clf;

subplot(2,1,1); hold on;
plot(tData,xData,'ko')
plot(tData,zSpline(1,:),'b-')
legend('data','spline');

subplot(2,1,2); hold on;
plot(tData,yData,'ko')
plot(tData,zSpline(2,:),'b-')
legend('data','spline');

repo: dscTutorials


text:
% extract the boundary vertices of all ROIs
% NOTE: requires Matlab R2014b

%%
% load_data;

%%
tic

figure;
xyz_ls=plot_ROI_3D(CONST,1:nCells);
close(gcf);

s=0.7; % shrink parameter, not affecting much
P=xyz_ls';
P(:,2)=-P(:,2); % special treatment for y
k_ls= boundary(P,s);

display('time to generate boundary:')
toc

%%
save(['data/boundary_k_s',num2str(s),'.mat'],'k_ls','s');
save('data/xyz_ls.mat','xyz_ls');

author: xiuyechen

repo: FishExplorer


text:
function [Q fcnEvals iter] = adaptiveLobatto(fcn, a, b, varargin)
% adaptiveLobatto - Numerically evaluate integral, adaptive Lobatto quadrature.
%
% function [Q fcnEvals iter] = adaptiveLobatto(fcn, a, b, varargin)
%
% (c) Matthias Conrad and Nils Papenberg (2007-08-03)
%
% Authors:
%   Matthias Conrad (e-mail: conrad@tiaco.de)
%   Nils Papenberg  (e-mail: papenber@math.uni-luebeck.de)
%
% Version:
%		Release date: 2008-08-12   Version: 1.2
%   MATLAB Version 7.5.0.338 (R2007b)
%
% Description:
%   The adaptive Lobatto algorithm programmed in an iterative not recursive
%   manner
%
% Input arguments:
%   fcn             - function to be integrated
%   a               - first point of interval
%   b               - final point of interval
%   #varargin       - further options of algorithm
%     tol           - tolerance accuracy of quadrature [ 1e-6 ]
%     parts         - initial number of partitions [ 2 ]
%     maxFcnEvals   - maximal number of function evaluations allowed [ 20000 ]
%     maxParts      - maximal number of partitions allowed [ 8000 ]
%
% Output arguments:
%   Q               - numerical integral of function fcn on [a,b]
%   fcnEvals        - number of function evaluations
%   iter            - number of iterations
%
% Details:
%   This function behavior is similar to of Matlab integrated function "quadv".
%
%   Example:
%     Q = adaptiveLobatto(@(x) [-cos(50*x); sin(x)], 0, pi, 'tol', 1e-6)
%
% References:
%   [1] Gander, W. & Gautschi, W. Adaptive Quadrature - Revisited
%       Eidgenoessische technische Hochschule Zuerich, 2000.

% check scalar limits of interval
if ~isscalar(a) || ~isscalar(b)
  error('Matlab:adaptiveLobatto:Limits',...
    'The limits of integration must be scalars.');
end

% default values
tol = 1e-6; parts = 2; maxFcnEvals = 20000; maxParts = 8000;

% rewrite default options if needed
for j = 1 : length(varargin) / 2
  eval([varargin{2 * j - 1},'=varargin{',int2str(2 * j),'};']);
end

% initial values, termination constant, parts of interval and integral value
m = parts; parts = 4 * parts + 1; Q = 0;
minH = eps(b - a) / 1024; maxResolution = 0; iter = 0;
poleWarning = 0;

% width constants
alpha = sqrt(2/3); beta = sqrt(1/5);

% check if interval has infinite boundaries, in case substitute function
if ~isfinite(a) || ~isfinite(b)
  warning('Matlab:adaptiveLobatto:infiniteInterval',...
    'The integral has an infinite interval; proceed with a substitution of function on finite interval.')
  if ~isfinite(a) && isfinite(b)
    [Q fcnEvals iter] = adaptiveLobatto(fcn, 0, b, varargin);
    fcn = @(t) infiniteLeft(t, fcn);
    a = 0; b = 1;
  elseif isfinite(a) && ~isfinite(b)
    [Q fcnEvals iter] = adaptiveLobatto(fcn, a, 0, varargin);
    fcn = @(t) infiniteRight(t, fcn);
    a = 0; b = 1;
  else
    fcn = @(t) infiniteBoth(t, fcn);
    a = - pi / 2; b = pi / 2;
  end
end

% initialize grid
t = linspace(a, b, m + 1);
A = t(1:end-1); B = t(2:end);

% widths and midpoints of intervals
H = diff(t)/2; J = (A + B) / 2;

% grid points
F = -alpha * H + J; D = -beta * H + J; C =  J;
E =  beta * H + J; G =  alpha * H + J;
t = [A; F; D; C; E; G; B]; t = t(:);

% function evaluations
y = fcn([A, F, D, C, E, G, B]); fcnEvals = 7 * m;

% avoid infinities at start point of interval
if any(~isfinite(y(:,1)))
  y(:,1) = fcn(a + eps(superiorfloat(a,b)) * (b - a));
  fcnEvals = fcnEvals + 1;
  poleWarning = 1;
end

% avoid infinities at end point of interval
if any(~isfinite(y(:, end)))
  y(:, end) = fcn(b - eps(superiorfloat(a,b)) * (b - a));
  fcnEvals = fcnEvals + 1;
  poleWarning = 1;
end

% poles at initial points
if ~isempty(find(~isfinite(max(abs(y))))), poleWarning = 1; end

% hand over function values
yA = y(:,     1 :   m); yF = y(:,   m+1 : 2*m); yD = y(:, 2*m+1 : 3*m);
yC = y(:, 3*m+1 : 4*m); yE = y(:, 4*m+1 : 5*m); yG = y(:, 5*m+1 : 6*m);
yB = y(:, 6*m+1 : end);

% dimension of parallel integration
n = size(yA,1);

% adaptive Lobatto iteration
while 1

  % number of iteration
  iter = iter + 1;

  % four point Lobatto formula
  Q1 = kron(H, ones(n,1)) / 6 .* (yA + 5 * (yD + yE) + yB);
  % seven point Kronrod formula
  Q2 = kron(H, ones(n,1)) / 1470 .* (77 * (yA + yB) + 432 * (yF + yG) + 625 * (yD + yE) + 672 * yC);

  % difference of Lobatto formulas
  diffQ = Q2 - Q1; diffQ(find(isnan(diffQ))) = 0;

  % intervals which do not fulfill termination criterion
  idx = find(max(abs(diffQ), [], 1) > tol);

  % intervals fulfill termination criterion
  idxQ = setdiff(1:length(A), idx);

  % check stop criterions
  STOP1 = isempty(idx); % check regular termination
  STOP2 = fcnEvals > maxFcnEvals;  % check maximal function evaluations
  STOP3 = 5 * length(idx) > maxParts; % check maximal partition

  % regular termination
  if STOP1
    Q = Q + sum(Q2, 2);
    break
  end

  % check if maximal resolution reached
  idxH = find(abs(H) < minH);
  if ~isempty(idxH)
    Q = Q + sum(Q2(idxH), 2);
    idx = setdiff(idx, idxH);
    idxQ = setdiff(idxQ, idxH);
    maxResolution = 1;
    % termination criterion
    if isempty(idx), break, end
  end

  % maximal function evaluations reached
  if STOP2
    warning('Matlab:adaptiveLobatto:MaxEvaluations',...
      'The maximal number of function evaluations reached; singularity likely.')
    Q = Q + sum(Q2, 2);
    break
  end

  % maximal partition reached
  if STOP3
    warning('Matlab:adaptiveLobatto:parts',...
      'The maximal number of parts reached.')
    Q = Q + sum(Q2, 2);
    break
  end

  % update quadrature value
  Q = Q + sum(Q2(:,idxQ) + diffQ(:,idxQ) / 15, 2);

  % number of intervals
  m = 6 * length(idx);

  % initialize t
  t = zeros(1, 6 * length(idx));

  % hand over new start points A
  t(1:6:end) = A(idx); t(2:6:end) = F(idx); t(3:6:end) = D(idx);
  t(4:6:end) = C(idx); t(5:6:end) = E(idx); t(6:6:end) = G(idx);
  A = t;

  % hand over new end points B
  t(1:6:end) = F(idx); t(2:6:end) = D(idx); t(3:6:end) = C(idx);
  t(4:6:end) = E(idx); t(5:6:end) = G(idx); t(6:6:end) = B(idx);
  B = t;

  y = zeros(n, 6 * length(idx));
  % hand over new start values A
  y(:,1:6:end) = yA(:,idx); y(:,2:6:end) = yF(:,idx); y(:,3:6:end) = yD(:,idx);
  y(:,4:6:end) = yC(:,idx); y(:,5:6:end) = yE(:,idx); y(:,6:6:end) = yG(:,idx);
  yA = y;

  % hand over new end values B
  y(:,1:6:end) = yF(:,idx); y(:,2:6:end) = yD(:,idx); y(:,3:6:end) = yC(:,idx);
  y(:,4:6:end) = yE(:,idx); y(:,5:6:end) = yG(:,idx); y(:,6:6:end) = yB(:,idx);
  yB = y;

  % widths and midpoints of intervals
  H = (B - A) / 2; J = (A + B) / 2;

  % calculate new mid points
  F = -alpha * H + J; D = -beta * H + J; C =  J;
  E =  beta * H + J; G =  alpha * H + J;

  % function evaluations
  y = fcn([F, D, C, E, G]); fcnEvals = fcnEvals + 5 * m;

  % poles at new points
  if ~isempty(find(~isfinite(max(abs(y))))), poleWarning = 1; end

  % hand over new midpoint values of F D C E and G
  yF = y(:,     1 :   m); yD = y(:,   m+1 : 2*m); yC = y(:, 2*m+1 : 3*m);
  yE = y(:, 3*m+1 : 4*m); yG = y(:, 4*m+1 : 5*m);

end

% display warnings
if any(~isfinite(Q))
  warning('Matlab:adaptiveLobatto:Infinite',...
    'The Quadrature of the function reached infinity or is Not-a-Number.')
end
if maxResolution
  warning('Matlab:adaptiveLobatto:MaxResolution',...
    'The maximal resolution of partial interval reached; singularity likely.')
end
if poleWarning
  warning('Matlab:adaptiveLobatto:PoleDetection',...
    'A detection of a pole; singularity likely.')
end

return

% substitute function interval [-inf, 0] on [0, 1]
function f = infiniteLeft(t, fcn)
f = fcn(log(t));
f = f ./ kron(ones(size(f,1),1), t);
return

% substitute function interval [0, inf] on [0, 1]
function f = infiniteRight(t, fcn)
f = fcn(-log(t));
f = f ./ kron(ones(size(f,1),1), t);
return

% substitute function interval [-inf, inf] on [-pi / 2, pi / 2]
function f = infiniteBoth(t, fcn)
f = fcn(tan(t));
f = f ./ kron(ones(size(f,1),1), cos(t).^2);
return


text:
function out = TSTL_predict(y, plen, NNR, stepSize, pmode, embedParams)
% TSTL_predict  Local constant iterative time-series prediction.
%
% References TSTOOL code 'predict', which does local constant iterative
% prediction for scalar data using fast nearest neighbour searching.
% There are four modes available for the prediction output.
%
% TSTOOL: http://www.physik3.gwdg.de/tstool/
%
%---INPUTS:
%
% y, scalar column vector time series
%
% plen, prediction length in samples or as proportion of time series length NNR,
%
% NNR, number of nearest neighbours
%
% stepSize, number of samples to step for each prediction
%
% pmode, prediction mode, four options:
%           (i) 0: output vectors are means of images of nearest neighbours
%           (ii) 1: output vectors are distance-weighted means of images
%                     nearest neighbours
%           (iii) 2: output vectors are calculated using local flow and the
%                    mean of the images of the neighbours
%           (iv) 3: output vectors are calculated using local flow and the
%                    weighted mean of the images of the neighbours
% embedParams, as usual to feed into BF_Embed, except that now you can set
%              to zero to not embed.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Check inputs, set defaults
% ------------------------------------------------------------------------------
doPlot = false; % plot outputs to figure (e.g., for debugging)

% (*) Prediction length, plen (the length of the output time series)
if nargin < 2 || isempty(plen)
    plen = 1; % output the same length (proportion)
end
% (proportion set after embedding, as the embedding will lose points
% according to the dimension of the space)

% (*) number of neighest neighbours, NNR
if nargin < 3 || isempty(NNR)
    NNR = 1; % use 1 nearest neighbour
end

% (*) stepSize (in samples)
if nargin < 4 || isempty(stepSize)
    stepSize = 2;
end

% (*) prediction mode, pmode:
if nargin < 5 || isempty(pmode)
    pmode = 0; % output vectors are means of the images of nearest neighbours
end

% (*) embedParams
if nargin < 6 || isempty(embedParams)
    embedParams = {'ac','fnnmar'};
    fprintf(1,'Using default embedding using autocorrelation for tau and Cao''s method for m\n');
end

% ------------------------------------------------------------------------------
%% Embed the scalar time series by time-delay method
% ------------------------------------------------------------------------------
% embedpn = BF_Embed(y,embedParams{1},embedParams{2},2);
% delay = embedpn(1);
% dim = embedpn(2);
if iscell(embedParams)
    s = BF_Embed(y,embedParams{1},embedParams{2},1); % last in
elseif embedParams == 0
    s = signal(y);
end

if ~isa(s,'signal')
    out = NaN; return
end

Ns = length(data(s));
if Ns < 50
    error('This is a very short time series! :(');
end
y = y(1:Ns); % for statistical purposes...
if plen > 0 && plen <= 1
    plen = floor(plen*Ns); % specify a proportion of the time series length
end

% ------------------------------------------------------------------------------
%% Run the code
% ------------------------------------------------------------------------------
try
    rs = predict(s, plen, NNR, stepSize, pmode);
catch
    error('TSTOOL''s predict function didn''t run correctly')
end

y_pred = data(rs);
y_pred1 = y_pred(:,1); % for this embedding dimension (?)

if doPlot
    figure('color','w'); box('on'); view(rs);

    figure('color','w'); box('on');
    hold off; plot(y,'k'), hold on; plot(y_pred1,'m'), hold off;
end

% ------------------------------------------------------------------------------
%% Compare the output to the properties of the true time series
% ------------------------------------------------------------------------------

% actual basic statistical properties
out.pred1mean = mean(y_pred1);
out.pred1std = std(y_pred1);
out.pred1maxc = abs(max(y_pred1)-max(y));
out.pred1maxc = abs(max(y_pred(:))-max(y));
out.pred1minc = abs(min(y_pred1)-min(y));
out.predminc = abs(min(y_pred(:))-min(y));
out.pred1rangec = abs(range(y_pred1)/range(y)-1);

% look at structure in cross correlation function, xcf
% (requires that prediction length the same as the time series itself)
[xcf, lags] = xcorr(y,y_pred1,'coeff');
% plot(lags,xcf);

out.maxabsxcf = max(abs(xcf)); % maximum of cross-correlation function; where it occurs
out.maxabsxcflag = lags(find(abs(xcf) == out.maxabsxcf,1,'first'));
out.maxxcf = max(xcf); % maximum positive cross-correlation
out.maxxcflag = lags(find(xcf == out.maxxcf,1,'first'));
out.meanxcf = mean(xcf);
out.minxcf = min(xcf);
out.stdxcf = std(xcf);

out.pred1_ac1 = CO_AutoCorr(y_pred1,1,'Fourier'); % autocorrelation at lag one of prediction
out.pred1ac1diff = abs(out.pred1_ac1 - CO_AutoCorr(y,1,'Fourier')); % difference in autocorrelations of prediction and original
out.pred1_ac2 = CO_AutoCorr(y_pred1,2,'Fourier'); % autocorrelation at lag one of prediction
out.pred1ac2diff = abs(out.pred1_ac2 - CO_AutoCorr(y,2,'Fourier')); % difference in autocorrelations of prediction and original
out.pred_tau_comp = CO_FirstCrossing(y_pred1,'ac',0,'continuous')/CO_FirstCrossing(y,'ac',0,'continuous'); % ratio of first zero crossings of autocorrelation function

% Autocorrelation structure:
acs_y = CO_AutoCorr(y,1:10,'Fourier');
acs_y_pred1 = CO_AutoCorr(y_pred1,1:10,'Fourier');
out.acs1_10_sumabsdiffpred1 = sum(abs(acs_y - acs_y_pred1));

% mean square residuals: this will likely be a bad measure (as may be out
% of sync)
out.pred1rmsres = sqrt(mean((y-y_pred1).^2));

% align at best positive cross-correlation and then look at residuals
if out.maxxcflag > 0
    y_lagged = y(out.maxxcflag:end);
    Nlag = length(y_lagged);
    y_pred1_lagged = y_pred1(1:Nlag);
elseif out.maxxcflag == 0
    y_lagged = y;
    y_pred1_lagged = y;
    Nlag = length(y);
else % negative
    y_pred1_lagged = y_pred1(-out.maxxcflag:end);
    Nlag = length(y_pred1_lagged);
    y_lagged = y(1:Nlag);
end

% hold off; plot(y_pred1_lagged);
% hold on; plot(y_lagged,'r');
out.Nlagxcorr = Nlag;
res = y_lagged - y_pred1_lagged;
out.bestpred1rmsres = sqrt(mean(res.^2)); % rms residuals
out.ac1bestres = CO_AutoCorr(res,1,'Fourier'); % autocorrelation of residuals

% now look at fraction of points that are within a threshold of each
% other...
fracresfn = @(x) sum(abs(res) < x)/Nlag;
out.fracres005 = fracresfn(0.05); %sum(abs(res)<0.05)/Nlag;
out.fracres01 = fracresfn(0.1); %sum(abs(res)<0.1)/Nlag;
out.fracres02 = fracresfn(0.2); %sum(abs(res)<0.2)/Nlag;
out.fracres03 = fracresfn(0.3); %sum(abs(res)<0.3)/Nlag;
out.fracres05 = fracresfn(0.5); %sum(abs(res)<0.5)/Nlag;

% now look at fraction of points within a circle of (time-measurement) radius
% near a real point in the time series
% this could be done using the closest neighbour of the simulation to the
% real time series

% Could compare different dimensions rather than just the first... find the
% best... etc.

% There are heaps of metrics you could use to compare... Ultimately may be
% able to implement model outputs as rows so as to compare across many
% measures...

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function varargout=corr(varargin)
% INTERNAL FUNCTION
%

[varargout{1:nargout}]=corr(varargin{:});
end

author: jmaih


text:
function msg = apm_arx(p,m,ny,nu,alpha,beta,gamma,name)

if nargin <= 6,
    error('apm_arx: Input ARX Model')
end
if nargin == 7,
    name = 'arx';
end
if nargin == 8,
    % trim whitespace
    filename = strtrim(name);
    % parse file name and extension
    [path,name,ext] = fileparts(name);
end
filename = [name '.apm'];

% check sizes for alpha
[alpha_ny,alpha_p] = size(alpha);
if (alpha_ny~=ny),
    disp(['input ny: ', num2str(ny)])
    disp(['alpha ny: ', num2str(alpha_ny)])
    error('ARX Size mismatch for alpha')
end
if (alpha_p~=p),
    disp(['input p: ', num2str(p)])
    disp(['alpha p: ', num2str(alpha_p)])
    error('ARX Size mismatch for alpha')
end

% check sizes for beta
[beta_nu,beta_m,beta_p] = size(beta);
if (beta_m~=m),
    disp(['input m: ', num2str(m)])
    disp(['beta m: ', num2str(beta_m)])
    error('ARX Size mismatch for beta')
end
if (beta_nu~=nu),
    disp(['input nu: ', num2str(nu)])
    disp(['beta nu: ', num2str(beta_nu)])
    error('ARX Size mismatch for beta')
end
if (beta_p~=p),
    disp(['input p: ', num2str(p)])
    disp(['beta p: ', num2str(beta_p)])
    error('ARX Size mismatch for alpha')
end

fid = fopen(filename,'w');
fprintf( fid,'\n');
fprintf( fid,'Objects \n');
fprintf( fid,['  ' name ' = arx\n']);
fprintf( fid,'End Objects \n');
fprintf( fid,'\n');
fprintf( fid,'Connections\n');
fprintf( fid,['  u[1:%d] = ' name '.u[1:%d]\n'],m,m);
fprintf( fid,['  y[1:%d] = ' name '.y[1:%d]\n'],p,p);
fprintf( fid,'End Connections\n');
fprintf( fid,'\n');
fprintf( fid,'Model \n');
fprintf( fid,'  Parameters \n');
if (m==1)
    fprintf( fid,'    u = 0\n');
else
    fprintf( fid,'    u[1:%d] = 0\n',m);
end
fprintf( fid,'  End Parameters \n');
fprintf( fid,'\n');
fprintf( fid,'  Variables \n');
if (p==1)
    fprintf( fid,'    y = 0\n');
else
    fprintf( fid,'    y[1:%d] = 0\n',p);
end
fprintf( fid,'  End Variables \n');
fprintf( fid,'\n');
fprintf( fid,'  Equations \n');
fprintf( fid,'    ! add any additional equations here \n');
fprintf( fid,'  End Equations \n');
fprintf( fid,'End Model \n');
fprintf( fid,'\n');
fprintf( fid,['File ' name '.txt\n']);
fprintf( fid,'  %d      ! m=number of inputs\n',m);
fprintf( fid,'  %d      ! p=number of outputs\n',p);
fprintf( fid,'  %d      ! nu=number of input terms\n',nu);
fprintf( fid,'  %d      ! ny=number of output terms\n',ny);
fprintf( fid,'End File\n');
fprintf( fid,'\n');
fprintf( fid,'! Alpha matrix (ny x p)\n');
fprintf( fid,['File ' name '.alpha.txt \n']);
for i = 1:ny
    for j = 1:p
        fprintf( fid,'  ');
        if j<=p-1
            fprintf( fid,'%d , ', alpha(i,j));
        else
            fprintf( fid,'%d ', alpha(i,j)); % no comma
        end
    end
    fprintf( fid,'\n');
end
fprintf( fid,'End File \n');
fprintf( fid,'\n');
fprintf( fid,'! Beta matrix (p x (nu x m))\n');
fprintf( fid,['File ' name '.beta.txt \n']);
for i = 1:p
    for j = 1:nu
        fprintf( fid,'  ');
        for k = 1:m
            if k<=m-1
                fprintf( fid,'%d , ', beta(j,k,i));
            else
                fprintf( fid,'%d ', beta(j,k,i));
            end
        end
        fprintf( fid,'\n');
    end
end
fprintf( fid,'End File \n');
fprintf( fid,'\n');
fprintf( fid,'! Gamma vector (p x 1)\n');
fprintf( fid,['File ' name '.gamma.txt \n']);
for i = 1:p
    fprintf( fid,'%d ', gamma(i));
    fprintf( fid,'\n');
end
fprintf( fid,'End File \n');
fclose(fid);

msg = ['Created ARX (Auto-Regressive eXogenous inputs) Model: ' filename];

end

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function fem=springs(fem,iter);
% SPRING is legacy code that moves internal points of a finite
%  element mesh such that the mesh quality is improved.  The movement
%  of nodes is accomplished with a Laplacian smoothing action, so
%  that each point not located along the boundary is moved toward
%  the center of mass of the polygon formed by the adjacent
%  triangles. The process is repeated according to the setting
%  of the ITER variable.
%
%       Calls laplace_smooth.m with specified input.
% FileName:  springs.m
% Written By:  unknown
% Date Last Modified:
%    May 9, 2007 -- Chris Massey, NRL Code 7322, Stennis Spc Cnt., MS
%                   Deleted code and replaced with a call to
%                   laplace_smooth.m.
%    May 15, 2007 -- Chris Massey, -- Fixed bug, should send over the
%                   node list not the element list.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin==1,
    iter=5;
end

fem = laplace_smooth(fem,1:length(fem.x),iter);

return

author: CHLNDDEV


text:
classdef YawTermCondition < AbstractEventTerminationCondition
    %YawTermCondition Summary of this class goes here
    %   Detailed explanation goes here

    properties
        steeringModel(1,1) AbstractSteeringModel = RollPitchYawPolySteeringModel.getDefaultSteeringModel();
        targetYawAngle(1,1) double = 0;
        bodyInfo KSPTOT_BodyInfo
    end

    methods
        function obj = YawTermCondition(targetYawAngle)
            obj.targetYawAngle = targetYawAngle;
        end

        function evtTermCondFcnHndl = getEventTermCondFuncHandle(obj)
            evtTermCondFcnHndl = @(t,y) obj.eventTermCond(t,y, obj.targetYawAngle, obj.steeringModel, obj.bodyInfo);
        end

        function initTermCondition(obj, initialStateLogEntry)
            obj.steeringModel = initialStateLogEntry.steeringModel;
            obj.bodyInfo = initialStateLogEntry.centralBody;
        end

        function name = getName(obj)
            name = sprintf('Yaw Angle (%.3f deg)', rad2deg(obj.targetYawAngle));
        end

        function tf = shouldBeReinitOnRestart(obj)
            tf = false;
        end

        function params = getTermCondUiStruct(obj)
            params = struct();

            params.paramName = 'Target Yaw Angle';
            params.paramUnit = 'deg';
            params.useParam = 'on';
            params.useStages = 'off';
            params.useTanks = 'off';
            params.useEngines = 'off';
            params.useStopwatches = 'off';

            params.value = rad2deg(obj.targetYawAngle);
            params.refStage = LaunchVehicleStage.empty(1,0);
            params.refTank = LaunchVehicleEngine.empty(1,0);
            params.refEngine = LaunchVehicleEngine.empty(1,0);
            params.refStopwatch = LaunchVehicleStopwatch.empty(1,0);
        end

        function optVar = getNewOptVar(obj)
            optVar = YawAngleTermCondOptimVar(obj);
        end

        function optVar = getExistingOptVar(obj)
            optVar = obj.optVar;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end
    end

    methods(Static)
        function termCond = getTermCondForParams(paramValue, stage, tank, engine, stopwatch)
            termCond = YawTermCondition((paramValue));
        end
    end

    methods(Static, Access=private)
        function [value,isterminal,direction] = eventTermCond(t,y, targetYawAngle, steeringModel, bodyInfo)
            ut = t;
            rVect = y(1:3);
            vVect = y(4:6);

            dcm = steeringModel.getBody2InertialDcmAtTime(ut, rVect, vVect, bodyInfo);
            [~,~,yawAngle] = computeEulerAnglesFromInertialBodyAxes(ut, rVect, vVect, bodyInfo, dcm(:,1), dcm(:,2), dcm(:,3));

            yawAngle = AngleZero2Pi(yawAngle);

            value = yawAngle - targetYawAngle;
            isterminal = 1;
            direction = 0;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function x = Adj_Evaluate_FT(y,shift,boxlen,center,w)
% Adj_Evaluate_FT -- 1d unequispaced adjoint Fourier transform
% (with windowing). Adjoint of Evaluate_FT.
% Usage:
%   x = Evaluate_FT(y,shift,boxlen,center,w)
% Inputs:
%   y	      vector of length 2n
%   shift     vector of shifts
%   boxlen    half length of interval around each value of shift
%   center    boolean variable: 0/1 = unbiased/biased FT
%   w         tapering window of size 2*boxlen
% Outputs:
%   x      vector of length: n
% Description:
%  For w = 1, evaluates the Adj FT of a signal irregularly
%  sampled in frequency omega(j,k) = shift(j) + k, -boxlen <= k
%  <boxlen If center = 0,
%
%  x(t) = sum_m exp(i 2pi omega(m) t/n) y(k), -n/2 <= t < n/2
%
%  If center = 1,
%
%  x(t) = sum_m exp(i 2pi omega(m) t/n) y(k), 0 <= t < n
%
%  For arbitrary windows, y is first tapered with w.
%
%  Adj_Evaluate_FT uses two different methods for computing the
%  unesquispaced Fourier transform. If the number of shifts is less
%  than 8, it uses an exact method, otherwise it uses the
%  accelerated Adjoint USFFT.
%  See Also
%    Evaluate_FT, Adj_USFFT, Adj_USFT_simple
%
% By Emmanuel candes, 2003-2004

if nargin < 5,
  w = ones(1,2*boxlen);
end

if nargin < 4,
  center = 0;
end

n  = length(y)/2;
boxcnt = length(shift);

if (boxcnt <= 16)
  x = Adj_USFT_simple(y,shift,boxlen,center,w);
else
  % Make frequency grid
  k = (-boxlen):(boxlen -1);
  omega = meshgrid(shift,k) + meshgrid(k,shift).';
  omega = 2*pi*(omega(:))./n;

  y = repmat(w(:),boxcnt,1) .* y; % Tapering
  x = Adj_USFFT(n,y,omega,16,6,center)./sqrt(n); % Take Adjoint USFFT
end

author: thomaskuestner


text:
classdef obspysource < dataretrieval.spatiotemporal_file
   %UNTITLED2 Summary of this class goes here
   %   Detailed explanation goes here

   properties
   end

   methods
      function data = retrieve(obj, where_, from_ , until_)
         error('unimplemented obspy data source');
      end

      function out_ = translate(obj, in_)

      end
   end

end

author: geoscience-community-codes

repo: GISMO