text:
function [model] = BCPF_MP(Y, varargin)
% Bayesian CP Factorization using Gaussian Mixture Priors for Image Completion
% Author : Qibin Zhao  2014
%
% -----------------------------------------------------------------------
%  [model] = BCPF_MP(Y, 'PARAM1', val1, 'PARAM2', val2, ...)
%
%  INPUTS
%     Y              - Input tensor
%     'obs'          - Binary (0-1) tensor indicating missing entries
%                      (0: missing; 1: observed)
%     'init'         - Initialization method
%                     - 'ml'  : SVD initilization (default)
%                     - 'rand': Random matrices
%     'maxRank'      - The initialization of rank (larger than true rank)
%     'dimRed'       - 1: Remove unnecessary components automaticly (default)
%                    - 0: Not remove
%     'maxiters'     - max number of iterations (default: 100)
%     'tol'          - lower band change tolerance for convergence dection
%                      (default: 1e-5)
%     'noise'        - whether noise is updated
%                        - 'on': update noise parameter (default)
%                        - 'off': fixed noise parameter (1e-5)
%     'predVar'      - Predictive distribution
%                         - 1:  compute and output
%                         - 0:  doesnot compute  (default)
%     'verbose'      - visualization of results
%                       - 0: no
%                       - 1: text (default)
%                       - 2: online display image
%                       - 3: show factors by image
%                       - 4: show factors by hinton plot (very slow)
%   OUTPUTS
%      model         - Model parameters and hyperparameters
% -----------------------------------------------------------------------
%
%   Example:
%
%     [model] = BCPF_MP(Y, 'obs', O, 'init', 'rand', 'maxRank', 10, 'dimRed', 1, 'maxiters', 100, ...
%                                'tol', 1e-6, 'verbose', 3);
%
% < Bayesian CP Factorization of Incomplete Image using Gaussian Mixture Priors >
% Copyright (C) 2014  Qibin Zhao
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
warning off; %#ok<WNOFF>
randn('state',1); rand('state',1); %#ok<RAND>
dimY = size(Y);
N = ndims(Y);

%% Set parameters from input or by using defaults
ip = inputParser;
ip.addParamValue('obs', ones(dimY), @(x) (isnumeric(x) || islogical(x)) );
ip.addParamValue('init', 'rand', @(x) (ismember(x,{'ml','rand'})));
ip.addParamValue('maxRank', max(dimY), @isscalar);
ip.addParamValue('maxiters', 100, @isscalar);
ip.addParamValue('tol', 1e-5, @isscalar);
ip.addParamValue('verbose', 1, @isscalar);
ip.addParamValue('noise', 'on', @(x)ismember(x,{'on','off'}));
ip.addParamValue('dimRed', 1, @isscalar);
ip.addParamValue('predVar', 0, @isscalar);
ip.addParamValue('nd', 1, @isscalar);
ip.parse(varargin{:});

O     = ip.Results.obs;
init  = ip.Results.init;
R   = ip.Results.maxRank;
maxiters  = ip.Results.maxiters;
tol   = ip.Results.tol;
verbose  = ip.Results.verbose;
DIMRED   = ip.Results.dimRed;
noise = ip.Results.noise;
predVar = ip.Results.predVar;
nd = ip.Results.nd;

%% Initialization
Y = tensor(Y.*O);
O = tensor(O);
nObs = sum(O(:));

a_gamma0     = 1e-6;
b_gamma0     = 1e-6;
if  strcmp(noise,'on')
    a_beta0      = 1e-6;
    b_beta0      = 1e-6;
else
    a_beta0      = 1e-1;
    b_beta0      = 1e-6;
end
gammas = ones(R,1);
beta = 1e4;
dscale =1;

W = cell(N,1);
W1 = cell(N,1);
oR = nObs/prod(dimY);
for n=1:N
    W{n} = zeros(dimY(n),dimY(n));
    for i=1:dimY(n)
        for j=1:dimY(n)
            W{n}(i,j) = exp(-2*oR*abs(i-j)^2);
        end
    end
    W1{n} = W{n} - diag(diag(W{n}));
    W1{n} = nd*5*bsxfun(@times, W1{n}, 1./sum(W1{n},2));
end

switch init,
    case 'ml'    % Maximum likelihood
        Z = cell(N,1);
        ZSigma = cell(N,1);
        if ~isempty(find(O==0))
            %   Y(find(O==0)) = sum(Y(:))/nObs;
            Y1 = Y;
            for n=1:N
                Y1 = ttm(Y1, W{n}, n);
            end
            Y(find(O==0)) = Y1(find(O==0));
        end
        for n = 1:N
            ZSigma{n} = (repmat(eye(R), [1 1 dimY(n)]));
            [U, S, V] = svd(double(tenmat(Y,n)), 'econ');
            if R <= size(U,2)
                Z{n} = U(:,1:R)*(S(1:R,1:R)).^(0.5);
            else
                Z{n} = [U*(S.^(0.5)) randn(dimY(n), R-size(U,2))];
            end
        end
        Y = Y.*O;
    case 'rand'   % Random initialization
        Z = cell(N,1);
        ZSigma = cell(N,1);
        for n = 1:N
            Z{n} = rand(dimY(n),R);
            ZSigma{n} = (repmat(eye(R), [1 1 dimY(n)]));
            if n<3
                %   Z{n} = cholcov(W{n})'*Z{n};
                Z{n} = Z{n} +  W1{n}*Z{n};
            end
        end
end

% --------- E(aa') = cov(a,a) + E(a)E(a')----------------
EZZT = cell(N,1);
for n=1:N
    EZZT{n} = (reshape(ZSigma{n}, [R*R, dimY(n)]))';
end

Fit =0;
LB = 0;
X = double(ktensor(Z));

%% Create figures
if verbose >2,
    scrsz = get(0,'ScreenSize');
    h1 = figure('Position',[scrsz(3)*0.2 scrsz(4)*0.3 scrsz(3)*0.6 scrsz(4)*0.4]);
    figure(h1);
    switch verbose,
        case 4,
            subplot(2,3,1); hintonDiagram(Z{1}); title('Mode-1'); ylabel('Length of #-mode');
            subplot(2,3,2); hintonDiagram(Z{2}); title('Mode-2'); xlabel('Latent dimensions');
            if N>=3, subplot(2,3,3); hintonDiagram(Z{3}); title('Mode-3'); end
        case 3,
            subplot(2,3,1); imagesc(Z{1}); title('Mode-1'); ylabel('Length of #-mode');
            subplot(2,3,2); imagesc(Z{2}); title('Mode-2'); xlabel('Latent dimensions');
            if N>=3, subplot(2,3,3); imagesc(Z{3}); title('Mode-3');end
    end
    subplot(2,3,4); bar(gammas); title('Posterior mean of \lambda'); xlabel('Latent components'); ylabel(''); axis tight;
    subplot(2,3,5); plot(LB, '-r.','LineWidth',1.5,'MarkerSize',10 ); title('Lower bound'); xlabel('Iteration');  grid on;
    subplot(2,3,6); plotGamma(a_beta0, a_beta0); title('Posterior pdf'); xlabel('Noise precision \tau');grid on;
    set(findall(h1,'type','text'),'fontSize',12);
    drawnow;
end
if verbose ==2;
    h3 = figure;
    temp = 255.*(X-min(X(:)))/(max(X(:))-min(X(:)));
    imshow(uint8(temp));
    title(['Iter.= '  num2str(0),',  Rank = ' num2str(R)],'FontSize', 13, 'color','b');
    tic;
    xlabel(['(BCPF-MP)  Time: ' num2str(round(toc)), ' seconds'],'FontSize', 13, 'color','b');
    drawnow;
end

%% Model learning
for it=1:maxiters,
    %% Update factor matrices
    Aw = diag(gammas);
    for n=1:N
        % compute E(Z_{\n}^{T} Z_{\n})
        ENZZT = reshape(khatrirao_fast(EZZT{[1:n-1, n+1:N]},'r')' * double(tenmat(O,n)'), [R,R,dimY(n)]);
        % compute E(Z_{\n})
        FslashY = khatrirao_fast(Z{[1:n-1, n+1:N]},'r')' * tenmat(Y.*O, n)';
        for i=1:dimY(n)
            ZSigma{n}(:,:,i) = (beta * ENZZT(:,:,i) + Aw )^(-1);
            Z{n}(i,:) = (beta * ZSigma{n}(:,:,i) * FslashY(:,i))';
        end
        if n<3
            Z{n} = Z{n} + W1{n}*Z{n};
        end
        EZZT{n} = (reshape(ZSigma{n}, [R*R, dimY(n)]) + khatrirao_fast(Z{n}',Z{n}'))';
    end

    %% Update latent tensor X
    X = double(ktensor(Z));

    %% Update hyperparameters gamma
    a_gammaN = (0.5*sum(dimY) + a_gamma0)*ones(R,1);
    b_gammaN = 0;
    for n=1:N
        b_gammaN = b_gammaN + diag(Z{n}'*Z{n}) + diag(sum(ZSigma{n},3));
    end
    b_gammaN = b_gamma0 + 0.5.* b_gammaN;
    gammas = a_gammaN./b_gammaN;

    %% update noise beta
    %  The most time and space consuming part
    if 0 % save time but large space needed
        EX2 =  O(:)' * khatrirao_fast(EZZT,'r') * ones(R*R,1);
    else  % save space but slow
        temp1 = cell(N,1);
        EX2 =0;
        for i =1:R
            for n=1:N
                temp1{n} = EZZT{n}(:,(i-1)*R+1: i*R);
            end
            EX2 = EX2 + O(:)' * khatrirao_fast(temp1,'r')* ones(R,1);
        end
    end
    err = Y(:)'*Y(:) - 2*Y(:)'*X(:) + EX2;
    if  strcmp(noise,'on')
        a_betaN = a_beta0 + 0.5*nObs;
        b_betaN = b_beta0 + 0.5*err;
    else
        a_betaN = a_beta0;
        b_betaN = b_beta0;
    end
    beta = a_betaN/b_betaN;
    Fit = 1 - sqrt(sum(err(:)))/norm(Y(:));

    %% Lower bound
    temp1 = -0.5*nObs*safelog(2*pi) + 0.5*nObs*(psi(a_betaN)-safelog(b_betaN)) - 0.5*(a_betaN/b_betaN)*err;
    temp22 =0;
    for n=1:N
        temp22= temp22 + Z{n}'*Z{n} + sum(ZSigma{n},3);
    end
    temp2 = -0.5*R*sum(dimY)*safelog(2*pi) + 0.5*sum(dimY)*sum(psi(a_gammaN)-safelog(b_gammaN)) -0.5*trace(diag(gammas)* temp22);
    temp3 = sum(-safelog(gamma(a_gamma0)) + a_gamma0*safelog(b_gamma0) -  b_gamma0.*(a_gammaN./b_gammaN) + (a_gamma0-1).*(psi(a_gammaN)-safelog(b_gammaN)));
    temp4 = -safelog(gamma(a_beta0)) + a_beta0*safelog(b_beta0) + (a_beta0-1)*(psi(a_betaN)-safelog(b_betaN)) - b_beta0*(a_betaN/b_betaN);
    temp5=0;
    for n=1:N
        for i=1:size(ZSigma{n},3)
            temp5 = temp5 + 0.5*safelog(det(ZSigma{n}(:,:,i))) + 0.5*R*(1+safelog(2*pi));
        end
    end
    temp6 = sum(safelog(gamma(a_gammaN)) - (a_gammaN-1).*psi(a_gammaN) -safelog(b_gammaN) + a_gammaN);
    temp7 = safelog(gamma(a_betaN)) - (a_betaN-1)*psi(a_betaN) -safelog(b_betaN) + a_betaN;
    LB(it) = temp1 + temp2 + temp3 + temp4 + temp5 + temp6 + temp7;

    %% Prune irrelevant dimensions?
    Zall = cell2mat(Z);
    comPower = diag(Zall' * Zall);
    comTol = sum(dimY)*eps(norm(Zall,'fro'));
    rankest = sum(comPower> comTol );
    if max(rankest)==0
        disp('Rank becomes 0 !!!');
        break;
    end
    if DIMRED==1  && it >=2,
        if R~= max(rankest)
            indices = comPower > comTol;
            gammas = gammas(indices);
            temp = ones(R,R);
            temp(indices,indices) = 0;
            temp = temp(:);
            for n=1:N
                Z{n} = Z{n}(:,indices);
                ZSigma{n} = ZSigma{n}(indices,indices,:);
                EZZT{n} = EZZT{n}(:, temp == 0);
            end
            R = max(rankest);
        end
    end

    %% Display progress
    if it>2
        LBRelChan = abs(LB(it) - 2*LB(it-1) + LB(it-2))/-LB(2);
    else
        LBRelChan = NaN;
    end
    if verbose,
        fprintf('Iter. %d: RelChan = %g, Fit = %g, R = %d \n', it, LBRelChan, Fit, rankest);
    end

    %% visualize online results
    if verbose >2 ,
        switch verbose,
            case 4,
                set(0,'CurrentFigure',h1);
                subplot(2,3,1); hintonDiagram(Z{1}); title('Mode-1'); ylabel('Length of #-mode');
                subplot(2,3,2); hintonDiagram(Z{2}); title('Mode-2'); xlabel('Latent dimensions');
                if N>=3, subplot(2,3,3); hintonDiagram(Z{3}); title('Mode-3'); end
            case 3,
                set(0,'CurrentFigure',h1);
                subplot(2,3,1); imagesc(Z{1}); title('Mode-1'); ylabel('Length of #-mode');
                subplot(2,3,2); imagesc(Z{2}); title('Mode-2'); xlabel('Latent dimensions');
                if N>=3, subplot(2,3,3); imagesc(Z{3}); title('Mode-3'); end
        end
        subplot(2,3,4); bar(gammas); title('Posterior mean of \lambda'); xlabel('Latent components'); ylabel(''); axis tight;
        subplot(2,3,5); plot(LB, '-r.','LineWidth',1.5,'MarkerSize',10 ); title('Lower bound'); xlabel('Iteration');  grid on;
        subplot(2,3,6); plotGamma(a_betaN, b_betaN); title('Posterior pdf'); xlabel('Noise precision \tau');grid on;
        set(findall(h1,'type','text'),'fontSize',12);
        drawnow;
    end
    if verbose==2
        set(0,'CurrentFigure',h3);
        figure(h3);
        %        temp = (X-min(X(:)))/(max(X(:))-min(X(:)));
        %        image(temp);
        %        axis off;
        imshow(uint8(X));
        title(['Iter.= '  num2str(it),',  Rank = ' num2str(max(rankest))],'FontSize', 13, 'color','b');
        xlabel(['(BCPF-MP)  Time: ' num2str(round(toc)), ' seconds'],'FontSize', 13, 'color','b');
        drawnow;
    end

    %% Convergence check
    if it>5 && abs(LBRelChan) < tol
        disp('\\\======= Converged===========\\\');
        break;
    end
end

%% Predictive distribution
if predVar==1
    Xvar =  tenzeros(size(Y));
    for n=1:N
        Xvar = tenmat(Xvar,n);
        Fslash = khatrirao_fast(Z{[1:n-1, n+1:N]},'r');
        if 1
            temp1 = double(tenmat(tensor(ZSigma{n}),3));
            temp2 = khatrirao_fast(Fslash', Fslash');
            Xvar(:,:) = Xvar(:,:) + temp1*temp2;
        else
            % ---  slow computation ------
            for i=1:size(Xvar,1)     %#ok
                Xvar(i,:) = Xvar(i,:) + diag(Fslash * ZSigma{n}(:,:,i) *Fslash')';
            end
            % ---  slow computation ------
        end
        Xvar = tensor(Xvar);
    end
    Xvar = Xvar + beta^(-1);
    Xvar = Xvar.*(2*a_betaN)/(2*a_betaN-2);
    Xvar = Xvar.*(dscale^2);
else
    Xvar =[];
end

%% Prepare the results
SNR = 10*log10(var(X(:))*beta);
X = ktensor(Z)*dscale;
X = arrange(X);

%% Output
model.X = X;
model.ZSigma = ZSigma;
model.gammas = gammas;
model.Fit = Fit;
model.SNR = SNR;
model.Xvar = double(Xvar);
model.TrueRank = rankest;
model.LowBound = max(LB);

function y = safelog(x)
x(x<1e-300)=1e-200;
x(x>1e300)=1e300;
y=log(x);


text:
function training_data = balanceReactionsInTrainingData(training_data)

if ~isfield(training_data, 'Ematrix') || isempty(training_data.Ematrix)
    [MW, Ematrix] = getMolecularWeight(training_data.nstd_inchi, 0);
    training_data.Ematrix = Ematrix(:, 2:end); % remove H, columns are [C, N, O, P, S, e-]
    conserved = training_data.Ematrix' * training_data.S;

    % need to check that all elements are balanced (except H, but including e-)
    % if only O is not balanced, add water molecules

    % check all reactions which can be checked (not NaN) and should be checked
    % (i.e. not formation or redox reactions)
    inds = find(~isnan(conserved(1,:)) .* training_data.balance');

    % first add water molecules to reactions that need it
    i_h2o = find(training_data.cids == 1);
    training_data.S(i_h2o, inds) = training_data.S(i_h2o, inds) - conserved(3, inds);

    % recalculate conservation matrix
    conserved = training_data.Ematrix' * training_data.S;

    inds_to_remove = inds(find(any(conserved(:, inds))));

    inds = setdiff(1:size(training_data.S, 2), inds_to_remove);
    training_data.S = training_data.S(:, inds);
    training_data.dG0_prime = training_data.dG0_prime(inds);
    training_data.T = training_data.T(inds);
    training_data.I = training_data.I(inds);
    training_data.pH = training_data.pH(inds);
    training_data.pMg = training_data.pMg(inds);
    training_data.weights = training_data.weights(inds);
    training_data.balance = false(size(inds));
end

fprintf('Successfully created balanced training-data structure: %d compounds and %d reactions\n',...
        size(training_data.S, 1), size(training_data.S, 2));

author: opencobra

repo: cobratoolbox


text:
function linpack_s_test16 ( )

%*****************************************************************************80
%
%% TEST16 tests SPOCO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 5;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST16\n' );
  fprintf ( 1, '  For a positive definite symmetric matrix,\n' );
  fprintf ( 1, '  SPOCO estimates the reciprocal condition number.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Set the matrix A.
%
  a(1:n,1:n) = 0.0;

  for i = 1 : n
    a(i,i) = 2.0;
    if ( 1 < i )
      a(i,i-1) = -1.0;
    end
    if ( i < n )
      a(i,i+1) = -1.0;
    end
  end
%
%  Estimate the condition.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimate the condition.\n' );

  [ a, rcond, z, info ] = spoco ( a, lda, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reciprocal condition  = %f\n', rcond );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function h = huniform ( p, varargin )

%*****************************************************************************80
%
%% HUNIFORM returns a uniform mesh size function.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real P(NP,ND), the point coordinates.
%
%    Input, VARARGIN, room for extra arguments.
%
%    Output, real H(NP,1), the mesh size function.
%
  np = size ( p, 1 );
  h = ones ( np, 1 );

  return
end


text:
function[varargout] = vsum(varargin)
%VSUM  Sum over non-NaN elements along a specified dimension.
%
%   Y=VSUM(X,DIM) takes the sum of all non-NaN elements of X along
%   dimension DIM.
%
%   [Y,NUM]=VSUM(X,DIM) also outputs the number of non-NaN data points NUM,
%   which has the same dimension as Y.
%
%   [Y1,Y2,...YN]=VSUM(X1,X2,...XN,DIM) also works.
%
%   VSUM(X1,X2,...XN,DIM);  with no output arguments overwrites the
%   original input variables.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2001--2016 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
  vsum_test,return
end

dim=varargin{end};

for i=1:length(varargin)-1
  [varargout{i},numi{i}]=vsum1(varargin{i},dim);
end

for i=length(varargin):nargout
  varargout{i}=numi{i-length(varargin)+1};
end

eval(to_overwrite(nargin-1))

function[y,num]=vsum1(x,dim)
ver=version;
%if str2num(ver(1)) >=9
try
    x=vswap(x,inf,nan);
    x=vswap(x,-inf,nan);
    x=vswap(x,inf+1i*inf,nan);

    y=sum(x,dim,'omitnan');
    if nargout==2
        num=sum(isfinite(x),dim,'omitnan');
    end
%else
catch
    y=sum(x,dim);
    if allall(isfinite(y))
        num=size(x,dim)+zeros(size(y));
    else
        nani=~isfinite(x);
        x(nani)=0;

        y=sum(x,dim);
        num=sum(~nani,dim);

        y(num==0)=nan;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[]=vsum_test
x1=[1 2 ; nan 4];
x2=[nan 6; nan 5];
ans1=[3 4]';
ans2=[6 5]';

vsum(x1,x2,2);
reporttest('VSUM output overwrite', aresame(x1,ans1) && aresame(x2,ans2))

x1=[1 2 ; nan 4];
ans1=[3 4]';
ans2=[2 1]';

[y1,y2]=vsum(x1,2);
reporttest('VSUM sum & num', aresame(y1,ans1) && aresame(y2,ans2))

x1=[1 2 ; 0 4];
ans1=[3 4]';
ans2=[2 2]';

[y1,y2]=vsum(x1,2);
reporttest('VSUM sum & num, no NaNs', aresame(y1,ans1) && aresame(y2,ans2))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [MVAbase, bus, gen, gencost, branch, f, success, et] = ...
                runuopf(casedata, mpopt, fname, solvedcase)
%RUNUOPF  Runs an optimal power flow with unit-decommitment heuristic.
%   [RESULTS, SUCCESS] = RUNUOPF(CASEDATA, MPOPT, FNAME, SOLVEDCASE)
%
%   Runs an optimal power flow (AC OPF by default) with a heuristic which
%   allows it to shut down "expensive" generators, optionally returning
%   a RESULTS struct and SUCCESS flag.
%
%   Inputs (all are optional):
%       CASEDATA : either a MATPOWER case struct or a string containing
%           the name of the file with the case data (default is 'case9')
%           (see also CASEFORMAT and LOADCASE)
%       MPOPT : MATPOWER options struct to override default options
%           can be used to specify the solution algorithm, output options
%           termination tolerances, and more (see also MPOPTION).
%       FNAME : name of a file to which the pretty-printed output will
%           be appended
%       SOLVEDCASE : name of file to which the solved case will be saved
%           in MATPOWER case format (M-file will be assumed unless the
%           specified name ends with '.mat')
%
%   Outputs (all are optional):
%       RESULTS : results struct, with the following fields:
%           (all fields from the input MATPOWER case, i.e. bus, branch,
%               gen, etc., but with solved voltages, power flows, etc.)
%           order - info used in external <-> internal data conversion
%           et - elapsed time in seconds
%           success - success flag, 1 = succeeded, 0 = failed
%           (additional OPF fields, see OPF for details)
%       SUCCESS : the success flag can additionally be returned as
%           a second output argument
%
%   Calling syntax options:
%       results = runuopf;
%       results = runuopf(casedata);
%       results = runuopf(casedata, mpopt);
%       results = runuopf(casedata, mpopt, fname);
%       results = runuopf(casedata, mpopt, fname, solvedcase);
%       [results, success] = runuopf(...);
%
%       Alternatively, for compatibility with previous versions of MATPOWER,
%       some of the results can be returned as individual output arguments:
%
%       [baseMVA, bus, gen, gencost, branch, f, success, et] = runuopf(...);
%
%   Example:
%       results = runuopf('case30');
%
%   See also RUNOPF, RUNDUOPF.

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

%%-----  initialize  -----
%% default arguments
if nargin < 4
    solvedcase = '';                %% don't save solved case
    if nargin < 3
        fname = '';                 %% don't print results to a file
        if nargin < 2
            mpopt = mpoption;       %% use default options
            if nargin < 1
                casedata = 'case9'; %% default data file is 'case9.m'
            end
        end
    end
end

%%-----  run the unit de-commitment / optimal power flow  -----
[r, success] = uopf(casedata, mpopt);

%%-----  output results  -----
if fname
    [fd, msg] = fopen(fname, 'at');
    if fd == -1
        error(msg);
    else
        if mpopt.out.all == 0
            printpf(r, fd, mpoption(mpopt, 'out.all', -1));
        else
            printpf(r, fd, mpopt);
        end
        fclose(fd);
    end
end
printpf(r, 1, mpopt);

%% save solved case
if solvedcase
    savecase(solvedcase, r);
end

if nargout == 1 || nargout == 2
    MVAbase = r;
    bus = success;
elseif nargout > 2
    [MVAbase, bus, gen, gencost, branch, f, et] = ...
        deal(r.baseMVA, r.bus, r.gen, r.gencost, r.branch, r.f, r.et);
% else  %% don't define MVAbase, so it doesn't print anything
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [output, layerOut] = FeatureTree_obj(visible, para, layer)
% remove target and cost function nodes if any
for i=1:length(layer)
    if strcmpi(layer{i}.name, 'target') || strcmpi(layer{i}.name, 'cross_entropy') || strcmpi(layer{i}.name, 'mse') || strcmpi(layer{i}.name, 'logistic')
        layer{i}.name = 'ignore';
    end
    if i > max(para.out_layer_idx)  % we don't compute layers that is after the last output layer
        layer{i}.name = 'ignore';
    end
end
para.NET.sentenceMinibatch = 1;
visible = visible.ShuffleData(para);
randomOrder = 0;

output = {}; layerOut = {};
for blk_i = 1:visible.nBlock
    minibatch = visible.PrepareMinibatch(para.precision, 1, para.NET.batchSize, blk_i,randomOrder);
    nMinibatch = size(minibatch,2);

    for utt_i = 1:nMinibatch
        PrintProgress(utt_i, nMinibatch, 100);
        batch_data = GetMinibatch(minibatch, utt_i, para.useGPU);

        % Use mode=3 to generate network output only
        if nargout>1
            [~, layerOut{end+1}, output{end+1}] = DNN_Cost10(layer, batch_data, para,3);
        else
            [~, ~, output{end+1}] = DNN_Cost10(layer, batch_data, para,3);
        end
    end
end
end

author: singaxiong

repo: SignalGraph


text:
function PrepareTestData()

global param;

sceneFolder = param.testScenes;
outputFolder = param.testData;

[sceneNames, scenePaths, numScenes] = GetFolderContent(sceneFolder, [], true);

MakeDir(outputFolder);

for i = 1 : numScenes

    count = fprintf('Started working on scene %d of %d', i, numScenes);

    %%% reading input data
    curExpo = ReadExpoTimes(scenePaths{i});
    [curImgsLDR, curLabel] = ReadImages(scenePaths{i});

    %%% processing data
    [inputs, label] = ComputeTestExamples(curImgsLDR, curExpo, curLabel);

    %%% writing data
%     label=[];
%     inputs = ((curImgsLDR{1}));
%     inputs = cat(3, inputs, (curImgsLDR{2}));
%     inputs = cat(3, inputs, (curImgsLDR{3}));
    WriteTestExamples(inputs, label, [outputFolder, '/', sceneNames{i}, '.h5']);
%     WriteTestExamples(inputs, label, [outputFolder, '/', 'TrainSequence', num2str(i) '.h5']);

    fprintf(repmat('\b', [1, count]));
end

fprintf('Done\n\n');

author: qingsenyangit

repo: AHDRNet


text:
function[varargout]=vdiff(varargin)
%VDIFF	Length-preserving first central difference.
%
%   DX=VDIFF(X,DIM) differentiates X along dimension DIM using the first
%   central difference; DX is the same size as X.
%
%   [D1,D2,...,DN]=VDIFF(X1,X2,...,XN,DIM) for multiple input variables
%   also works.
%
%   VDIFF(X1,X2,...,DIM); with no output arguments overwrites the
%   original input variables.
%
%   DXDT=VDIFF(DT,...) optionally uses scalar timestep DT to approximate
%   a time derivative, i.e. DXDT equals DX divided by DT.
%   _____________________________________________________________________
%
%   First and last points
%
%   The first and last points must be treated differently, as the central
%   difference is not defined there.  Three different methods can be used.
%
%   VDIFF(...,STR) specifies which method to use.
%
%        'endpoint'  uses the first forwards / first backwards difference
%                    at the first and last point, respectively.
%        'periodic'  treats the array as being periodic along dimension DIM,
%                    so that the central difference is defined at endpoints.
%        'nans'      fills in the first and last values with NANs.
%
%   The default behavior is 'endpoint'.
%   _____________________________________________________________________
%
%   'vdiff --t' runs some tests.
%
%   Usage:  x=vdiff(x,dim);
%           x=vdiff(dt,x,dim);
%           x=vdiff(dt,x,dim,'periodic');
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2000--2015 J.M. Lilly --- type 'help jlab_license' for details

%   I am so irritated by diff

if strcmpi(varargin{1}, '--t')
  vdiff_test,return
end

if length(varargin{1})==1
    dt=varargin{1};
    varargin=varargin(2:end);
else
    dt=1;
end

if ischar(varargin{end})
    str=varargin{end};
    varargin=varargin(1:end-1);
else
    str='endpoints';
end

n=varargin{end};
varargin=varargin(1:end-1);

for i=1:length(varargin)
    varargout{i}=vdiff1(varargin{i},n,str);
    if dt~=1
       varargout{i}=varargout{i}/dt;
    end
end

eval(to_overwrite(length(varargin)))

function[y]=vdiff1(x,n,str)

if ~isempty(x)
    %y=vshift(x,1,n)/2-vshift(x,-1,n)/2;

    %/*********************************************************************
    %This block is just a slightly faster way to take the first central difference
    sizex=size(x);
    %Correct for the fact that sizex ignores singleton dimensions
    if n>length(sizex)
        for i=length(sizex)+1:n
            sizex(i)=1;
        end
    end
    sizey=sizex;
    sizey(n)=sizex(n)-2;

%     dx=diff(x,1,n)/2;
%     array=zeros(length(size(x)),1);
%     array(n)=-1;
%     y=vindex(dx,1:size(x,n)-2,n)+vindex(dx,2:size(x,n)-1,n);
%
    y=zeros(sizey);
    i1=1:2:size(x,n);
    if iseven(size(x,n))
        i2=i1+1;
    else
        i2=i1(1:end-1)+1;
    end

    %This is just the same as .... y=vshift(x,1,n)/2-vshift(x,-1,n)/2;
    y=vindexinto(y,diff(vindex(x,i1,n),1,n)/2,i1(1:end-1),n)+vindexinto(y,diff(vindex(x,i2,n),1,n)/2,i2(1:end-1),n);
    %\*********************************************************************

    if strcmpi(str(1:3),'end')
        dx1=vindex(x,2,n)-vindex(x,1,n);
        dxend=vindex(x,size(x,n),n)-vindex(x,size(x,n)-1,n);
        y=cat(n,dx1,cat(n,y,dxend));
    elseif strcmpi(str(1:3),'nan')
        dx1=nan*vindex(x,1,n);
        y=cat(n,dx1,cat(n,y,dx1));
    elseif strcmpi(str(1:3),'per')
        dx1=(vindex(x,2,n)-vindex(x,1,n))/2+(vindex(x,1,n)-vindex(x,size(x,n),n))/2;
        dxend=(vindex(x,size(x,n),n)-vindex(x,size(x,n)-1,n))/2+(vindex(x,1,n)-vindex(x,size(x,n),n))/2;
        y=cat(n,dx1,cat(n,y,dxend));
    end
else
    y=[];
end

function[]=vdiff_test

y1=(1:4)';
y2=2*(1:4)';
[x1,x2]=vdiff(y1,y2,1);
bool=aresame(x1,[1 1 1 1]').*aresame(x2,2*[1 1 1 1]');
reporttest('VDIFF', bool)
vdiff(y1,y2,1);
bool=aresame(y1,[1 1 1 1]').*aresame(y2,2*[1 1 1 1]');
reporttest('VDIFF output overwrite', bool)

dt=pi;
y1=(1:4)';
y2=2*(1:4)';
[x1,x2]=vdiff(pi,y1,y2,1);
bool=aresame(x1,[1 1 1 1]'./dt).*aresame(x2,2*[1 1 1 1]'./dt,1e-10);
reporttest('VDIFF with non-unit time step', bool)

rng(1);
xxx=randn(100000,2);
tic;y=vshift(xxx,1,1)/2-vshift(xxx,-1,1)/2;etime1=toc;
tic;y2=vdiff(xxx,1,'periodic');etime2=toc;
reporttest('VDIFF new and former algorithms match, row dimension, even length', aresame(y,y2,1e-10))
disp(['New algorithm was ' num2str(etime2./etime1) ' times faster than former algorithm.'])

xxx=randn(2,100000);
tic;y=vshift(xxx,1,2)/2-vshift(xxx,-1,2)/2;etime1=toc;
tic;y2=vdiff(xxx,2,'periodic');etime2=toc;
reporttest('VDIFF new and former algorithms match, column dimension, even length', aresame(y,y2,1e-10))
disp(['New algorithm was ' num2str(etime2./etime1) ' times faster than former algorithm.'])

rng(1);
xxx=randn(100001,2);
tic;y=vshift(xxx,1,1)/2-vshift(xxx,-1,1)/2;etime1=toc;
tic;y2=vdiff(xxx,1,'periodic');etime2=toc;
reporttest('VDIFF new and former algorithms match, row dimension, odd length', aresame(y,y2,1e-10))
disp(['New algorithm was ' num2str(etime2./etime1) ' times faster than former algorithm.'])

xxx=randn(2,100001);
tic;y=vshift(xxx,1,2)/2-vshift(xxx,-1,2)/2;etime1=toc;
tic;y2=vdiff(xxx,2,'periodic');etime2=toc;
reporttest('VDIFF new and former algorithms match, column dimension, odd length', aresame(y,y2,1e-10))
disp(['New algorithm was ' num2str(etime2./etime1) ' times faster than former algorithm.'])


text:
function opt = mpoption_info_glpk(selector)
%MPOPTION_INFO_GLPK  Returns MATPOWER option info for GLPK.
%
%   DEFAULT_OPTS = MPOPTION_INFO_GLPK('D')
%   VALID_OPTS   = MPOPTION_INFO_GLPK('V')
%   EXCEPTIONS   = MPOPTION_INFO_GLPK('E')
%
%   Returns a structure for GLPK options for MATPOWER containing ...
%   (1) default options,
%   (2) valid options, or
%   (3) NESTED_STRUCT_COPY exceptions for setting options
%   ... depending on the value of the input argument.
%
%   This function is used by MPOPTION to set default options, check validity
%   of option names or modify option setting/copying behavior for this
%   subset of optional MATPOWER options.
%
%   See also MPOPTION.

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

if nargin < 1
    selector = 'D';
end
if have_feature('glpk')
    switch upper(selector)
        case {'D', 'V'}     %% default and valid options
            opt = struct(...
                'glpk',           struct(...
                    'opts',            [], ...
                    'opt_fname',       '' ...
                ) ...
            );
        case 'E'            %% exceptions used by nested_struct_copy() for applying
            opt = struct(...
                'name',         { 'glpk.opts' }, ...
                'check',        0 ...
                );
%                 'copy_mode',    { @glpk_options } ...
        otherwise
            error('mpoption_info_glpk: ''%s'' is not a valid input argument', selector);
    end
else
    opt = struct([]);       %% GLPK is not available
end


text:
function stimList = getRandom(stimList)
% Randomizes rows of a matrix, preserving dependencies between columns
%
% :Usage:
% ::
%
%     stimList = getRandom(stimList)
%
% :Input:
%   a col. vector or matrix of stimulus conditions,
%
%   e.g. [1 1 1 1 2 2 2 2 3 3 4 4]'
%
%output: a randomized permutation of this vector or matrix
% all columns are resorted with the same order
%
% ..
%    Tor Wager, created 2 / 01, last modified 2/25/04 to sort whole matrix
%    Modified 8/7/2015 by Tor to speed and simplify code
% ..

n = size(stimList, 1);

wh = randperm(n)';

stimList = stimList(wh, :);

% Old
%[n, k] = size(stimList);
% randvector = rand(n, 1);						% create random vector
% stimList(:, k + 1) = randvector;
% stimList = sortrows(stimList, k);				% sort the rows by random seed
% stimList = stimList(:, 1:end-1);

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function label = DBSCAN(data, eps, minPts)
% label = DBSCAN(data, epsilon, minPts)
%   Main part of DBSCAN (Density-Based Spatial Clustering of Application with
% Noise)clustering algorithm.
% ----
% Args:
%   data: data to be clustered (n * p)
%   eps: distance threshold for finding neighbors
%   minPts: minimum required number of neighbor observations for one core object
% ----
% Returns:
%   label: corresponding class for each observation

% Initialization
core_object = [];
count = 1;
k = 1; % label index
vmask = ones(size(data, 1), 1); % variable used to record whether one observation has been visited (1: not visited, 0: visited)
label = zeros(size(data, 1), 1); % pre-allocate result

% Find core objects
for i = 1:size(data, 1)
    x = data(i, :); % current observation of interests
    [~, neighbors] = Find_Neighbor(data, x, eps);

    % Check if the number of observations in neighbors of x is larger than minPts
    if((size(neighbors,1)-1) >= minPts)

        % Add the index of x to the set of core object
        core_object(count, 1) = i;
        count = count + 1;
    end
end

fprintf('---- DBSCAN finds a total number of %i core objects ----\n', size(core_object, 1));

while(~isempty(core_object))

    % Construct one queue for further use
    queue = [];
    qcount = 1;

    % Update core_object, observations have been visited should be deleted
    tcount = 1;
    list = [];
    for i = 1:size(core_object, 1)
        if(vmask(core_object(i, 1),1) == 0)
            list(tcount ,1) = i;
            tcount = tcount + 1;
        end
    end
    core_object(list,:) = [];

    if(isempty(core_object))
        break;
    end

    % Randomly choose one core object
    index = core_object(randi(size(core_object, 1)), 1);
    queue(qcount, 1) = index; % add the selected core object to the queue
    qcount = qcount + 1;
    vmask(index, 1) = 0; % update vmask
    label(index, 1) = k; % assign result
    core_object(core_object == index, :) = []; % remove the selected core object from the set

    while(~isempty(queue))
        pivot = queue(1, 1);
        queue(1, :) = []; % remove the first element in queue
        qcount = qcount - 1;

        % Find neighbors for current observation of interests
        [tindex, tneighbors] = Find_Neighbor(data, data(pivot, :), eps);
        if(size(tneighbors, 1) >= minPts)
            for i = 1 : size(tindex, 1)
                % If one neighbor has not been visited, add it to the queue
                if(vmask(tindex(i, 1),1) == 1)
                    queue(qcount, 1) = tindex(i, 1);
                    qcount = qcount + 1;
                    vmask(tindex(i, 1),1) = 0; % update vmask
                    label(tindex(i, 1),1) = k; % assign result
                end
            end
        end
    end

    fprintf('---- Finish finding observations for class %i ----\n', k);
    k = k + 1; % next class
end
end

function d = Eculidean_Distance(x1, x2)
% Calculate eculidean distance between x1 and x2.
%
% Args:
%   x1: observation 1
%   x2: observation 2
%
% Returns:
%   d: Eculidean distance between x1 and x2

d = sqrt((x1 - x2) * (x1 - x2)');
end

function [index, neighbors] = Find_Neighbor(data, x, eps)
% Find neighbors for x, where neighbors are defined to be observations in data
% with eculidean distance from x smaller than epsilon.
%
% Args:
%   data: dataset
%   x: current observation of interests
%   eps: pre-defined distance threshold
%
% Returns:
%   index: the index of these neighbors in data
%   neighbors: neighbor set of x

neighbors = []; % pre-allocate neighbor result
count = 1;
for i = 1 : size(data,1)
    d = Eculidean_Distance(data(i,:), x);
    if(d <= eps)
        index(count, 1) = i;
        neighbors(count,:) = data(i,:);
        count = count + 1;
    end
end
end

author: xuyxu

repo: Clustering

save_path: github-repos/MATLAB/xuyxu-Clustering


text:
function [HotellingT2] = HotellingT2(X,alpha)
%Hotelling T-Squared testing procedures for multivariate samples.
%
%   Syntax: function [HotellingT2] = HotellingT2(X,alpha)
%
%     Inputs:
%          X - multivariate data matrix.
%      alpha - significance level (default = 0.05).
%
%     Outputs:
%          It depends of the Hotelling's T-Squared multivariate test of interest,
%          being able to be:
%
%            |-One-sample
%            |                          |-Homoskedasticity (to test)
%            |            |-Independent |
%            |            |             |-Heteroskedasticity (to test)
%            |-Two-sample |
%                         |
%                         |-Dependent
%
%          Each case calls to a corresponding function that contains a complete
%          explanation.
%
%  Created by A. Trujillo-Ortiz and R. Hernandez-Walls
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             atrujo@uabc.mx
%             And the special collaboration of the post-graduate students of the 2002:2
%             Multivariate Statistics Course: Karel Castro-Morales, Alejandro Espinoza-Tenorio,
%             Andrea Guia-Ramirez.
%
%  Copyright (C) December 2002.
%

if nargin < 2,
    alpha = 0.05; %(default)
end;

if nargin < 1,
   error('Requires at least one input argument.');
end;

sam = input('Do you have one multivariate sample (1) or two multivariate samples (2)?: ');
if sam == 1;
   T2Hot1(X,alpha)
else
   id = input('They are independent (1) or dependent (2)?: ');
   if id == 1;
      disp('The covariance matrix homogeneity will be testing.:');
      MBoxtest(X,alpha);
      disp(' ')
      dc = input('Do they were significant? (y/n): ','s');
      if dc == 'y'
         T2Hot2ihe(X,alpha);
      else
         T2Hot2iho(X,alpha);
      end;
   else
      T2Hot2d(X,alpha);
   end;
end;

return;

author: Sable

repo: mcbench-benchmarks


text:
function psat_values_test ( )

%*****************************************************************************80
%
%% PSAT_VALUES_TEST demonstrates the use of PSAT_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PSAT_VALUES_TEST:\n' );
  fprintf ( 1, '  PSAT_VALUES stores values of\n' );
  fprintf ( 1, '  the saturation pressure of water\n' );
  fprintf ( 1, '  as a function of temperature.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      T            PSAT(T)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, tc, psat ] = psat_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', tc, psat );

  end

  return
end


text:
%EUL2TR Convert Euler angles to homogeneous transform
%
% T = EUL2TR(PHI, THETA, PSI, OPTIONS) is an SE(3) homogeneous
% transformation matrix (4x4) with zero translation and rotation equivalent
% to the specified Euler angles. These correspond to rotations about the Z,
% Y, Z axes respectively. If PHI, THETA, PSI are column vectors (Nx1) then
% they are assumed to represent a trajectory and R is a three-dimensional
% matrix (4x4xN), where the last index corresponds to rows of PHI, THETA,
% PSI.
%
% R = EUL2R(EUL, OPTIONS) as above but the Euler angles are taken from the
% vector (1x3)  EUL = [PHI THETA PSI]. If EUL is a matrix (Nx3) then R is a
% three-dimensional matrix (4x4xN), where the last index corresponds to
% rows of RPY which are assumed to be [PHI,THETA,PSI].
%
% Options::
%  'deg'      Angles given in degrees (radians default)
%
% Note::
% - The vectors PHI, THETA, PSI must be of the same length.
% - The translational part is zero.
%
% See also EUL2R, RPY2TR, TR2EUL, SE3.eul.

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

function T = eul2tr(phi, varargin)

    R = eul2r(phi, varargin{:});
    T = r2t(R);

author: petercorke

repo: spatialmath-matlab


text:
function J=polarU2DCrossGrad(uList,systemType)
%%POLARU2DCROSSGRAD Given the direction cosine value u in 2D, obtain the
%              derivative of the polar azimuth angle with respect to u.
%
%INPUTS: uList A 1XnumPoints (for only u) or a 2XnumPoints (if full unit
%              vectors are given) set of direction cosines in 2D.
%   systemType An optional parameter specifying the axis from which the
%              angles are measured. Possible values are
%              0 (The default if omitted or an empty matrix is passed) The
%                azimuth angle is counterclockwise from the x axis.
%              1 The azimuth angle is measured clockwise from the y axis.
%
%OUTPUTS: J A 1XnumPoints set of derivatives of the azimuth angle with
%           respect to u evaluated at the given u values.
%
%EXAMPLE:
%Here, we verify that the derivatives returned by this function are about
%equal to those returned via numeric differentiation (forward
%differencing).
% points=[0.1,0.2,-0.2,0,-0.9];%u points
% systemType=1;
% epsVal=1e-9;
%
% az=u2PolAng2D(points,systemType);
% az1=u2PolAng2D(points+epsVal,systemType);
% JNumDiff=(az1-az)/epsVal;
% J=polarU2DCrossGrad(points,systemType);
%
% max(abs(JNumDiff-J))
%One will see that the difference is on the order of 3e-8, which is a good
%agreement.
%
%June 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(systemType))
    systemType=0;
end

hasV=size(uList,1)>1;

N=size(uList,2);

J=zeros(1,N);
for curPoint=1:N
    u=uList(1,curPoint);

    if(hasV)
        v=uList(2,curPoint);
    else
        v=sqrt(1-u^2);
    end

    switch(systemType)
        case 0
            J(curPoint)=-1/v;
        case 1
            J(curPoint)=1/v;
        otherwise
            error('Invalid system type specified.')
    end
end

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [feature_pos_out, transformed_points_out, A_out, error] = em2_affine(...
    events, ...
    feature_pos, ...
    template_points, ...
    template_weights, ...
    A_in, ...
    flow, ...
    params, ...
    fig)
%EM2_AFFINE Estimates an affine transform between a set of events and a set
%of template points.
%
% EM2_AFFINE estimates an affine warp (A, b) between a set of events with
% known optical flow, and a set of 2D template points. This EM step runs
% from the events alone, as in:
% Alex Zihao Zhu, Nikolay Atanasov and Kostas Daniilidis.
% "Event-based Feature Tracking with Probabilistic Data Association",
% IEEE International Conference on Robotics and Automation (ICRA), 2017.
%
% Syntax:  [feature_pos_out, transformed_points_out, A_out, error] = EM2_AFFINE(...
%            events, ...
%            feature_pos, ...
%            template_points, ...
%            template_weights, ...
%            A_in, ...
%            flow, ...
%            params, ...
%            fig)
%
% Inputs:
%    events           - 4xN, each column is (x,y,t,p).
%    feature_pos      - 2x1, pixel position of the feature.
%    template_points  - 2xM, set of points to align the events with.
%    template_weights - 1xM, weights for each template point.
%    A_in             - 2x2, initialization for the affine warp.
%    flow_init        - 2x1, flow for the event window.
%    params           - parameters, defined in get_params().
%    fig              - figure handle for plotting.
%
% Outputs:
%    feature_pos_out        - 2x1, feature pos after alignment, feature_pos - b.
%    transformed_points_out - 2xN, input events in the feature window
%                               shifted by the flow A[x; y] + dt * flow + b.
%    A_out                  - 2x2, estimate for the affine warp.
%    error                  - Error code, defined in TrackingErrors.
%
% See also GET_PARAMS
%
% Author: Alex Zihao Zhu, University of Pennsylvania
% Email: alexzhu(at)seas.upenn.edu
% Copyright 2018 University of Pennsylvania
% Alex Zihao Zhu, Nikolay Atanasov, Kostas Daniilidis
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS, CONTRIBUTORS, AND THE
% TRUSTEES OF THE UNIVERSITY OF PENNSYLVANIA "AS IS" AND ANY EXPRESS OR
% IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
% OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT OWNER, CONTRIBUTORS OR THE TRUSTEES OF
% THE UNIVERSITY OF PENNSYLVANIA 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.

%% Initialization
% Estimated canonical events
transformed_points_out = [];
% New estimated patch position if estimated
feature_pos_out = [nan; nan];
errs = zeros(params.em2_params.max_iters, 1);
num_iter = 0;
error = TrackingErrors.ICPFailed;
scatter_plot_handle = [];

target_time = events(3, 1);
A = A_in;
A_out = A;
pos_shift = [0;0];

event_window = [];

time_shifted_points = bsxfun(@minus, events(1:2, :), feature_pos) + ...
    bsxfun(@times, flow, (target_time - events(3,:)));

if isempty(template_points)
    return
end

kdtree = KDTreeSearcher(...
    template_points' / (sqrt(2) * params.em2_params.sigma), ...
    'Distance', 'euclidean');

while true
    if num_iter > params.em2_params.max_iters
        error = TrackingErrors.ICPMaxIters;
        return
    end

    transformed_points = bsxfun(@plus, A * time_shifted_points, pos_shift);

    if isempty(event_window)
        event_window = ...
            transformed_points(1, :) >= -params.window_size/2 & ...
            transformed_points(2, :) >= -params.window_size/2 & ...
            transformed_points(1, :) <= params.window_size/2 & ...
            transformed_points(2, :) <= params.window_size/2;

        n_in_window = sum(event_window);
        if n_in_window < params.min_events_for_em
            error = TrackingErrors.ICPFailed;
            return
        end

        time_shifted_points = time_shifted_points(:, event_window);
        transformed_points = transformed_points(:, event_window);
        %         if n_in_window > 200
        %             shifted_events = shifted_events(:, 1:200);
        %             transformed_events = transformed_events(:, 1:200);
        %             n_in_window = 200;
        %         end

        weights = zeros(size(template_points, 2), size(transformed_points, 2));
    end

    [neighbors_cell, distances_cell] = rangesearch(...
        kdtree, ...
        transformed_points' / (sqrt(2) * params.em2_params.sigma), ...
        params.em2_params.max_distance);

    distancesstacked = cell2mat(cellfun(...
        @transpose, distances_cell, 'UniformOutput', false))';

    if isempty(distancesstacked)
        error = TrackingErrors.ICPFailed;
        return
    end

    num_neighbors = cellfun('length', neighbors_cell);
    template_correspondences = cell2mat(cellfun(...
        @transpose, neighbors_cell, 'UniformOutput', false))';

    transformed_correspondences = repelem(1:n_in_window, num_neighbors);

    weightsstacked = exp(-distancesstacked);

    valid_inds = sub2indc(transformed_correspondences, template_correspondences, size(weights));

    % It's cheaper to multiply to 0 to reset the weights matrix than to
    % reinitialize it using zeros.
    weights = weights * 0;
    weights(valid_inds) = weightsstacked;

    if ~isempty(template_weights)
        weights = bsxfun(@times, template_weights, weights);
    end

    weight_sum = sum(weights, 1);
    valid_weights = (weight_sum > 0);

    weights = bsxfun(@rdivide, weights, weight_sum + 1e-10);

    weightsstacked = weights(valid_inds);

    %% Simultaneously minimize over affine warp and translation.
    weighted_template_points = template_points*weights;
    weighted_template_points = weighted_template_points(:, valid_weights);

    current_events = time_shifted_points(:, valid_weights);

    D = [current_events(1:2, :); ones(1, size(current_events, 2))];
    X = weighted_template_points;

    transform = X/D;

    A = transform(1:2, 1:2);
    scale = sqrt(abs(det(A)));
    pos_shift = (transform(:, 3));

    if scale < params.em2_params.min_scale
        error = TrackingErrors.ICPFailed;
        return
    end

    %% Calculate errors, plot debug information.
    err = weightsstacked*distancesstacked'/ size(current_events, 2);
    errs(num_iter+1) = err;

    if params.debug
        set(0, 'CurrentFigure', fig)
        subplot(2,1,1)
        plot(errs(errs > 0), 'b')
        title('Change in error (convergence criterion)')
        xlim([0 params.em2_params.max_iters])
        subplot(2,1,2)
        if (~isempty(scatter_plot_handle))
            delete(scatter_plot_handle)
        end
        scatter_plot_handle = scatter(transformed_points(1, :), transformed_points(2, :),'r.');
        hold on
        scatter(template_points(1, :), template_points(2, :),template_weights*10,'b.');
        hold off
        axis equal
        axis([-params.window_size/2-5, ...
            params.window_size/2+5, ...
            -params.window_size/2-5, ...
            params.window_size/2+5])
        axis ij
        title('EM 2')

        pause(0.01)
    end

    if num_iter
        derr = err - errs(num_iter);
        if derr > -params.em2_params.min_err && derr <= 0
            break;
        end
    end
    num_iter = num_iter + 1;
end

A_out = A;
% Update the feature position based on the shift.
feature_pos_out = feature_pos - pos_shift;

% Compute the template events to be used in later iterations.
dt = events(3, end)-events(3, 1);
time_shifted_points = bsxfun(@minus, events(1:2, :), ...
    feature_pos + flow * dt) + ...
    bsxfun(@times, flow, (events(3, end) - events(3, :)));

transformed_points = A_out*time_shifted_points;

% Make the window size a bit bigger for the template.
window_size = round(params.window_size * 1.5);

event_window = ...
    transformed_points(1, :) >= -window_size/2 & ...
    transformed_points(2, :) >= -window_size/2 & ...
    transformed_points(1, :) < window_size/2 & ...
    transformed_points(2, :) < window_size/2;

transformed_points_out = transformed_points(:, event_window);

if params.debug
    pause(0.5)
end

% Measure how many transformed points are considered 'outliers' when
% matched against the templates. Outliers are defined as points without any
% template points within 2 pixels. Note that this function still requires
% more work to be properly robust.
[~, distances_cell] = knnsearch(...
    kdtree, ...
    transformed_points_out'/(sqrt(2)*params.em2_params.sigma), ...
    'K', 1);
distances_cell = distances_cell * 2 * params.em2_params.sigma^2;
percent_outliers = sum(distances_cell > 2.0) / length(distances_cell);

if percent_outliers > params.em2_params.max_outlier_thresh
    error = TrackingErrors.ICPError;
    feature_pos_out = [nan; nan];
    return;
end

end

author: daniilidis-group


text:
function [F1c,F2c,F1,F2,C1,C2,VT]=truncatedOctahedronMesh(varargin)

% function [F1c,F2c,F1,F2,C1,C2,VT]=truncatedOctahedronMesh(r,nCopies)
% ------------------------------------------------------------------------
%
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2019/02/08 Created
% 2019/10/13
% ------------------------------------------------------------------------

%% Parse input

switch nargin
    case 1
        r=varargin{1};
        nCopies=2;
    case 2
        r=varargin{1};
        nCopies=varargin{2};
end

if isempty(nCopies)
    nCopies=2;
end

if numel(nCopies)==1
    nCopies=nCopies*ones(1,3);
end

%%
%Get truncated octahedron
[Fs,Vs,~]=truncatedOctahedron(r);
Fs1=Fs{1}; %hexagon faces
Fs2=Fs{2}; %quadrilateral faces

%Rotate
[R,~]=euler2DCM([0 0 0.25*pi]); %Define rotation
Vs=Vs*R; %Rotate coordinates

%Derive face centre points for offsets
% Vn=patchCentre(Fs,Vs);
Xi=Vs(:,1); Yi=Vs(:,2); Zi=Vs(:,3);
Vn=[mean(Xi(Fs2),2) mean(Yi(Fs2),2) mean(Zi(Fs2),2)];

%%
nTotal=prod(nCopies); %Total number of copies
offsetDirs=[1 2 6]; %Offset direction N.B. varying these affects the offsets/signs below

%Create cell indices for vertices
indC=ones(size(Vs,1),1)*(1:1:nTotal);
indC=indC(:);
[I,J,K] = ind2sub(nCopies,indC);
I=I-1; J=J-1; K=K-1;

%Create cell indices for quad faces
indF1=ones(size(Fs1,1),1)*(1:1:nTotal);
indF1=indF1(:);
indF1=(indF1-1);

indF2=ones(size(Fs2,1),1)*(1:1:nTotal);
indF2=indF2(:);
indF2=(indF2-1);

%Defining the quad faces matrix
F1=repmat(Fs1,nTotal,1)+size(Vs,1).*indF1(:,ones(1,size(Fs1,2)));
F2=repmat(Fs2,nTotal,1)+size(Vs,1).*indF2(:,ones(1,size(Fs2,2)));

%Defining offsets
sK=iseven(K); %Shift is adjusted according to z coordinate to create a "cube"
D1=I*2*Vn(offsetDirs(1),:)+(sK*Vn(offsetDirs(1),:)); % X offsets
D2=J*2*Vn(offsetDirs(2),:)+(sK*Vn(offsetDirs(2),:)); % Y offsets
D3=K*Vn(offsetDirs(3),:); % Z offsets

D=D1+D2+D3;

%Defining vertices matrix
VT=repmat(Vs,nTotal,1)+D;

%Merge points
[F1,VT,~,ind2]=mergeVertices(F1,VT);
F2=ind2(F2);

C1=indF1+1; %Index or color number
C2=indF2+1; %Index or color number

%Split up face matrix in to cell groups
F1c=mat2cell(F1,size(Fs1,1)*ones(1,nTotal),size(Fs1,2));
F2c=mat2cell(F2,size(Fs2,1)*ones(1,nTotal),size(Fs2,2));

end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% listdlg2() - listdlg function clone with coloring and help for
%               eeglab().
%
% Usage: same as listdlg()
%
% Author: Arnaud Delorme, CNL / Salk Institute, La Jolla, 16 August 2002
%
% See also: inputdlg2(), errordlg2(), supergui(), inputgui()

% Copyright (C) Arnaud Delorme, CNL / Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [vals, okornot, strval] = listdlg2(varargin);

if nargin < 2
   help listdlg2;
   return;
end;
for index = 1:length(varargin)
	if iscell(varargin{index}), varargin{index} = { varargin{index} }; end;
	if isstr(varargin{index}), varargin{index} = lower(varargin{index}); end;
end;
g = struct(varargin{:});

try,  g.promptstring;  catch, g.promptstring = ''; end;
try,  g.liststring;    catch, error('''liststring'' must be defined'); end;
try,  g.selectionmode; catch, g.selectionmode = 'multiple'; end;
try,  g.listsize;      catch, g.listsize = []; end;
try,  g.initialvalue;  catch, g.initialvalue = []; end;
try,  g.name;          catch, g.name = ''; end;

fig = figure('visible', 'off');
set(gcf, 'name', g.name);
if isstr(g.liststring)
    allstr =  g.liststring;
else
    allstr = '';
    for index = 1:length(g.liststring)
        allstr = [ allstr '|' g.liststring{index} ];
    end;
    allstr = allstr(2:end);
end;

geometry = {[1] [1 1]};
geomvert = [min(length(g.liststring), 10) 1];
if ~strcmpi(g.selectionmode, 'multiple') | ...
        (iscell(g.liststring) & length(g.liststring) == 1) | ...
        (isstr (g.liststring) & size  (g.liststring,1) == 1 & isempty(find(g.liststring == '|')))
	if isempty(g.initialvalue), g.initialvalue = 1; end;
    minval = 1;
	maxval = 1;
else
    minval = 0;
    maxval = 2;
end;
listui = {{ 'Style', 'listbox', 'tag', 'listboxvals', 'string', allstr, 'max', maxval, 'min', minval } ...
		  { 'Style', 'pushbutton', 'string', 'Cancel', 'callback', ['set(gcbf, ''userdata'', ''cancel'');'] }  ...
		  { 'Style', 'pushbutton', 'string', 'Ok'    , 'callback', ['set(gcbf, ''userdata'', ''ok'');'] } };

if ~isempty(g.promptstring)
	geometry = {[1] geometry{:}};
	geomvert = [1 geomvert];
	listui = { { 'Style', 'text', 'string', g.promptstring } listui{:}};
end;
[tmp tmp2 allobj] = supergui( fig, geometry, geomvert, listui{:} );

% assign value to listbox
% must be done after creating it
% ------------------------------
lstbox = findobj(fig, 'tag', 'listboxvals');
set(lstbox, 'value', g.initialvalue);

if ~isempty(g.listsize)
	pos = get(gcf, 'position');
	set(gcf, 'position', [ pos(1:2) g.listsize]);
end;
h = findobj( 'parent', fig, 'tag', 'listboxvals');

okornot = 0;
strval = '';
vals = [];
figure(fig);
drawnow;
waitfor( fig, 'userdata');
try,
	vals = get(h, 'value');
	strval = '';
    if iscell(g.liststring)
        for index = vals
            strval = [ strval ' ' g.liststring{index} ];
        end;
    else
        for index = vals
            strval = [ strval ' ' g.liststring(index,:) ];
        end;
    end;
	strval = strval(2:end);
	if strcmp(get(fig, 'userdata'), 'cancel')
		okornot = 0;
	else
		okornot = 1;
	end;
	close(fig);
    drawnow;
end;

author: PatternRecognition

repo: OpenBMI


text:
function [T, M, S] = dtiTsqTestStat(g1, g2, Y, mask)

% Computes voxel-wise Hotelling T^2 statistics for two groups from a data array of
% diffusion tensors in dt6 format.
%
%   [T, M, S] = dtiTsqTestStat(g1, g2, DT_ARRAY, [MASK])
%
% Input:
%   g1, g2      List of indices that correspond to each group out of 1:N
%                   E.g: g1 = 1:7, g2 = 8:14, N = 14
%   DT_ARRAY    Data array of size XxYxZxpxN (or nxpxN), where X, Y, Z are the volume
%                   dimensions and N is the number of subjects.
%                   p is the vector dimension (p = 6 for DT6 data).
%                   (n is the number of voxels).
%   MASK        Optional XxYxZ binary array. Values of M and S are computed
%                   where mask = 1; in other voxels, M and S are set to 0.
%                   Default is entire volume.
%
% Output:
%   T           XxYxZx1 array of test statistics (0 where mask = 0)
%   M           XxYxZxpx2 array of mean vectors for both groups (0 where mask = 0)
%   S           XxYxZxpxp array of covariances (0 where mask = 0)
%
% Utilities:    ndfun.m, dtiSplitTensor.m, dti33to6.m
%
% WARNING: If using Pentium 4, eliminate NaN's from array before running
% (processor bug).
%
% Copyright by Armin Schwartzman, 2005

% HISTORY:
%   2004.06.23 ASH (armins@stanford.edu) wrote it.
%

% Check inputs
if (ndims(Y)==2 | ndims(Y)==3),
    Ind = 1;    % Data in indexed nx6xN format
    Y = shiftdim(Y, -2);
else
    Ind = 0;    % Data in XxYxZx6xN format
end
if (ndims(Y)<4 | ndims(Y)>5),
    error('Wrong input format');
end
if (~exist('mask')),
    mask = ones([size(Y,1) size(Y,2) size(Y,3)]);
end

% Computations
N1 = length(g1);
N2 = length(g2);
N  = N1 + N2;

p = size(Y, 4);
Y = permute(Y, [4 5 1 2 3]); % permutation required by ndfun
Yavg1 = mean(Y(:,g1,:,:,:), 2);
Yavg2 = mean(Y(:,g2,:,:,:), 2);
davg = Yavg1 - Yavg2;
d1 = Y(:,g1,:,:,:)-repmat(Yavg1, [1 N1 1 1 1]);
d2 = Y(:,g2,:,:,:)-repmat(Yavg2, [1 N2 1 1 1]);
M = cat(2, Yavg1(:,1,:,:,:), Yavg2(:,1,:,:,:));
clear Y*
S = ndfun('mult', d1, permute(d1, [2 1 3:5])) + ndfun('mult', d2, permute(d2, [2 1 3:5]));
outmask = repmat(shiftdim(~mask, -2), [p p 1 1 1]);
S(outmask) = 0;
outmask = repmat(eye(p), [1 1 size(mask)]) & outmask;
S(outmask) = 1;
b = ndfun('backslash', S, davg);
b = ndfun('mult', permute(davg, [2 1 3:5]), b);
T = b * N1*N2/(N1+N2) * (N1+N2-p-1)/p;

% Adjust output
T = permute(T, [3 4 5 1 2]);
M = permute(M, [3 4 5 1 2]);
S = permute(S, [3 4 5 1 2]);
if Ind,
    T = shiftdim(T, 2);
    M = shiftdim(M, 2);
    S = shiftdim(S, 2);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% t_mrd
%
% Illustrates how to initiate mrDiffusion with scripts.
%
% There are several examples of loading data or adjusting the view by using
% the window handles and refreshing the screen.
%
% Additional mrDiffusion tutorial scripts are named t_mrd<XXXXX>.
%
% See also: t_mrdFibers, t_mrdTensor, t_mrdViewFibers, t_mrdDWI, dtiSet, dtiGet, dwiGet,
% dwiSet, dwiCreate
%
% (c) Stanford VISTA Team, 2011

%% You can initiate a mrDiffusion window with no data

% We use dtiF to store the figure number and dtiH to store the window
% handles.
[dtiF, dtiH] = mrDiffusion;

% Close the figure.  This  renders the dtiH
% (handles to objects in the window) dead.
close(dtiF);
clear dtiH

% You can also get the handles this way
dtiF = mrDiffusion;
dtiH = guidata(dtiF);

% When you make a change to the data in the figure, by adjusting the dtiH
% (window handles), you reset the data by the call
guidata(dtiF,dtiH);

% Let's start again
close(dtiF);

%% Vistadata contains a sample data set with a calculated dt6

% Load the sample data set and set the window to invisible
dataDir = fullfile(mrvDataRootPath,'diffusion','sampleData');
dt6Name = fullfile(dataDir,'dti40','dt6.mat');
[dtiF, dtiH] = mrDiffusion('off',dt6Name);

% Make the window visible
figure(dtiF)

% Normally, we interact with the diffusion data using dtiGet/Set or
% dwiGet/Set.  In this case, we have opened only a diffusion tensor imaging
% structure through this window.  Hence, we interact with the data using
% dtiGet.  We can find out which image is displayed in the background using
dtiGet(dtiH,' background name')
d = dtiGet(dtiH,' background data');
showMontage(d)

% Close the window, clear the handles
close(dtiF)

%% You can open the data and leave the window visible this way
[dtiF,dtiH] = mrDiffusion('on',dt6Name);

% Reset
close(dtiF)

%% Load some fibers

% Open the window
[dtiF,dtiH] = mrDiffusion('on',dt6Name);

% Load a sample fiber group
fgName = fullfile(mrvDataRootPath,'diffusion','sampleData','fibers','leftArcuate.pdb');
fg = mtrImportFibers(fgName);

% Attach the fiber group to the handles
dtiH = dtiSet(dtiH,'current fiber group',fg);

% We interact with the window handles using the Matlab set/get functions
% applied to the handles.
set(dtiH.cbShowFibers,'Value',1)
guidata(dtiF,dtiH)

% After setting and attaching, you refresh the window this way
dtiH = dtiRefreshFigure(dtiH);

%% Make the fibers invisible
set(dtiH.cbShowFibers,'Value',0)
guidata(dtiF,dtiH)
dtiH = dtiRefreshFigure(dtiH);

% See the script t_mrdViewFibers for examples of viewing the fibers with
% the mrMesh viewer

%% Read mrDiffusion figure handle settings

overlayThresh = get(dtiH.slider_overlayThresh, 'Value');
overlayAlpha = str2double(get(dtiH.editOverlayAlpha, 'String'));
curOvNum = get(dtiH.popupOverlay,'Value');

%% Find a mrDiffusion window
dtiF = dtiGet([],'main figure');
close(dtiF)

%% End


text:
function new_files = remove_disdaq_vols(img_files, num_vols_per_run, num_disdaq_vols, varargin)
% :Usage:
% ::
%
%     new_files = remove_disdaq_vols(img_files, num_vols_per_run, num_disdaq_vols, ['overwrite', 0|1], ['FSLOUTPUTTYPE', fsl_output_type], ['strict', 0|1])
%
% :Inputs:
%
%   **img_files:**
%        cellstr of files to remove disdaqs from - each cell represents a run
%
%   **num_vols_per_run:**
%        vector of volume counts per run, *not* including disdaq vols
%
%   **num_disdaq_vols:**
%        constant describing how many data points to remove from the beginning of each run
%
%   **'overwrite':**
%        if set, will overwriting images in place - defaults to 0
%
%   **'FSLOUTPUTTYPE':**
%        outputfile type - defaults to 'NIFTI'
%
%   **'strict':**
%        if set, will error out unless data given to it is
%        exactly proper length - defaults to 1 - turn off ONLY with good reason
%
% :Output:
%
%   **new_files:**
%        input files, but with a 'd' prepended, unless 'overwriting' was specified
%
% NB: Not set up for 3d files yet!!!
%
% :Examples:
% ::
%
%    % for an experiment
%    num_vols_per_run = [124 140 109];  % NOT including disdaqs
%    num_disdaq_vols = 4;

    global FSLDIR;
    scn_setup();

    % Flags
    strict = 1;
    overwriting = 0;
    return_char = 0;
    fsl_output_type = 'NIFTI';

    % Check inputs
    if ~isempty(varargin)
        for i = 1:length(varargin)
            if ischar(varargin{i})
                switch(varargin{i})
                    case {'overwrite' 'overwriting'}
                        overwriting = varargin{i+1};
                    case 'strict'
                        strict = varargin{i+1};
                    case {'fsl_output_type', 'FSLOUTPUTTYPE'}
                        fsl_output_type = varargin{i+1};
                end
            end
        end
    end

    if(ischar(img_files))
        disp('Cellstr not passed in, changing...')
        img_files = cellstr(img_files);
        return_char = 1;
    end

    % Use FSL to remove disdaqs
    new_files = cell(size(img_files));
    for i=1:length(img_files)
        current_file = img_files{i};
        num_frames = scn_num_volumes(current_file);

        if(strict && (num_frames ~= num_vols_per_run(i) + num_disdaq_vols(i)))
            error('The number of vols (%d) + disdaqs (%d) don''t add up to the length (%d) of the file (''%s'') passed in.', num_vols_per_run, num_disdaq_vols, num_frames, current_file);
        else
            if(overwriting)
                output_file = current_file;
            else
                [d f e] = fileparts(current_file);
                output_file = fullfile(d, ['d' f]);
            end

            fslroi_command = sprintf('export FSLDIR=%s && . %s/etc/fslconf/fsl.sh && FSLOUTPUTTYPE=%s && %s/bin/fslroi %s %s %d %d', ...
                FSLDIR, FSLDIR, fsl_output_type, FSLDIR, current_file, output_file, num_disdaq_vols, num_vols_per_run);
            disp(fslroi_command);
            system(fslroi_command);
        end

        new_files{i} = [output_file e];
    end

    % Convert output back to char is only a single file was passed in
    if(return_char)
        new_files = char(new_files);
    end
end


text:
function accuracy = rpca_src(TrainSet, TestSet, train_num, test_num, class_num, lambda, options)
% Robust PCA-based sparse representation classification (RobustPCA-SRC) algorithm
%
% Inputs:
%       TrainSet            train sets of size dxn, where d is dimension and n is number of sets
%       train_num           numner of train sets
%       class_num           numner of classes
%       lambda              regularization paramter
% Output:
%       accuracy            classification accurary
%
%
% Created by H.Kasai on July 07, 2017

    % extract options
    if ~isfield(options, 'verbose')
        verbose = false;
    else
        verbose = options.verbose;
    end

    if ~isfield(options, 'eigenface')
        eigenface = true;
    else
        eigenface = options.eigenface;
    end

    if ~isfield(options, 'eigenface_dim')
        eigenface_dim = train_num;
    else
        eigenface_dim = options.eigenface_dim;
    end

    if ~isfield(options, 'solver_max_iter')
        solver_max_iter = 100;
    else
        solver_max_iter = options.solver_max_iter;
    end

    % define problem definitions
    solver_lambda = 1/sqrt(max(size(TrainSet.X)));
    solver_lambda = solver_lambda/3;
    mask = logical(zeros(size(TrainSet.X)));
    problem = robust_pca(TrainSet.X, mask, solver_lambda);
    problem = robust_pca(TrainSet.X, mask, solver_lambda);

    % perform robust PCA
    solver_options.max_iter = solver_max_iter;
    solver_options.verbose = verbose;
    solver_options.mu = 10*solver_lambda/3;
    [w, ~] = admm_robust_pca(problem, solver_options);
    TrainSet.X = w.L;

    problem = robust_pca(TestSet.X, mask, solver_lambda);
    [w, ~] = admm_robust_pca(problem, solver_options);
    TestSet.X = w.L;

    % generate eigenface
    if eigenface
        [disc_set, ~, ~] = Eigenface_f(TrainSet.X, eigenface_dim);

        % project on subspace
        TrainSet.X  =  disc_set' * TrainSet.X;
        TestSet.X   =  disc_set' * TestSet.X;
    end

    % normalize data to l2-norm
    [TrainSet.X, ~] = data_normalization(TrainSet.X, TrainSet.y, 'std');
    [TestSet.X, ~] = data_normalization(TestSet.X, TestSet.y, 'std');

    % perform a standard SRC
    if 1
        % use src function
        options.eigenface = false;
        accuracy = src(TrainSet, TestSet, train_num, test_num, class_num, lambda, options);
    else

        % prepare class array
        classes = unique(TrainSet.y);

        % prepare predicted label array
        identity = zeros(1, test_num);

        for i = 1 : test_num

            y = TestSet.X(:, i);

            % calculate sparse code
            xp = l1_ls(TrainSet.X, y, lambda, 1e-3, 1);

            % prepare residual array
            residuals = zeros(1, class_num);

            % calculate residual for each class
            for j = 1 : class_num
                idx = find(TrainSet.y == classes(j));
                residuals(j) = norm(y-TrainSet.X(:,idx)*xp(idx))/sum(xp(idx).*xp(idx));
                %residuals(j) = norm(y - TrainSet.X(:,idx)*xp(idx));
            end

            % calculate the predicted label with minimum residual
            [dis, label] = min(residuals);
            identity(i) = label;

            if verbose
                correct = (label == TestSet.y(1, i));
                fprintf('# RobustPCA-SRC: test:%03d, predict class: %03d --> ground truth :%03d (%d)\n', i, label, TestSet.y(1, i), correct);
            end

        end

        % calculate accuracy
        correct_num = sum(identity == TestSet.y);
        accuracy = correct_num/test_num;
    end
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function test_suite = test_regression_additive1

%   Run specific demo, save values and compare the results to the expected.
%   Works for both xUnit Test Framework package by Steve Eddins and for
%   the built-in Unit Testing Framework (as of Matlab version 2013b).
%
%   See also
%     TEST_ALL, DEMO_REGRESSION_ADDITIVE1
%
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2014 Tuomas Sivula

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  % Check which package to use
  if exist('initTestSuite', 'file')
    initTestSuite_path = which('initTestSuite');
  else
    initTestSuite_path = '';
  end
  if ~isempty(initTestSuite_path) && ...
     exist([initTestSuite_path(1:end-15) 'runtests'], 'file')
    % xUnit package
    initTestSuite;
  else
    % Built-in package
    % Use all functions except the @setup
    tests = localfunctions;
    tests = tests(~cellfun(@(x)strcmp(func2str(x),'setup'),tests));
    test_suite = functiontests(tests);
  end
end

% -------------
%     Tests
% -------------

function testRunDemo(testCase)
  % Run the correspondin demo and save the values. Note this test has to
  % be run at lest once before the other test may succeed.
  rundemo(getName())
end

function testPredictionsFIC(testCase)
  verifyVarsEqual(testCase, getName(), ...
    {'Eft_fic', 'Varft_fic'}, {[], @mean}, ...
    'AbsTol', 0.1)
end

function testPredictionsPIC(testCase)
  verifyVarsEqual(testCase, getName(), ...
    {'Eft_pic', 'Varft_pic'}, ...
    'AbsTol', 0.1)
end

function testPredictionsSparse(testCase)
  verifyVarsEqual(testCase, getName(), ...
    {'Eft_csfic', 'Varft_csfic'}, ...
    'AbsTol', 0.1)
end

% ------------------------
%     Helper functions
% ------------------------

function testCase = setup
  % Helper function to suply empty array into variable testCase as an
  % argument for each test function, if using xUnit package. Not to be
  % used with built-in test framework.
  testCase = [];
end

function name = getName
  % Helperfunction that returns the name of the demo, e.g. 'binomial1'.
  name = mfilename;
  name = name(6:end);
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function MOV = PQavgMOVB (MOVC, Nchan, Nwup)
% Time average MOV precursors

% P. Kabal $Revision: 1.1 $  $Date: 2003/12/07 13:34:46 $

Fs = 48000;
NF = 2048;
Nadv = NF / 2;
Fss = Fs / Nadv;
tdel = 0.5;
tex = 0.050;

% BandwidthRefB, BandwidthTestB
[MOV(0+1), MOV(1+1)] = PQ_avgBW (MOVC.BW);

% Total NMRB, RelDistFramesB
[MOV(2+1), MOV(10+1)] = PQ_avgNMRB (MOVC.NMR);

% WinModDiff1B, AvgModDiff1B, AvgModDiff2B
N500ms = ceil (tdel * Fss);
Ndel = max (0, N500ms - Nwup);
[MOV(3+1), MOV(6+1), MOV(7+1)] = PQ_avgModDiffB (Ndel, MOVC.MDiff);

% RmsNoiseLoudB
N50ms = ceil (tex * Fss);
Nloud = PQloudTest (MOVC.Loud);
Ndel = max (Nloud + N50ms, Ndel);
MOV(8+1) = PQ_avgNLoudB (Ndel, MOVC.NLoud);

% ADBB, MFPDB
[MOV(4+1), MOV(9+1)] = PQ_avgPD (MOVC.PD);

% EHSB
MOV(5+1) = PQ_avgEHS (MOVC.EHS);

%-----------------------------------------
function EHSB = PQ_avgEHS (EHS)

[Nchan, Np] = size (EHS.EHS);

s = 0;
for (j = 0:Nchan-1)
    s = s + PQ_LinPosAvg (EHS.EHS(j+1,:));
end
EHSB = 1000 * s / Nchan;

%-----------------------------------------
function [ADBB, MFPDB] = PQ_avgPD (PD)

global PQopt

c0 = 0.9;
if (isempty (PQopt))
    c1 = 1;
else
    c1 = PQopt.PDfactor;
end

N = length (PD.Pc);
Phc = 0;
Pcmax = 0;
Qsum = 0;
nd = 0;
for (i = 0:N-1)
    Phc = c0 * Phc + (1 - c0) * PD.Pc(i+1);
    Pcmax = max (Pcmax * c1, Phc);

    if (PD.Pc(i+1) > 0.5)
        nd = nd + 1;
        Qsum = Qsum + PD.Qc(i+1);
    end
end

if (nd == 0)
    ADBB = 0;
elseif (Qsum > 0)
    ADBB = log10 (Qsum / nd);
else
    ADBB = -0.5;
end

MFPDB = Pcmax;

%-----------------------------------------
function [TotalNMRB, RelDistFramesB] = PQ_avgNMRB (NMR)

[Nchan, Np] = size (NMR.NMRavg);
Thr = 10^(1.5 / 10);

s = 0;
for (j = 0:Nchan-1)
    s = s + 10 * log10 (PQ_LinAvg (NMR.NMRavg(j+1,:)));
end
TotalNMRB = s / Nchan;

s = 0;
for (j = 0:Nchan-1)
    s = s + PQ_FractThr (Thr, NMR.NMRmax(j+1,:));
end
RelDistFramesB = s / Nchan;

%-----------------------------------------
function [BandwidthRefB, BandwidthTestB] = PQ_avgBW (BW)

[Nchan, Np] = size (BW.BWRef);

sR = 0;
sT = 0;
for (j = 0:Nchan-1)
    sR = sR + PQ_LinPosAvg (BW.BWRef(j+1,:));
    sT = sT + PQ_LinPosAvg (BW.BWTest(j+1,:));
end
BandwidthRefB  = sR / Nchan;
BandwidthTestB = sT / Nchan;

%-----------------------------------------
function [WinModDiff1B, AvgModDiff1B, AvgModDiff2B] = PQ_avgModDiffB (Ndel, MDiff)

NF = 2048;
Nadv = NF / 2;
Fs = 48000;

Fss = Fs / Nadv;
tavg = 0.1;

[Nchan, Np] = size (MDiff.Mt1B);

% Sliding window average - delayed average
L = floor (tavg * Fss);     % 100 ms sliding window length
s = 0;
for (j = 0:Nchan-1)
    s = s + PQ_WinAvg (L, MDiff.Mt1B(j+1,Ndel+1:Np-1+1));
end
WinModDiff1B = s / Nchan;

% Weighted linear average - delayed average
s = 0;
for (j = 0:Nchan-1)
    s = s + PQ_WtAvg (MDiff.Mt1B(j+1,Ndel+1:Np-1+1), MDiff.Wt(j+1,Ndel+1:Np-1+1));
end
AvgModDiff1B = s / Nchan;

% Weighted linear average - delayed average
s = 0;
for (j = 0:Nchan-1)
    s = s + PQ_WtAvg (MDiff.Mt2B(j+1,Ndel+1:Np-1+1), MDiff.Wt(j+1,Ndel+1:Np-1+1));
end
AvgModDiff2B = s / Nchan;

%-----------------------------------------
function RmsNoiseLoudB = PQ_avgNLoudB (Ndel, NLoud)

[Nchan, Np] = size (NLoud.NL);

% RMS average - delayed average and loudness threshold
s = 0;
for (j = 0:Nchan-1)
    s = s + PQ_RMSAvg (NLoud.NL(j+1,Ndel+1:Np-1+1));
end
RmsNoiseLoudB = s / Nchan;

%-----------------------------------
% Average values values, omitting values which are negative
function s = PQ_LinPosAvg (x)

N = length(x);

Nv = 0;
s = 0;
for (i = 0:N-1)
    if (x(i+1) >= 0)
        s = s + x(i+1);
        Nv = Nv + 1;
    end
end

if (Nv > 0)
    s = s / Nv;
end

%----------
% Fraction of values above a threshold
function Fd = PQ_FractThr (Thr, x)

N = length (x);

Nv = 0;
for (i = 0:N-1)
    if (x(i+1) > Thr)
        Nv = Nv + 1;
    end
end

if (N > 0)
    Fd = Nv / N;
else
    Fd = 0;
end

%-----------
% Sliding window (L samples) average
function s = PQ_WinAvg (L, x)

N = length (x);

s = 0;
for (i = L-1:N-1)
    t = 0;
    for (m = 0:L-1)
        t = t + sqrt (x(i-m+1));
    end
    s = s + (t / L)^4;
end

if (N >= L)
    s = sqrt (s / (N - L + 1));
end

%----------
% Weighted average
function s = PQ_WtAvg (x, W)

N = length (x);

s = 0;
sW = 0;
for (i = 0:N-1)
    s = s + W(i+1) * x(i+1);
    sW = sW + W(i+1);
end

if (N > 0)
    s = s / sW;
end

%----------
% Linear average
function LinAvg = PQ_LinAvg (x)

N = length (x);
s = 0;
for (i = 0:N-1)
    s = s + x(i+1);
end

LinAvg = s / N;

%----------
% Square root of average of squared values
function RMSAvg = PQ_RMSAvg (x)

N = length (x);
s = 0;
for (i = 0:N-1)
    s = s + x(i+1)^2;
end

if (N > 0)
    RMSAvg = sqrt(s / N);
else
    RMSAvg = 0;
end

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
function [Epb Eav Etl hst] = entropia_f(I,Iclass)

% Entropy (E) of intensity image

% 06/08/2011    - Version 4.0F

% Formula: E = -sum(p.*log2(p));

% Author:           Aristidis D. Vaiopoulos
% Acknowledgement:  Mathworks MATLAB entropy function

% If image class is not forced by user, detect and use native class
if nargin == 1;
    Iclass = class(I);
end
% Construct the bin values and convert image according to Iclass
switch Iclass
   case {'logical'}
%       cs = 1;
      cl = 2^1;
      I = logical(I);
      bv = [0 1];
   case {'uint8'}
%       cs = 2;
      cl = 2^8;
      I = im2uint8(I);
      bv = 0:(2^8)-1;
   case {'int16'}
%       cs = 3;
      cl = 2^16;
      I = im2int16(I);
      bv = -( (2^16)/2 ):( (2^16)/2-1 );
   case {'uint16'}
%       cs = 4;
      cl = 2^16;
      I = im2uint16(I);
      bv = 0:(2^16)-1;
   case {'single','double'}
%       cs = 5;
      cl = 2^16;
      I = im2uint16(I);
      bv = 0:(2^16)-1;
    otherwise
%       cs = -1;
      disp('-Unsupported data class.')
    return;

end
% Transpose bin values
bv = bv';

% Find the number of bands
bands = size(I);
if length(bands) == 3
    bands = bands(1,3);
else
    bands = 1;
end
% Calculate histogram counts per band
p = zeros(cl,bands);
for b = 1:bands
    p(:,b) = imhist(I(:,:,b),cl);
end

% Give histogram
if nargout == 4
    hst = [bv p];
end
% normalize p so that sum(p) is one.
p = p ./ (numel(I)/bands);
% normalize p for whole image.
pt = sum(p,2)/bands;
% logarithmization
lp = log2(p);
lpt= log2(pt);
% nullify -Inf values due to logarithmization of zeros in p
lp(isinf(lp)) = 0;
lpt(isinf(lpt)) = 0;
% Entropy per band
Epb = -sum(p.*lp,1)';
% Average Entropy
Eav = mean(Epb);
% Total Entropy (whole image)
Etl = -sum(pt.*lpt);

end


text:
% pop_firma() - Filter data using moving average FIR filter
%
% Usage:
%   >> [EEG, com] = pop_firma(EEG); % pop-up window mode
%   >> [EEG, com] = pop_firma(EEG, 'forder', order);
%
% Inputs:
%   EEG       - EEGLAB EEG structure
%   'forder'  - scalar filter order. Mandatory even
%
% Outputs:
%   EEG       - filtered EEGLAB EEG structure
%   com       - history string
%
% Author: Andreas Widmann, University of Leipzig, 2005
%
% See also:
%   firfilt, plotfresp

%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) 2005 Andreas Widmann, University of Leipzig, widmann@uni-leipzig.de
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [EEG, com] = pop_firma(EEG, varargin)

    com = '';
    if nargin < 1
        help pop_firma;
        return;
    end
    if isempty(EEG.data)
        error('Cannot process empty dataset');
    end

    if nargin < 2
        drawnow;
        uigeom = {[1 1 1] [1] [1 1 1]};
        uilist = {{'style' 'text' 'string' 'Filter order (mandatory even):'} ...
                  {'style' 'edit' 'string' '' 'tag' 'forderedit'} {} ...
                  {} ...
                  {} {} {'Style' 'pushbutton' 'string' 'Plot filter responses' 'callback' {@complot, EEG.srate}}};
        result = inputgui(uigeom, uilist, 'pophelp(''pop_firma'')', 'Filter the data -- pop_firma()');
        if length(result) == 0, return; end

        if ~isempty(result{1})
            args = [{'forder'} {str2num(result{1})}];
        else
            error('Not enough input arguments');
        end
    else
        args = varargin;
    end

    % Convert args to structure
    args = struct(args{:});

    % Filter coefficients
    b = ones(1, args.forder + 1) / (args.forder + 1);

    % Filter
    disp('pop_firma() - filtering the data');
    EEG = firfilt(EEG, b);

    % History string
    com = sprintf('%s = pop_firma(%s', inputname(1), inputname(1));
    for c = fieldnames(args)'
        if ischar(args.(c{:}))
            com = [com sprintf(', ''%s'', ''%s''', c{:}, args.(c{:}))];
        else
            com = [com sprintf(', ''%s'', %s', c{:}, mat2str(args.(c{:})))];
        end
    end
    com = [com ');'];

% Callback plot filter properties
function complot(obj, evt, srate)
    args.forder = str2num(get(findobj(gcbf, 'tag', 'forderedit'), 'string'));
    if isempty(args.forder)
        error('Not enough input arguments');
    end
    b = ones(1, args.forder + 1) / (args.forder + 1);
    H = findobj('tag', 'filter responses', 'type', 'figure');
    if ~isempty(H)
        figure(H);
    else
        H = figure;
        set(H, 'color', [.93 .96 1], 'tag', 'filter responses');
    end
    plotfresp(b, 1, [], srate);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


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

%%
clear; close all; clc;

%% Syntax
% |[maxVal,maxInd]=spmax(A,B,dim,nanflag,logicRelevant,nanOut);|

%% Description
% This function is like the max function but is designed for sparse arrays.
% In particular it allows one to "ignore zeros" in the determination of the
% maxima.

%% Examples
%

%%
% Create example matrix
i=[2 1 1 2  2 3  3 4 4 5  5 5 6 6 7 8];
j=[1 1 2 3  4 5  6 7 8 9 10 11 12 13 13 13];
s=[-1 3 1 2 -1 1 -2 5 5 -1 0 2 3 10 11 NaN];
siz=max([i(:);j(:)]+1)*ones(1,2);
A=sparse(i,j,s,siz(1),siz(2),numel(s));
A=A+A';

full(A) % View matrix

L=sparse(i,j,1,siz(1),siz(2),numel(s));
logicRelevant=(L+L')>0;

%%
% Compute maxima allong a certain direction (while omit nan is default)

amaxRows=spmax(A,[],1);
full(amaxRows)

amaxColumns=spmax(A,[],2);
full(amaxColumns)

%%
% Including nans

amaxRows=spmax(A,[],1,'includenan');
full(amaxRows)

amaxColumns=spmax(A,[],2,'includenan');
full(amaxColumns)

%%
% Computing maxima across all desired relevant entries (including
% "relevant/real zeros")

amaxRows=spmax(A,[],1,'omitnan',logicRelevant);
full(amaxRows)

amaxColumns=spmax(A,[],2,'omitnan',logicRelevant);
full(amaxColumns)

%%
% Computin maxima across all desired relevant entries and output NaN where
% the sparse array only contains "non-relevant or non-real" zeros.

nanOut=1;

amaxRows=spmax(A,[],1,'omitnan',logicRelevant,nanOut);
full(amaxRows)

amaxColumns=spmax(A,[],2,'omitnan',logicRelevant,nanOut);
full(amaxColumns)

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


text:
% Sample codes for the task-related component analysis (TRCA)-based steady
% -state visual evoked potential (SSVEP) detection method [1]. The filter
% bank analysis [2] can also be combined to the TRCA-based algorithm.
%
% Dataset (sample.mat):
%   A 40-target SSVEP dataset recorded from a single subject. The stimuli
%   were generated by the joint frequency-phase modulation (JFPM) [3]
%     - Stimulus frequencies    : 8.0 - 15.8 Hz with an interval of 0.2 Hz
%     - Stimulus phases         : 0pi, 0.5pi, 1.0pi, and 1.5pi
%     - # of channels           : 9 (1: Pz, 2: PO5,3:  PO3, 4: POz, 5: PO4,
%                                    6: PO6, 7: O1, 8: Oz, and 9: O2)
%     - # of recording blocks   : 6
%     - Data length of epochs   : 5 [seconds]
%     - Sampling rate           : 250 [Hz]
%
% See also:
%   train_trca.m
%   test_trca.m
%   filterbank.m
%   itr.m
%
% Reference:
%   [1] M. Nakanishi, Y. Wang, X. Chen, Y.-T. Wang, X. Gao, and T.-P. Jung,
%       "Enhancing detection of SSVEPs for a high-speed brain speller using
%        task-related component analysis",
%       IEEE Trans. Biomed. Eng, 65(1): 104-112, 2018.
%   [2] X. Chen, Y. Wang, S. Gao, T. -P. Jung and X. Gao,
%       "Filter bank canonical correlation analysis for implementing a
%        high-speed SSVEP-based brain-computer interface",
%       J. Neural Eng., 12: 046008, 2015.
%   [3] X. Chen, Y. Wang, M. Nakanishi, X. Gao, T. -P. Jung, S. Gao,
%       "High-speed spelling with a noninvasive brain-computer interface",
%       Proc. Int. Natl. Acad. Sci. U. S. A, 112(44): E6058-6067, 2015.
%
% Masaki Nakanishi, 22-Dec-2017
% Swartz Center for Computational Neuroscience, Institute for Neural
% Computation, University of California San Diego
% E-mail: masaki@sccn.ucsd.edu

%% Clear workspace
clear all
close all
clc
help tutorial_trca

%% Set paths

addpath('../src');

%% Parameter for analysis (Modify according to your analysis)

% Filename
filename = '../data/sample.mat';

% Data length for target identification [s]
len_gaze_s = 0.5;

% Visual latency being considered in the analysis [s]
len_delay_s = 0.13;

% The number of sub-bands in filter bank analysis
num_fbs = 5;

% 1 -> The ensemble TRCA-based method, 0 -> The TRCA-based method
is_ensemble = 0;

% 100*(1-alpha_ci): confidence intervals
alpha_ci = 0.05;

%% Fixed parameter (Modify according to the experimental setting)

% Sampling rate [Hz]
fs = 250;

% Duration for gaze shifting [s]
len_shift_s = 0.5;

% List of stimulus frequencies
list_freqs = [8:1:15 8.2:1:15.2 8.4:1:15.4 8.6:1:15.6 8.8:1:15.8];

% The number of stimuli
num_targs = length(list_freqs);

% Labels of data
labels = [1:1:num_targs];

%% Preparing useful variables (DONT'T need to modify)

% Data length [samples]
len_gaze_smpl = round(len_gaze_s*fs);

% Visual latency [samples]
len_delay_smpl = round(len_delay_s*fs);

% Selection time [s]
len_sel_s = len_gaze_s + len_shift_s;

% Confidence interval
ci = 100*(1-alpha_ci);

%% Performing the TRCA-based SSVEP detection algorithm

fprintf('Results of the ensemble TRCA-based method.\n');

% Preparing data
load(filename);
[~, num_chans, ~, num_blocks] = size(eeg);
segment_data = len_delay_smpl+1:len_delay_smpl+len_gaze_smpl;
eeg = eeg(:, :, segment_data, :);

% Estimate classification performance
for loocv_i = 1:1:num_blocks

    % Training stage
    traindata = eeg;
    traindata(:, :, :, loocv_i) = [];
    model = train_trca(traindata, fs, num_fbs);

    % Test stage
    testdata = squeeze(eeg(:, :, :, loocv_i));
    estimated = test_trca(testdata, model, is_ensemble);

    % Evaluation
    is_correct = (estimated==labels);
    accs(loocv_i) = mean(is_correct)*100;
    itrs(loocv_i) = itr(num_targs, mean(is_correct), len_sel_s);
    fprintf('Trial %d: Accuracy = %2.2f%%, ITR = %2.2f bpm\n',...
        loocv_i, accs(loocv_i), itrs(loocv_i));

end % loocv_i

% Summarize
[mu, ~, muci, ~] = normfit(accs, alpha_ci);
fprintf('Mean accuracy = %2.2f %% (%2d%% CI: %2.2f - %2.2f %%)\n',...
    mu, ci, muci(1), muci(2));

[mu, ~, muci, ~] = normfit(itrs, alpha_ci);
fprintf('Mean ITR = %2.2f bpm (%2d%% CI: %2.2f - %2.2f bpm)\n\n',...
    mu, ci, muci(1), muci(2));

author: mnakanishi

repo: TRCA-SSVEP


text:
function X_rec = recoverData(Z, U, K)
%RECOVERDATA Recovers an approximation of the original data when using the
%projected data
%   X_rec = RECOVERDATA(Z, U, K) recovers an approximation the
%   original data that has been reduced to K dimensions. It returns the
%   approximate reconstruction in X_rec.
%

% You need to return the following variables correctly.
X_rec = zeros(size(Z, 1), size(U, 1));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the approximation of the data by projecting back
%               onto the original space using the top K eigenvectors in U.
%
%               For the i-th example Z(i,:), the (approximate)
%               recovered data for dimension j is given as follows:
%                    v = Z(i, :)';
%                    recovered_j = v' * U(j, 1:K)';
%
%               Notice that U(j, 1:K) is a row vector.
%

for i = 1 : size(Z, 1),
    for j = 1 : size(U,1),
        v = Z(i, :)';
        rec = v' * U(j, 1 : K)';
        X_rec(i, j) = rec;
    end
end

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

end

author: UtkarshPathrabe


text:
function coeffs3D = vals2coeffs(F)
%VALS2COEFFS   Convert tensor of values to tensor of coefficients.
%
% See also CHEBFUN3/COEFFS2VALS.

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

[m, n, p] = size(F);

%% Step 1:
% Mode-1 unfolding of F to get a matrix of size m x n*p:
F1 = chebfun3.unfold(F, 1);

% Apply 1D "vals2coeffs" in the X direction i.e., the 1st dimension of F1:
F1 = chebtech2.vals2coeffs(F1);

% Tensorize F1 back to its original m x n x p size:
F1 = chebfun3.fold(F1, [m, n, p], 1, [2 3]);
% This is simply F1 = reshape(V1, m, n, p);

%% Step 2:
% Mode-2 unfolding of (the tensorized) F1 to get a matrix of size m x n*p:
F2 = chebfun3.unfold(F1, 2);

% Apply 1D "vals2coeffs" in the Y direction, i.e. to the 1st direction of
% F2:
F2 = chebtech2.vals2coeffs(F2);

% Tensorize F2 back to its original m x n x p size:
F2 = chebfun3.fold(F2, [m, n, p], 2, [1 3]);

%% Step 3:
% Mode-3 unfolding of (the tensorized) F2 to get a matrix of size p x m*n:
F3 = chebfun3.unfold(F2, 3);

% Now, vals2coeffs is applied in the Z direction.
F3 = chebtech2.vals2coeffs(F3);

% Reshape F3 back to the original m x n x p size:
coeffs3D = chebfun3.fold(F3, [m, n, p], 3, [1 2]);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun