text:
function [no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method,isovalues)
%
% [no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method,isovalues)
%
% converting a 3D volumetric image to surfaces
%
% author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
%
% input:
%	 img: a volumetric binary image; if img is empty, vol2surf will
%	      return user defined surfaces via opt.surf if it exists
%	 ix,iy,iz: subvolume selection indices in x,y,z directions
%	 opt: function parameters
%	   if method is 'cgalsurf' or 'cgalpoly':
%	     opt=a float number>1: max radius of the Delaunay sphere(element size)
%	     opt.radbound: same as above, max radius of the Delaunay sphere
%	     opt.distbound: maximum deviation from the specified isosurfaces
%	     opt(1,2,...).radbound: same as above, for each levelset
%	   if method is 'simplify':
%	     opt=a float number<1: compression rate for surf. simplification
%	     opt.keepratio=a float less than 1: same as above, same for all surf.
%	     opt(1,2,..).keepratio: setting compression rate for each levelset
%	   opt(1,2,..).maxsurf: 1 - only use the largest disjointed surface
%				0 - use all surfaces for that levelset
%          opt(1,2,..).side: - 'upper': threshold at upper interface
%                              'lower': threshold at lower interface
%	   opt(1,2,..).maxnode: - the maximum number of surface node per levelset
%	   opt(1,2,..).holes: user specified holes interior pt list
%	   opt(1,2,..).regions: user specified regions interior pt list
%	   opt(1,2,..).surf.{node,elem}: add additional surfaces
%	   opt(1,2,..).{A,B}: linear transformation for each surface
%	   opt.autoregion: if set to 1, vol2surf will try to determine
%              the interior points for each closed surface automatically
%	 dofix: 1: perform mesh validation&repair, 0: skip repairing
%	 method: - if method is 'simplify', iso2mesh will first call
%		   binsurface to generate a voxel-based surface mesh and then
%		   use meshresample/meshcheckrepair to create a coarser mesh;
%		 - if method is 'cgalsurf', iso2mesh will call the surface
%		   extraction program from CGAL to make surface mesh
%		 - if method is not specified, 'cgalsurf' is assumed by default
%	 isovalues: a list of isovalues where the levelset is defined
%
% output:
%	 no:  list of nodes on the resulting suface mesh, 3 columns for x,y,z
%	 el:  list of trianglular elements on the surface, [n1,n2,n3,region_id]
%	 regions: list of interior points for all sub-region, [x,y,z]
%	 holes:   list of interior points for all holes, [x,y,z]
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

fprintf(1,'extracting surfaces from a volume ...\n');

el=[];
no=[];

if(isstruct(opt) & isfield(opt,'holes'))
    holes=opt.holes;
else
    holes=[];
end
if(isstruct(opt) & isfield(opt,'regions'))
    regions=opt.regions;
else
    regions=[];
end
maxlevel=0;

if(~isempty(img))

    img=img(ix,iy,iz);
    dim=size(img);
    newdim=dim+[2 2 2];
    newimg=zeros(newdim);
    newimg(2:end-1,2:end-1,2:end-1)=img;

    if(nargin<8)
        maxlevel=max(newimg(:));
        isovalues=1:maxlevel;
    else
        isovalues=unique(sort(isovalues));
        maxlevel=length(isovalues);
    end

    for i=1:maxlevel
      if(i<maxlevel)
          levelmask=int8(newimg>=isovalues(i) & newimg<isovalues(i+1));
      else
          levelmask=int8(newimg>=isovalues(i));
      end
      [levelno,levelel]=binsurface(levelmask);
      if(~isempty(levelel))
          if(isstruct(opt) & isfield(opt,'autoregion'))
              if(opt.autoregion)
                  seeds=surfseeds(levelno,levelel);
              else
                  seeds=surfinterior(levelno,levelel);
              end
          else
              seeds=surfinterior(levelno,levelel);
          end
          if(~isempty(seeds))
              disp([sprintf('region %d centroid :',i) sprintf('\t%f %f %f\n', seeds')]);
              if(~isempty(regions))
                  regions(end+1:end+size(seeds,1),:)=seeds;
              else
                  regions=seeds;
              end
          end
      end
    end

    for i=1:maxlevel
        fprintf(1,'processing threshold level %d...\n',i);

        if(nargin>=7 & strcmp(method,'simplify'))

          [v0,f0]=binsurface(newimg>=isovalues(i)); % not sure if binsurface works for multi-value arrays
          % with binsurface, I think the following line is not needed anymore
          %  v0(:,[1 2])=v0(:,[2 1]); % isosurface(V,th) assumes x/y transposed
          if(dofix)  [v0,f0]=meshcheckrepair(v0,f0);  end

          if(isstruct(opt) & length(opt)==maxlevel) keepratio=opt(i).keepratio;
          elseif (isstruct(opt) & length(opt)==1) keepratio=opt.keepratio;
          else keepratio=opt;  end;

          % first, resample the surface mesh with cgal
          fprintf(1,'resampling surface mesh for level %d...\n',i);
          [v0,f0]=meshresample(v0,f0,keepratio);

          % iso2mesh is not stable for meshing small islands,remove them (max 3x3x3 voxels)
          f0=removeisolatedsurf(v0,f0,3);

          if(dofix) [v0,f0]=meshcheckrepair(v0,f0); end

        elseif(nargin<7 | strcmp(method,'cgalsurf') | strcmp(method,'cgalpoly'))
          if(isstruct(opt) & length(opt)==maxlevel) radbound=opt(i).radbound;
          elseif (isstruct(opt) & length(opt)==1) radbound=opt.radbound;
          else radbound=opt;  end;

          maxsurfnode=40000;  % maximum node numbers for each level
          if(isstruct(opt) & length(opt)==maxlevel)
              if(isfield(opt(i),'maxnode')) maxsurfnode=opt(i).maxnode; end
          elseif (isstruct(opt) & length(opt)==1 )
              if(isfield(opt(1),'maxnode'))
                 maxsurfnode=opt.maxnode;
              end
          end

	  distbound=radbound;
          if(isstruct(opt) & length(opt)==maxlevel)
              if(isfield(opt(i),'distbound')) distbound=opt(i).distbound; end
          elseif (isstruct(opt) & length(opt)==1 )
              if(isfield(opt(1),'distbound')) distbound=opt.distbound; end
          end
	  surfside='';
          if(isstruct(opt) & length(opt)==maxlevel)
              if(isfield(opt(i),'side')) surfside=opt(i).side; end
          elseif (isstruct(opt) & length(opt)==1)
              if(isfield(opt(1),'side')) surfside=opt(1).side; end
          end
	  if(~isempty(surfside))
	     newimg0=newimg;
	  end
          if(strcmp(surfside,'upper'))
	      newimg(find(newimg<=isovalues(i)-1e-9))=isovalues(i)-1e-9;
	  elseif(strcmp(surfside,'lower'))
	      newimg(find(newimg>=isovalues(i)+1e-9))=isovalues(i)+1e-9;
	  end
          perturb=1e-4*abs(max(isovalues));
          if(all(newimg>isovalues(i)-perturb)) perturb=-perturb;  end
          [v0,f0]=vol2restrictedtri(newimg,isovalues(i)-perturb,regions(i,:),...
                     sum(newdim.*newdim)*2,30,radbound,distbound,maxsurfnode);

	  if(~isempty(surfside))
            newimg=newimg0;
	    clear newimg0;
	  end
        else
            error('method can only be one of "cgalsurf", "cgalpoly" or "simplify".');
        end

        % if use defines maxsurf=1, take only the largest closed surface
        if(isstruct(opt))
            if( (isfield(opt,'maxsurf') && length(opt)==1 && opt.maxsurf==1) | ...
                    (length(opt)==maxlevel && isfield(opt(i),'maxsurf') && opt(i).maxsurf==1))
                    f0=maxsurf(finddisconnsurf(f0));
            end
        end

        % if a transformation matrix/offset vector supplied, apply them

        if(isstruct(opt) & length(opt)==maxlevel)
          if(isfield(opt(i),'A') & isfield(opt(i),'B'))
            v0=(opt(i).A*v0'+repmat(opt(i).B(:),1,size(v0,1)))';
          end
        elseif (isstruct(opt) & length(opt)==1)
          if(isfield(opt,'A') & isfield(opt,'B'))
            v0=(opt.A*v0'+repmat(opt.B(:),1,size(v0,1)))';
          end
        end

        % if user specified holelist and regionlist, append them
        if(isstruct(opt)  & length(opt)==maxlevel)
        if(isfield(opt(i),'hole'))
                holes=[holes;opt(i).hole]
        end
        if(isfield(opt(i),'region'))
                regions=[regions;opt(i).region]
        end
        end

        if(i==0)
        el=[f0 (i+1)*ones(size(f0,1),1)];
        no=v0;
        else
        el=[el;f0+length(no) (i+1)*ones(size(f0,1),1)];
        no=[no;v0];
        end
    end

    %some final fix and scaling
    no(:,1:3)=no(:,1:3)-1; % because we padded the image with a 1 voxel thick null layer in newimg

    no(:,1)=no(:,1)*(max(ix)-min(ix)+1)/dim(1)+(min(ix)-1);
    no(:,2)=no(:,2)*(max(iy)-min(iy)+1)/dim(2)+(min(iy)-1);
    no(:,3)=no(:,3)*(max(iz)-min(iz)+1)/dim(3)+(min(iz)-1);

end  % if not isempty(img)

if(isstruct(opt) & isfield(opt,'surf'))
   for i=1:length(opt.surf)
	opt.surf(i).elem(:,4)=maxlevel+i;
        el=[el;opt.surf(i).elem+length(no)];
        no=[no;opt.surf(i).node];
   end
end

fprintf(1,'surface mesh generation is complete\n');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%
% Bayesian Optimization of Combinatorial Structures
%
% Copyright (C) 2018 R. Baptista & M. Poloczek
%
% BOCS 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.
%
% BOCS 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 BOCS.  If not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) 2018 MIT & University of Arizona
% Authors: Ricardo Baptista & Matthias Poloczek
% E-mails: rsb@mit.edu & poloczek@email.arizona.edu
%

% Script runs discrete optimization algorithms for the
% binary quadratic programming problem. The results are
% compared for different values of the \lambda tuning
% parameter.

clear; close all; clc
addpath(genpath('../algorithms'))
addpath(genpath('../stat_model'))
addpath(genpath('../test_problems/Quadratic'))
addpath(genpath('../tools'))

%% Setup parameters

% Setup fixed parameters
n_vars  = 10;
n_proc  = 1;
test_name = 'quad';

% Number of runs and optimization iterations
n_func     = 50;
n_runs     = 10;
n_init     = 20;
evalBudget = 120;

% problem parameters (Monte Carlo samples)
alpha_vect = logspace(0,2,3);

% Variance prior parameters (Inverse Gamma)
aPr    = 2;
bPr    = 1;

% Regularization parameters
lambda_vals = [0, 1e-4, 1e-2, 1];
lambda_str  = {'0', '1em4', '1em2', '1'};

% Set additive regularization function
reg_term = @(x) sum(x,2);

%% Generate Test Cases

% setup objective functions
inputs_all = cell(n_func, length(alpha_vect), n_runs);

for a=1:length(alpha_vect)

    % compute decay function
    K = @(s,t) exp(-1*(s-t)^2/alpha_vect(a));
    decay = zeros(n_vars,n_vars);

    for i=1:n_vars
        for j=1:n_vars
            decay(i,j) = K(i,j);
        end
    end

    for t1=1:n_func

        % Generate random quadratic model
        % apply exponential decay to Q
        Q = randn(n_vars, n_vars);
        Qa = Q.*decay;

        for t2=1:n_runs

            % Set inputs struct for each problem
            inputs_all{t1,a,t2} = struct;
            inputs_all{t1,a,t2}.n_vars      = n_vars;

            % Save other definitions
            inputs_all{t1,a,t2}.evalBudget  = evalBudget;
            inputs_all{t1,a,t2}.n_runs      = n_runs;
            inputs_all{t1,a,t2}.n_init      = n_init;
            inputs_all{t1,a,t2}.lambda_vals = lambda_vals;

            % Set priors for estimator
            inputs_all{t1,a,t2}.aPr         = aPr;
            inputs_all{t1,a,t2}.bPr         = bPr;

            % Setup true objective function
            inputs_all{t1,a,t2}.model = @(x) diag(x*Qa*x');

            % Save regularization term
            inputs_all{t1,a,t2}.reg_term = @(x) reg_term(x);

            % Generate initial samples for statistical models
            inputs_all{t1,a,t2}.x_vals = sample_models(n_init, n_vars);
            inputs_all{t1,a,t2}.y_vals = inputs_all{t1,a,t2}.model(inputs_all{t1,a,t2}.x_vals);

        end
    end
end

inputs_all = reshape(inputs_all, n_func*length(alpha_vect)*n_runs, 1);

% Make folder
mkdir(['../results/' test_name])

% Save test cases
save(['../results/' test_name '/all_tests'])

% Run cases
run_cases(inputs_all, lambda_vals, test_name, n_proc);

% -- END OF FILE --

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
function paths = searchFiles(submodules, options)
%
% Syntax:
%    paths = searchFiles(submodules, options)
%
% Usage:
%    paths = searchFiles()
%    paths = searchFiles(submodules)
%    paths = searchFiles(submodules, options)
%
% Examples:
%    paths = searchFiles()
%    paths = searchFiles({}, 'update')
%
%
global paths  %#ok<REDEFGI>

if ~exist('submodules','var') || isempty(submodules)
    submodules = parseGitSubmodulesFile();
end
if ~exist('options','var')
    options = '';
end
paths = {};

submodulenames = {};
submodulepaths = {};
urls = {};
kk = 1;
for ii = 1:size(submodules,1)
    urls{ii}            = submodules{ii,1};
    submodulepaths{ii}  = submodules{ii,3};
    if optionExists_startup(options,'update')
        removeFolderContents(submodulepaths{ii});
    elseif ~isemptyFolder(submodulepaths{ii})
        continue;
    end
    [~, submodulenames{kk}] = fileparts(urls{ii});
    kk = kk+1;
end
if isempty(submodulenames)
    return;
end

hf = figure('numbertitle', 'off', 'menubar','none', 'toolbar','none', 'name','Files Paths');
p = get(hf, 'position');
set(hf, 'position',[p(1), p(2), 1.5*p(3), p(4)]);

xpos  = .01;
ypos  =  .70;
xsize =  .20;
ysize =  .07;
xgapsize = .02;
ygapsize = .03;

xpos0 = 12*xpos;
ypos0 = ypos+.1*ypos;
xsize0 = 3.5*xsize;
ysize0 = 2*ysize;
fprintf('Display instructions at [%1.f, %1.f, %1.f, %1.f]\n', [xpos0, ypos0, xsize0, ysize0]);
uicontrol('parent',hf, 'style','text', 'string','Provide paths of downloaded submodules', 'units','normalized', 'position',[xpos0, ypos0, xsize0, ysize0], ...
    'fontsize',11, 'fontweight','bold', 'horizontalalignment','center');

for ii = 1:length(submodulenames)
    fprintf('%d. File %s\n', ii, submodulenames{ii});

    yposi = ypos - (ii-1)*(ysize+ygapsize);
    if yposi<0
        break;
    end

    xpos1 = .7*xpos;
    xsize1 = xsize;
    ysize1 = ysize;
    k = .98;
    fprintf('Display label at [%1.f, %1.f, %1.f, %1.f]\n', [xpos1, k*yposi, xsize1, ysize1]);
    htxt(ii)  = uicontrol('parent',hf, 'style','text',       'string',[submodulenames{ii}, ' :'], 'units','normalized', 'position',[xpos1, k*yposi, xsize1, ysize1], ...
        'fontsize',10, 'fontweight','bold', 'horizontalalignment','right');

    xpos2 = xpos1+xsize1+xgapsize/2;
    xsize2 = 2.5*xsize;
    ysize2 = ysize;
    fprintf('Display edit box at [%1.f, %1.f, %1.f, %1.f]\n', [xpos2, yposi, xsize2, ysize2]);
    hedit(ii) = uicontrol('parent',hf, 'style','edit',       'string','',        'units','normalized', 'position',[xpos2, yposi, xsize2, ysize2], ...
        'fontsize',9, 'fontweight','normal', 'horizontalalignment','left', 'tag',sprintf('edit%sPath',submodulenames{ii}));

    xpos3 = xpos2+xsize2+xgapsize;
    xsize3 = .8*xsize;
    ysize3 = ysize;
    fprintf('Display browse pushbutton at [%1.f, %1.f, %1.f, %1.f]\n', [xpos3, yposi, xsize3, ysize3]);
    hbttn(ii) = uicontrol('parent',hf, 'style','pushbutton', 'string','Browse',  'units','normalized', 'position',[xpos3, yposi, xsize3, ysize3], ...
        'fontsize',10, 'fontweight','bold','callback',{@browseBttnCallback, hedit(ii), submodulenames{ii}});

    fprintf('\n');
end

yposi = ypos - 2*ii*(ysize+ygapsize);
if yposi<0
    return;
end
xpos4  = .35;
xsize4 = 1.2*xsize;
ysize4 = 1.5*ysize;
fprintf('Display SUBMIT button [%1.f, %1.f, %1.f, %1.f]\n', [xpos4, yposi, xsize4, ysize4]);
uicontrol('parent',hf, 'style','pushbutton', 'string','SUBMIT', 'units','normalized', 'position',[xpos4, yposi, xsize4, ysize4], ...
    'fontsize',11, 'fontweight','bold', 'callback',{@submitBttnCallback, hedit, submodules});

waitForGui_startup(hf);

% ------------------------------------------------------
function pname = browseBttnCallback(~, ~, hEdit, submodulename)
pname = '';
[filename, pathname] = uigetfile('*.zip', sprintf('Select file for submodule "%s"', submodulename));
if filename==0
    return
end
pause(.01)

pname = filesepStandard_startup([pathname, '/', filename]);
if isempty(strfind(pname, submodulename)) %#ok<*STREMP>
    msgbox(sprintf('Selected file  "%s"  does not match the submodule  "%s". Please download and provide the path name for submodule  "%s"...', ...
        filename, submodulename, submodulename))
    return;
end
set(hEdit, 'string',pname)

% ------------------------------------------------------
function submitBttnCallback(hObject, ~, hEdits, submodules)
global paths

hf = get(hObject,'parent');

for ii = 1:length(hEdits)
    filepath = get(hEdits(ii), 'string');
    if ~ispathvalid_startup(filepath)
        continue;
    end
    [p,f,e] = fileparts(filepath);
    filenameUnziped = filesepStandard_startup([p,'/',f], 'nameonly:dir');

    if strcmp(e,'.zip')
        if ~ispathvalid_startup(filepath,'file')
            continue
        end
        try
            if ispathvalid_startup(filenameUnziped,'dir')
                rmdir(filenameUnziped,'s')
            end
            unzip(filepath, fileparts(filenameUnziped(1:end-1)));
        catch
        end
    end

    if ~ispathvalid_startup(filenameUnziped,'dir')
        continue
    end

    % Copy downloaded and unzipped submodule folder contents to
    % corresponding submodule folder in the parent repo
    for jj = 1:size(submodules,1)
        submodulepath = submodules{jj,3};
        [~, submodulename] = fileparts(submodules{jj,1});
        if isempty(strfind(filenameUnziped, submodulename))
            continue
        end
        fprintf('Copying %s to %s\n', [filenameUnziped, '/*'], submodulepath);
        copyFolderContents(filenameUnziped, submodulepath);
        rmdir(filenameUnziped,'s');
        paths{ii,1} = filenameUnziped;
        break;
    end
end

delete(hf)

author: BUNPC


text:
% Extract instantaneous phase from a signal
%
% This function extracts instantaneous phase of a signal based on Hilbert transform.
%
%  USAGE
%   phase_deg = general.extractPhase(eeg)
%   eeg             Vector that contains a signal
%   phase_deg       Vector that contains extracted phase. Vector values are in degrees.
%
function [phase_deg] = extractPhase(eeg)
%     phase_rad = phase(hilbert(eeg));
    phase_rad = angle(hilbert(eeg)); % faster than phase
    phase_deg = rad2deg(mod(phase_rad, 2*pi));
end

author: kavli-ntnu


text:
function opts = updateMovie(S, dt, p, opts, t, v, compGrid, plotGrid)
%UPDATEMOVIE   Update the movie when solving a PDE specified by a SPINOPSPHERE.

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

% Set-up:
nVars = S.numVars;
Clim = opts{1};
dataToPlot = opts{3};
ll = compGrid{1};
tt = compGrid{2};
N = 2*(size(ll, 1) - 1);
lll = plotGrid{1};
ttt = plotGrid{2};
Nplot = 2*(size(lll, 1) - 1);

for k = 1:nVars

    % Extract each variable:
    idx = (k-1)*N + 1;
    vv = dataToPlot(v(idx:idx+N-1,:));
    vv = [vv, vv(:,1)]; %#ok<*AGROW> add repeated values (periodic endpoints)
    vv = [vv; vv(1,:)];
    vv = vv([N/2+1:N 1], :); % extract values that correspond to theta in [0,pi]

    % Change axes if necessary:
    if ( nargout == 1 )
        minvnew = min(vv(:));
        maxvnew = max(vv(:));
        if ( maxvnew > Clim(2*(k-1) + 2) )
            vscalenew = max(abs(minvnew), maxvnew);
            Clim(2*(k-1) + 2) = maxvnew + .1*vscalenew;
        end
        if ( minvnew < Clim(2*(k-1) + 1) )
            vscalenew = max(abs(minvnew), maxvnew);
            Clim(2*(k-1) + 1) = minvnew - .1*vscalenew;
        end
    end

    % Interpolate each variable on a finer grid:
    if ( Nplot > N )
        vvv = interp2(ll, tt, vv, lll, ttt, 'spline');
    else
        vvv = vv;
    end

    % Update each variable:
    set(p{1,k}, 'cdata', vvv)
    set(p{1,k}.Parent, 'clim', [Clim(2*(k-1) + 1), Clim(2*(k-1) + 2)])

    % Update title:
    titleString = sprintf('n = m = %i (DoFs = %i), dt = %1.1e, t = %.4f', N, ...
        nVars*N^2, dt, t);
    set(p{2,k}, 'String', titleString)
    drawnow

end

% Update outputs:
opts{1} = Clim;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
if(exist('D:\Datasets\UNBC/', 'file'))
    UNBC_dir = 'E:\Datasets\UNBC/';
elseif(exist('E:\Datasets\UNBC/', 'file'))
    UNBC_dir = 'E:\Datasets\UNBC/';
else
    fprintf('UNBC location not found (or not defined)\n');
end

all_recs = {'042-ll042', '043-jh043', '047-jl047', '048-aa048', '049-bm049',...
            '052-dr052', '059-fn059', '064-ak064', '066-mg066', '080-bn080',...
            '092-ch092', '095-tv095', '096-bg096', '097-gf097', '101-mg101',...
            '103-jk103', '106-nm106', '107-hs107', '108-th108', '109-ib109',...
            '115-jy115', '120-kz120', '121-vw121', '123-jh123', '124-dn124'};

devel_recs = all_recs(1:5:25);
train_recs = setdiff(all_recs, devel_recs);

all_aus = [4, 6, 7, 9, 10, 12, 20, 25, 26, 43];

author: TadasBaltrusaitis

repo: OpenFace


text:
function X = full(X)
%FULL Convert to a (dense) tensor.
%
%   FULL(X) returns X, i.e., does nothing.
%
%   See also TENSOR, SPTENSOR/FULL, KTENSOR/FULL, TTENSOR/FULL.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

return;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [y,varname]=xml2cell(filename_or_xml_string)

%XML2CELL reads non-MbML compliant xmlfile into matlab nested cell arrays
%
% Syntax [y,varname]=xmlfile2cell(filename_or_xml_string)
%
% Description:
%   1. Convert any non-MbML xml into MbML compliant string
%   2. Stores individual structures as nested cell arrays
%
% If it cannot be garanteed that non-MbML compliant XML has an
% internal referential consistency for convenient conversion to structures,
% this function builds the m-variable object model (MOM) as nested
% cell arrays. This approach ignores the possible dimensionality of the
% object and stores each entry, in a single cell, nested at the
% appropriate level.
%
% Note 1 : if your XML string is MbML compliant use XML2MAT instead
% Note 2 : if your XML structure has consistent dimensionality use XML2STRUCT instead
%
% See also: XML2STRUCT, XML2MAT
%
% Jonas Almeida, almeidaj@musc.edu, 19 May 2003, MAT4NAT Tbox

% is this a xml string or xml filename ?
if filename_or_xml_string(1)=='<'
    y=filename_or_xml_string;
else
    y=strrep(file2str(filename_or_xml_string),'''','''''');
end

% convert first to MbML compliant string and then onto an m-variable
[y,varname]=xml2mat(mbmling(y,1));

author: Sable

repo: mcbench-benchmarks


text:
function [ ncycle, t, index ] = perm_to_cycle ( n, p )

%*****************************************************************************80
%
%% PERM_TO_CYCLE converts a permutation from array to cycle form.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the number of values being permuted.
%    N must be positive.
%
%    Input, integer P(N), describes the permutation using a
%    single array.  For each index I, I -> P(I).
%
%    Output, integer NCYCLE, the number of cycles.
%    1 <= NCYCLE <= N.
%
%    Output, integer T(N), INDEX(N), describes the permutation
%    as a collection of NCYCLE cycles.  The first cycle is
%    T(1) -> T(2) -> ... -> T(INDEX(1)) -> T(1).
%

%
%  Check.
%
  missing = perm_check ( n, p );

  if ( missing ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PERM_TO_CYCLE - Fatal error!\n' );
    fprintf ( 1, '  The input array is illegal.\n' );
    fprintf ( 1, '  Missing element = ', missing );
    error ( 'PERM_TO_CYCLE - Fatal error!' );
  end
%
%  Initialize.
%
  ncycle = 0;
  index(1:n) = 0;
  t(1:n) = 0;
  nset = 0;
%
%  Find the next unused entry.
%
  for i = 1 : n

    if ( 0 < p(i) )

      ncycle = ncycle + 1;
      index(ncycle) = 1;

      nset = nset + 1;
      t(nset) = p(i);
      p(i) = - p(i);

      while ( 1 )

        j = t(nset);

        if ( p(j) < 0 )
          break
        end

        index(ncycle) = index(ncycle) + 1;

        nset = nset + 1;
        t(nset) = p(j);
        p(j) = - p(j);

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [outputs,model,accuracies,accuracies_locs,bestlambda] = hkl_fold(nfold,typefold,X,Y,lambdas,loss,kernel,kernel_params,varargin);
%%%%%%%%%%%%%%%%%%%%%%%
% HKL K-FOLD
%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% required parameter
% nfold                 5 or 10
% typefold              'same' or 'scaled' -> using same lambda or same lambda/n
% X                     input data: n x p matrix (n=number of observations)
%						or kernel matrices ( n(n+1)/2 x p x q single)
% Y                     responses ( n x 1 matrix )
%                       NB: for classification in {0,1} (and not in {-1,1})
% lambdas               regularization parameters (might be vector or single number)
%						better to start first with large values of lambdas
% loss                  loss function ('square','logistic')
% kernel                kernel to be decomposed (see list below)
% kernel_params         kernel parameters
% varargin              see hkl.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% find out if display is suppressed...
f = strcmp('display',varargin(1:2:end));
if any(f)
    display = varargin{find(f)*2};
else
    display = 1;
end

if display fprintf('all folds\n'); end
% first perform learning on the full dataset%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[outputs,model,accuracies] = hkl(X,Y,lambdas,loss,kernel,kernel_params,varargin{:});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% centering and scaling is done separatetly for each fold
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% cut into pieces
n = length(Y);

% cut into folds
folds = cell(1,nfold);
compfolds = cell(1,nfold);
for ifold = 1:nfold
    if ifold<nfold
        folds{ifold} = floor( (ifold-1)*n/nfold + 1): floor( ifold*n/nfold );
    else
        folds{ifold} = floor( (ifold-1)*n/nfold + 1):n;
    end
    compfolds{ifold}=1:n;
    compfolds{ifold}(folds{ifold}) = [];

end

% define lambdas on the fold
switch typefold,
    case 'same'
        lambdas_fold = lambdas;
    case 'scaled'
        lambdas_fold = lambdas / ( 1- 1/nfold) ;
end

for ifold = 1:nfold
    if display fprintf('fold %d\n',ifold); end
    Yloc = Y(compfolds{ifold});
    Ytestloc = Y(folds{ifold});

    if ~strcmp(kernel,'base kernels') || ~strcmp(kernel,'base kernels-mkl') || ~strcmp(kernel,'base kernels-bimkl')
        Xloc = X(compfolds{ifold},:);
        Xtestloc = X(folds{ifold},:);
    else
        error('not implemented yet')

    end

    [outputs_loc,model_loc,accuracies_loc] = hlp_diskcache('predictivemodels',@hkl,Xloc,Yloc,lambdas_fold,loss,kernel,kernel_params,varargin{:},'Xtest',Xtestloc,'Ytest',Ytestloc);

    accuracies_locs{ifold} = accuracies_loc;
end

temp = [];
for ifold=1:nfold, temp = [ temp; accuracies_locs{ifold}.testing_error]; end;

[a,bestlambda] = min(mean(temp,1));
try
estimated_best_error = accuracies.testing_error(bestlambda)
best_error = min(accuracies.testing_error)
catch,end

repo: REST


text:
function [coef]=comp_dgtreal_fb(f,g,a,M)
%COMP_DGTREAL_FB  Filter bank DGT
%   Usage:  c=comp_dgt_fb(f,g,a,M,boundary);
%
%   This is a computational routine. Do not call it directly.

%   See help on DGT.

%   AUTHOR : Peter L. Søndergaard.

% Calculate the parameters that was not specified.
L=size(f,1);
N=L/a;
gl=length(g);
W=size(f,2);      % Number of columns to apply the transform to.
glh=floor(gl/2);  % gl-half
M2=floor(M/2)+1;

% Conjugate the window here.
g=conj(fftshift(g));

coef=zeros(M,N,W,assert_classname(f,g));

% Replicate g when multiple columns should be transformed.
gw=repmat(g,1,W);

% ----- Handle the first boundary using periodic boundary conditions. ---
for n=0:ceil(glh/a)-1

    % Periodic boundary condition.
    fpart=[f(L-(glh-n*a)+1:L,:);...
           f(1:gl-(glh-n*a),:)];

  fg=fpart.*gw;

  % Do the sum (decimation in frequency, Poisson summation)
  coef(:,n+1,:)=sum(reshape(fg,M,gl/M,W),2);

end;

% ----- Handle the middle case. ---------------------
for n=ceil(glh/a):floor((L-ceil(gl/2))/a)

  fg=f(n*a-glh+1:n*a-glh+gl,:).*gw;

  % Do the sum (decimation in frequency, Poisson summation)
  coef(:,n+1,:)=sum(reshape(fg,M,gl/M,W),2);
end;

% ----- Handle the last boundary using periodic boundary conditions. ---
for n=floor((L-ceil(gl/2))/a)+1:N-1

    % Periodic boundary condition.
    fpart=[f((n*a-glh)+1:L,:);... %   L-n*a+glh elements
           f(1:n*a-glh+gl-L,:)];  %  gl-L+n*a-glh elements

    fg=fpart.*gw;

    % Do the sum (decimation in frequency, Poisson summation)
    coef(:,n+1,:)=sum(reshape(fg,M,gl/M,W),2);
end;

% --- Shift back again to make it a frequency-invariant system. ---
for n=0:N-1
  coef(:,n+1,:)=circshift(coef(:,n+1,:),n*a-glh);
end;

coef=fftreal(coef);
coef=reshape(coef,M2,N,W);

%c=c(1:M2,:);

% Simple code using a lot of circshifts.
% Move f initially so it lines up with the initial fftshift of the
% window
%f=circshift(f,glh);
%for n=0:N-1
  % Do the inner product.
  %fg=circshift(f,-n*a)(1:gl,:).*gw;

  % Periodize it.
  %fpp=zeros(M,W);
  %for ii=0:gl/M-1
    %  fpp=fpp+fg(ii*M+1:(ii+1)*M,:);
    %end;
%  fpp=sum(reshape(fg,M,gl/M,W),2);

  % Shift back again.
%  coef(:,n+1,:)=circshift(fpp,n*a-glh); %),M,1,W);

%end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
classdef KernelPCA < handle
    %{
        Kernel Principal component analysis (KPCA)

        Version 2.2, 14-MAY-2021
        Email: iqiukp@outlook.com
    -------------------------------------------------------------------
    %}

    properties
        data
        label
        numComponents
        explainedLevel
        kernelFunc = Kernel('type', 'gaussian', 'gamma', 0.5)
        lambda
        coefficient %　principal component coefficients
        score %　principal component scores.
        cumContribution % cumulative contribution rate
        newData % Transform the mapping data back to original space
        T2
        SPE
        T2Limit
        SPELimit
        numSPEAlarm
        numT2Alarm
        accuracySPE
        accuracyT2
        eigenvalueTolerance = 1e-8 % tolerance for eigenvalues
        alpha = 1 % hyperparameter of the ridge regression that learns the reconstruction
        theta = 0.7 % experience parameter of fault diagnosis
        significanceLevel = 0.95
        display = 'on'
        temporary
        diagnosis = []
        runningTime
    end

    properties (Dependent)
        numSamples
        numFeatures
    end

    methods
        function obj = KernelPCA(parameter)
            name_ = fieldnames(parameter);
            for i = 1:size(name_, 1)
                obj.(name_{i, 1}) = parameter.(name_{i, 1});
            end
            KernelPCAOption.checkInputForDiagnosis(obj);
        end

        function obj = train(obj, data)
            tStart = tic;
            obj.data = data;
            obj.label = ones(obj.numSamples, 1);

            % compute the kernel matrix
            K = obj.kernelFunc.computeMatrix(obj.data, obj.data);
            % centralize the kernel matrix
            unit = ones(obj.numSamples, obj.numSamples)/obj.numSamples;
            K_c = K-unit*K-K*unit+unit*K*unit;
            % compute the eigenvalues and eigenvectors
            [obj.coefficient, obj.lambda] = obj.computeEigenvalue(K_c);
            %
            obj.score = K_c* obj.coefficient(:, 1:obj.numComponents);
            obj.newData = obj.reconstruct;
            obj.temporary.K = K;
            obj.temporary.K_c = K_c;
            obj.temporary.unit = unit;
            obj.computeControlLimit;

            % compute accuracy
            T2AlarmIndex = find(obj.T2 > obj.T2Limit);
            SPEAlarmIndex = find(obj.SPE > obj.SPELimit);
            obj.numSPEAlarm = length(SPEAlarmIndex);
            obj.numT2Alarm = length(T2AlarmIndex);

            label_ = obj.label;
            label_(SPEAlarmIndex) = -1;
            obj.accuracySPE = sum(label_ == obj.label)/obj.numSamples;
            label_ = obj.label;
            label_(T2AlarmIndex) = -1;
            obj.accuracyT2 = sum(label_ == obj.label)/obj.numSamples;
            %
            obj.runningTime = toc(tStart);
            if strcmp(obj.display, 'on')
                KernelPCAOption.displayTrain(obj)
            end
        end

        function results = test(obj, data, varargin)
            % test the model from the given data
            tStart = tic;
            results.evaluation = 'off';
            if nargin == 3
                results.evaluation = 'on';
                testLabel = varargin{1};
            end
            Kt = obj.kernelFunc.computeMatrix(data, obj.data);
            % centralize the kernel matrix
            unit = ones(size(data, 1), obj.numSamples)/obj.numSamples;
            Kt_c = Kt-unit*obj.temporary.K-Kt*obj.temporary.unit+unit*obj.temporary.K*obj.temporary.unit;
            %
            results.numSamples = size(data, 1);
            results.data = data;
            results.score = Kt_c*obj.coefficient(:, 1:obj.numComponents);

            % compute Hotelling's T2 statistic
            results.T2 = diag(results.score/diag(obj.lambda(1:obj.numComponents))*results.score');
            % compute the squared prediction error (SPE)
            results.SPE = sum((Kt_c*obj.coefficient).^2, 2)-sum(results.score.^2 , 2);

            % compute accuracy
            results.T2AlarmIndex = find(results.T2 > obj.T2Limit);
            results.SPEAlarmIndex = find(results.SPE > obj.SPELimit);
            if strcmp(results.evaluation, 'on')
                label_ = ones(size(results.data, 1), 1);
                label_(results.SPEAlarmIndex) = -1;
                results.accuracySPE = sum(label_ == testLabel)/results.numSamples;
                label_ = ones(size(results.data, 1), 1);
                label_(results.T2AlarmIndex) = -1;
                results.accuracyT2 = sum(label_ == testLabel)/results.numSamples;
            end
            results.numSPEAlarm = length(results.SPEAlarmIndex);
            results.numT2Alarm = length(results.T2AlarmIndex);
            results.temporary.Kt = Kt;
            results.runningTime = toc(tStart);

            if strcmp(obj.display, 'on')
                KernelPCAOption.displayTest(results)
            end

            % fault diagnosis
            if strcmp(obj.diagnosis.switch, 'on')
                results = obj.diagnose(results);
            end
        end

        function newData = reconstruct(obj)
            % Transform the mapping data back to original space.
            % References
            % ----------
            % Bakır G H, Weston J, Schölkopf B. Learning to find pre-images[J].
            % Advances in neural information processing systems, 2004, 16: 449-456.

            K_1 =  obj.kernelFunc.computeMatrix(obj.score, obj.score);
            K_1_ = K_1;
            for i = 1:obj.numSamples
                K_1(i, i) = K_1(i, i)+obj.alpha;
            end
            dual_coef = mldivide(K_1, obj.data);
            K_2 =  K_1_;
            newData = K_2*dual_coef;
        end

        function [coefficient, lambda] = computeEigenvalue(obj, K_c)
            % compute the eigenvalues and eigenvectors
            rng('default')
            [V, D, ~] = svd(K_c/obj.numSamples, 'econ');
            % ill-conditioned matrix
            if ~(isreal(V)) || ~(isreal(D))
                V = real(V);
                D = real(D);
            end
            lambda_ = diag(D);
            obj.cumContribution = cumsum(lambda_/sum(lambda_));

            if isempty(obj.numComponents)
                obj.numComponents = obj.numFeatures;
            else
                if obj.numComponents >= 1
                    obj.numComponents = obj.numComponents;
                else
                    obj.explainedLevel = obj.numComponents;
                    obj.numComponents = find(obj.cumContribution >= obj.numComponents, 1, 'first');
                end
            end
            lambda = lambda_;
            try
                coefficient = V./sqrt(obj.numSamples*lambda_)';
            catch
                coefficient = zeros(obj.numSamples, obj.numSamples);
                for i = 1:obj.numSamples
                    coefficient(:, i) = V(:, i)/sqrt(obj.numSamples*lambda_(i, 1));
                end
            end
        end

        function computeControlLimit(obj)
            % compute the squared prediction error (SPE)
            temp1 = obj.temporary.K_c*obj.coefficient;
            temp2 = obj.temporary.K_c*obj.coefficient(:, 1:obj.numComponents);
            obj.SPE = sum(temp1.^2, 2)-sum(temp2.^2, 2);
            obj.T2 = diag(obj.score/diag(obj.lambda(1:obj.numComponents))*obj.score');

            % compute the T2 limit (the F-Distribution)
            k = obj.numComponents*(obj.numSamples-1)/(obj.numSamples-obj.numComponents);
            obj.T2Limit = k*finv(obj.significanceLevel, obj.numComponents, obj.numSamples-obj.numComponents);

            % compute the SPE limit (the Chi-square Distribution)
            a = mean(obj.SPE);
            b = var(obj.SPE);
            g = b/2/a;
            h = 2*a^2/b;
            obj.SPELimit = g*chi2inv(obj.significanceLevel, h);
        end

        function results = diagnose(obj, results, varargin)
            % falut diagnosis
            tStart = tic;
            fprintf('\n')
            fprintf('*** Fault diagnosis ***\n')
            fprintf('Fault diagnosis start...\n')
            results.diagnosis = obj.diagnosis;
            data_ = results.data;
            results.diagnosis.data = data_(results.diagnosis.start:results.diagnosis.end, :);
            % contribution plots of train data
            if ~exist('.\data', 'dir')
                mkdir data;
            end
            file_ = dir('.\data\*.mat');
            name_ = {file_(1:length(file_)).name}';
            if ismember('diagnosis.mat', name_)
                load('.\data\diagnosis.mat', 'tmp_')
                tmp__ = KernelPCAOption.saveCheckObj(obj);
                if isequal(tmp__, tmp_)
                    load('.\data\diagnosis.mat', 'T2CpsTrain', 'SPECpsTrain')
                else
                    [T2CpsTrain, SPECpsTrain] = obj.computeContribution(results, 'train');
                    tmp_ = KernelPCAOption.saveCheckObj(obj);
                    save('.\data\diagnosis.mat', 'T2CpsTrain', 'SPECpsTrain', 'tmp_')
                end
            else
                [T2CpsTrain, SPECpsTrain] = obj.computeContribution(results, 'train');
                tmp_ = KernelPCAOption.saveCheckObj(obj);
                save('.\data\diagnosis.mat', 'T2CpsTrain', 'SPECpsTrain', 'tmp_')
            end

            % contribution plots of test data
            [T2CpsTest, SPECpsTest] = obj.computeContribution(results, 'test');
            % normalize the contribution plots
            T2CpsTrainMu = mean(T2CpsTrain);
            T2CpsTrainStd = std(T2CpsTrain);

            try
                T2Cps = bsxfun(@rdivide, bsxfun(@minus, T2CpsTest, T2CpsTrainMu), T2CpsTrainStd);
            catch
                mu_array = repmat(T2CpsTrainMu, size(T2CpsTest,1), 1);
                st_array = repmat(T2CpsTrainStd, size(T2CpsTest,1), 1);
                T2Cps = (T2CpsTest-mu_array)./st_array;
            end
            SPECpsTrainMu = mean(SPECpsTrain);
            SPECpsTrainStd = std(SPECpsTrain);
            try
                SPECps = bsxfun(@rdivide, bsxfun(@minus, SPECpsTest, SPECpsTrainMu), SPECpsTrainStd);
            catch
                mu_array = repmat(SPECpsTrainMu, size(SPECpsTest,1), 1);
                st_array = repmat(SPECpsTrainStd, size(SPECpsTest,1), 1);
                SPECps = (SPECpsTest-mu_array)./st_array;
            end

            % store the results
            results.diagnosis.T2Cps = T2Cps;
            results.diagnosis.SPECps = SPECps;

            %
            T2Cps_ = mean(abs(T2Cps), 1);
            results.diagnosis.meanT2Cps = T2Cps_/sum(T2Cps_, 2);

            SPECps_ = mean(abs(SPECps), 1);
            results.diagnosis.meanSPECps = SPECps_/sum(SPECps_, 2);

            %
            [value, index] = sort(results.diagnosis.meanT2Cps, 'descend');
            results.diagnosis.faultVariabeT2.value = value;
            results.diagnosis.faultVariabeT2.index = index;

            [value, index] = sort(results.diagnosis.meanSPECps, 'descend');
            results.diagnosis.faultVariabeSPE.value = value;
            results.diagnosis.faultVariabeSPE.index = index;

            results.diagnosis.runningTime = toc(tStart);
            if strcmp(obj.display, 'on')
                KernelPCAOption.displayDiagnose(results)
            end
        end

        function [T2Cps, SPECps] = computeContribution(obj, result, type)

            %  Compute the Contribution Plots (CPs)
            %
            %  Reference
            %  [1]  Deng X, Tian X. A new fault isolation method based on unified
            %       contribution plots[C]//Proceedings of the 30th Chinese Control
            %       Conference. IEEE, 2011: 4280-4285.
            % -------------------------------------------------------------------
            %  Thanks for the code provided by Rui.
            % --------------------------------------------------------------------

            data_ = obj.data;
            switch type
                case 'train'
                    Kt = obj.temporary.K;
                    Y = data_;
                case 'test'
                    Kt = result.temporary.Kt;
                    Y = result.diagnosis.data;
            end

            K = obj.temporary.K;
            M = size(data_, 1);
            [Mt, d] = size(Y);

            A_T2 = obj.coefficient(:, 1:obj.numComponents)*...
                diag(obj.lambda(1:obj.numComponents))^(-1)*...
                obj.coefficient(:, 1:obj.numComponents)';

            A_SPE = obj.coefficient(:, 1:obj.numComponents)*...
                obj.coefficient(:, 1:obj.numComponents)';
            newY = Y*obj.theta;

            % initialization
            Knew = zeros(Mt, M);
            Knew_d1 = zeros(1, M);
            Knew_d2 = zeros(Mt, M);
            T2Cps = zeros(Mt, d);
            SPECps = zeros(Mt, d);
            Knew_s = zeros(Mt, M);
            sigma = sqrt(1/2/obj.kernelFunc.gamma);

            % compute contribution of statistic
            for i = 1:Mt
                for j = 1:d
                    for k = 1:M
                        Knew(i, k) = Kt(i, k);
                        Knew_d1(k) = Knew(i, k)*2*obj.theta*(newY(i, j)-data_(k, j))/(-sigma^2); % derivative
                        Knew_d2(i, k) = -2*Knew_d1(k);
                    end
                    Knew_d1_s = Knew_d1-ones(1, M)*mean(Knew_d1);
                    Knew_s(i, :) = Knew(i, :)-ones(1, M)*K/M-Knew(i, :)*ones(M) ...
                        /M+ones(1, M)/M*K*ones(M)/M;
                    % contribution of T2
                    T2Cps(i, j) = Y(i, j)*(Knew_d1_s*A_T2*Knew_s(i, :)' ...
                        +Knew_s(i, :)*A_T2*Knew_d1_s');
                    % contribution of SPE
                    SPECps(i, j)= Y(i, j)*mean(Knew_d2(i, :))-Y(i, j) ...
                        *(Knew_d1_s*A_SPE*Knew_s(i, :)'+Knew_s(i, :)*A_SPE*Knew_d1_s');
                end
            end
        end

        function numSamples = get.numSamples(obj)
            numSamples= size(obj.data, 1);
        end

        function numFeatures = get.numFeatures(obj)
            numFeatures= size(obj.data, 2);
        end
    end
end

author: iqiukp

repo: KPCA-MATLAB

save_path: github-repos/MATLAB/iqiukp-KPCA-MATLAB


text:
classdef GridMap < handle
    properties
        boundary;       % size of Metric Map(leftdownPoint + rightupPoint)
        blocks;         % coordinate of obstacles on the Metric Map
        res;            % resolution
        margin;
        maxIndex;       % max index of ijk of GridMap
        occ_map;        % flag obstacles on the GridMap
    end
    methods
        function obj = GridMap()
        end

        function load_map(obj, filename, xy_res, z_res, margin)
            fid = fopen(filename);
            tline = fgets(fid);
            obj.blocks = [];
            while ischar(tline)
                % skip empty line and comments
                if numel(tline) > 1 & tline(1)~='#'
                    % convert char array to string
                    if strncmpi(tline, 'boundary', 8)
                        boundarys = strread(tline(9:end));
                    end
                    if strncmpi(tline, 'block', 5)
                        block = strread(tline(6:end));
                        assert(size(block,2) == 9);
                        obj.blocks = [obj.blocks; block];
                    end
                end
                tline = fgets(fid);
            end
            obj.boundary.ld = boundarys(1:3);
            obj.boundary.ru = boundarys(4:6);
            obj.res = [xy_res xy_res z_res];
            obj.margin = margin;
            obj.maxIndex = ceil(abs( (obj.boundary.ru - obj.boundary.ld)./obj.res ));
        end

        % obstacles = block + margin
        function flag_obstacles(obj)
            obj.occ_map = false(obj.maxIndex);
            for i = 1 : size(obj.blocks,1)
                block = obj.blocks(i, :);
                leftdownPoint = [min(block(1), block(4)) min(block(2), block(5)) min(block(3), block(6))];
                rightUpPoint = [max(block(1), block(4)) max(block(2), block(5)) max(block(3), block(6))];
                obj.blocks(i, 1:3) = leftdownPoint;
                obj.blocks(i, 4:6) = rightUpPoint;
                xyzs = obj.points_to_idx(leftdownPoint - obj.margin);
                xyze = obj.points_to_idx(rightUpPoint + obj.margin);
                obj.occ_map(xyzs(1):xyze(1), xyzs(2):xyze(2), xyzs(3):xyze(3)) = 1;
            end
        end

        function ijk = points_to_idx(obj, xyz)
            assert(size(xyz,1) > 0);
            ijk = floor((xyz - obj.boundary.ld)./obj.res + 1);
            ijk = min(max(ijk, 1), obj.maxIndex);
        end

        function xyz = idx_to_points(obj, ijk)
            assert(size(ijk,1) > 0);
            xyz = obj.boundary.ld + (ijk - 1 + 0.5).* obj.res;
        end

        function ret = idOverflowCheck(obj, i, j, k)
            ret = (i < 1 || i > obj.maxIndex(1) || j < 1 || j > obj.maxIndex(2) || k < 1 || k > obj.maxIndex(3));
        end

        % check points collide
        function C = collide(obj, points)
            ijk = obj.points_to_idx(points);
            C = false(size(ijk, 1), 1);
            for i = 1 : size(ijk, 1)
                id = ijk(i,:);
                C(i) =  obj.occ_map(id(1), id(2), id(3)) ~= 0;
            end
        end

        %% generate intermedia points
        function pts = generateiIntermediaPoints(obj, start_pos, end_pos)
            %% generate Tiny direction vector
            dd = end_pos - start_pos;
            dir = obj.res(1)*(dd/norm(dd, 2))/(100);

            %% generate intermedia points
            lenxy = sqrt(sum((start_pos(1:2) - end_pos(1:2)).^2));
            lenz = sqrt(sum((start_pos(3) - end_pos(3)).^2));
            n_pts = int32(max(lenxy/obj.res(1), lenz/obj.res(3))) + 3;
            pts = [linspace(start_pos(1), end_pos(1), n_pts);
                   linspace(start_pos(2), end_pos(2), n_pts);
                   linspace(start_pos(3), end_pos(3), n_pts)]';

            %% mid-point + Tiny direction vector
            for j = 2 : size(pts,1) - 1
                m1 = mod(pts(j,1), obj.res(1));
                m2 = mod(pts(j,2), obj.res(1));
                m3 = mod(pts(j,3), obj.res(1));
                if(m1*m2*m3 == 0)
                    pts(j,:) = pts(j,:) + dir;
                end
            end
        end
    end
end

author: LenaShengzhen

repo: AerialRobotics


text:
%  sabs smooth approximation of abs
%
%  ::
%
%
%  Args:
%
%  Returns:
%     :
%
%  Note:
%
%  Example:
%
%     See also:
%

author: jmaih


text:
function out=fieldnames(a)
% FIELDNAMES Get POLYNOMIAL object property names.
%
% NAMES=FIELDNAMES(P) returns a cell array of strings containing the
% names of the properties associated with uncertain polynomial object, P.
%
% See also GET, SET, POLYNOMIAL

% 11/12/2002: PJS  Initial Coding
%    This function allows tab completion of a public data fields.

aprop = PVget(a,'PropNames');
out = aprop.GPropNames;

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [ n_data, m, n, fmn ] = i4_rise_values ( n_data )

%*****************************************************************************80
%
%% I4_RISE_VALUES returns values of the integer rising factorial function.
%
%  Discussion:
%
%    The integer rising factorial function is sometimes symbolized by (m)_n.
%
%    The definition is
%
%      (m)_n = (m-1+n)! / (m-1)!
%            = ( m ) * ( m + 1 ) * ( m + 2 ) ... * ( m - 1 + n )
%            = Gamma ( m + n ) / Gamma ( m )
%
%    We assume 0 <= N <= M.
%
%    In Mathematica, the function can be evaluated by:
%
%      Pochhammer[m,n]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 December 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer M, N, the arguments of the function.
%
%    Output, integer FMN, the value of the function.
%
  n_max = 8;

  fmn_vec = [ ...
     1, 5, 30, 210, 1680, ...
     15120, 151200, 1, 10, 4000, ...
     110, 6840, 840, 970200, 5040 ];

  m_vec = [ ...
    5, 5, 5, 5, 5, ...
    5, 5, 50, 10, 4000, ...
    10, 18, 4, 98, 1 ];

  n_vec = [ ...
    0, 1, 2, 3, 4, ...
    5, 6, 0, 1, 1, ...
    2, 3, 4, 3, 7 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    m = 0;
    n = 0;
    fmn = 0;
  else
    m = m_vec(n_data);
    n = n_vec(n_data);
    fmn = fmn_vec(n_data);
  end

  return
end


text:
% Prepare testing data to quicken testing
init;

model = exp_params.model;

% Get test image path
im_dir = [VOC07PATH 'JPEGImages/'];
fid = fopen([VOC07PATH 'ImageSets/Main/test.txt']);
test_imgs = textscan(fid, '%s');
test_imgs = test_imgs{1};
num_test = size(test_imgs, 1);

% Network to use
caffe_params.model = exp_params.model;
caffe_params.model_file = exp_params.model_file;
caffe_params.def_file = exp_params.def_file;
caffe_params.device = exp_params.device;

% Load bounding boxes
load('data/SelectiveSearchVOC2007test.mat');
bbox = boxes(:);

% Initialization
max_img = 10;             % maximum number of images in the .mat file
tests = cell(1);
counter = 1;
file_counter = 1;

for ii=1:num_test
    disp(['Image: ' num2str(ii) '/' num2str(num_test)]);
    im_path = [im_dir test_imgs{ii} '.jpg'];
    im = imread(im_path);

    %boxes = selective_search(im, 'single');
    index = find(strcmp(images, test_imgs{ii}));
    boxes = bbox{index};
    num_boxes = size(boxes, 1);
    disp(['Number of boxes ' num2str(num_boxes)]);

    tic
    test = cell(num_boxes, 1);
    for jj=1:num_boxes
        box = boxes(jj, :);
        region = im(box(1):box(3), box(2):box(4), :);
        rep = extract_caffe_feature(region, caffe_params);

        test{jj} = struct;
        test{jj}.feature = mean(rep, 2);
        test{jj}.bbox = box;
    end
    tests{counter} = test;
    counter = counter + 1;
    toc

    if counter > max_img
        disp(['Save to file']);
        save(['data/' model '/test_' num2str(file_counter) '.mat'], 'tests', '-v7.3');
        counter = 1;
        file_counter = file_counter + 1;
        clear tests; tests = cell(max_img, 1);
    end
end

save(['data/' model '/test_' num2str(file_counter) '.mat'], 'tests', '-v7.3');


text:
% PURPOSE: subroutine for pop_continuousartdet.m
%
% Author: Javier Lopez-Calderon
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2011

function [WinRej chanrej] = movesegments(WinRej, chanrej, winoffsetsam, pnts)
fprintf('\nWARNING: Marked segments will be displaced in %g samples.\n\n', winoffsetsam);
WinRej = WinRej + winoffsetsam;
WinRej(WinRej<1) = 1;
WinRej(WinRej>pnts) = pnts;
[WinRej chanrej] = discardshortsegments(WinRej, chanrej, 0, 0);

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function g = uminus(f)
%- Unary minus of a BALLFUNV.
%   -F returns the BALLFUNV negated componentwise.
%
%   UMINUS(F) is called by the syntax -F.
%
% See also UPLUS.

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

if isempty( f )
    g = ballfunv();
    return
end

F = f.comp;
g = [ -F{1} ; -F{2} ; -F{3} ];
end


text:
function canlab_toolbox_setup
% canlab_toolbox_setup
% ------------------------------------------------------------------------
% This function is run with no input arguments. It checks for a number of
% CANlab repositories, and gives you the option to download them from
% Github, add them to the Matlab path with subfolders, and save the Matlab
% path. Run it from the directory ABOVE your individual repositories. e.g.,
% if you have a local directory for your Github repositories, run it from
% there. This function assumes you have all CANlab repositories in the same master
% folder.
%
% The following repositories are included. This is a subset of the
% repositories available, but contains must of the publicly useful ones:
% ------------------------------------------------------------------------
% CANlab core tools - Object-oriented interactive analysis and tools used in other functions
% CANlab help examples - Example scripts and html output, second-level analysis batch scripts
% MKDA meta-analysis toolbox - Coordinate-based meta-analysis of neuroimaging data
% FMRI_simulations - Misc tools for simulating and testing power and image data
% Private masks and signatures - Selected models shared on request
% M3 Mediation toolbox - Single- and Multi-level mediation tools
% Neuroimaging_Pattern_Masks - meta-analysis masks, signatures, parcellations
% Robust regression toolbox - Voxel-wise robust regression for neuroimaging data
% CanlabPrivate - Private repository for lab members and collaborators
%
% July 2018, Tor Wager

main_repository_dir = pwd;

% If you drag and drop the script in, Matlab changes to the folder with the script automatically,
% then changes back.  If so, assume main dir is 2 levels above this.
% If you have it on your path and run it by typing the script name only, it
% will not change directories.
[dd, ff] = fileparts(main_repository_dir);
if strcmp(ff, 'CanlabCore') && exist(fullfile(pwd, 'canlab_toolbox_setup.m'), 'file')
    % We are in CanlabCore subdirectory

    main_repository_dir = fileparts(fileparts(main_repository_dir));
end

dashes = '-----------------------------------------';
fprintf('\n%s\nChecking and/or installing CANlab tools in:\n%s\n%s\n', dashes, main_repository_dir, dashes);

repo_descrip = {'CANlab core tools - Object-oriented interactive analysis and tools used in other functions' ...
    'CANlab help examples - Example scripts and html output, second-level analysis batch scripts' ...
    'MKDA meta-analysis toolbox - Coordinate-based meta-analysis of neuroimaging data' ...
    'FMRI_simulations - Misc tools for simulating and testing power and image data' ...
    'Private masks and signatures - Selected models shared on request' ...
    'M3 Mediation toolbox - Single and Multi-level mediation tools' ...
    'Neuroimaging_Pattern_Masks - meta-analysis masks, signatures, parcellations' ...
    'Robust regression toolbox - Voxel-wise robust regression for neuroimaging data' ...
    'CanlabPrivate - Private repository for lab members and collaborators'};

% Find these files on disk, which will tell us where toolbox folders are:

key_files_to_find = {'fmri_data.m' ...
    'a2_second_level_toolbox_check_dependencies.m' ...
    'i_density.m' ...
    'power_figure3_num_comparisons.m' ...
    'apply_nps.m' ...
    'mediation_brain.m' ...
    'apply_all_signatures.m' ...
    'robust_results_batch.m' ...
    'power_calc.m'};

group_url = 'https://github.com/canlab/';

repo_names = {'CanlabCore' ...
    'CANlab_help_examples' ...
    'Canlab_MKDA_MetaAnalysis' ...
    'FMRI_simulations' ...
    'MasksPrivate' ...
    'MediationToolbox' ...
    'Neuroimaging_Pattern_Masks' ...
    'RobustToolbox' ...
    'CanlabPrivate'};

n_repos = length(key_files_to_find);

% ------------------------------------------------------------------------
% List repos
% ------------------------------------------------------------------------

disp('Looking for CANlab repositories and adding to Matlab path:');
disp(' ');

for i = 1:n_repos

    fprintf('%s\t%s\n', repo_descrip{i}, repo_names{i});

end

fprintf('\n%s\n', dashes);

% ------------------------------------------------------------------------
% Find each file
% Get enclosing folder
% Add with subfolders to path
% ------------------------------------------------------------------------

was_missing = true(1, n_repos);

for i = 1:n_repos

    fprintf('\n%s\n', dashes);

    fprintf('\nRepository: %s\nSite:%s\n', repo_descrip{i}, repo_names{i});

    myrepodir = fullfile(main_repository_dir, repo_names{i});

    % Future update:
    % Check if already on path. If so, add with subfolders to be sure all
    % are added. <could use regexp to find base dir and use that>
    % ------------------------------------------------------------------------
    %
    %     checkfile = which([repo_names{i} filesep '@fmri_data' filesep 'fmri_data.m']);
    %     toolboxdir = fileparts(fileparts(checkfile));
    %
    %     if ~exist(toolboxdir, 'dir')
    %
    %         disp('Cannot find CANlabCore toolbox');
    %         disp('You need the CANlab CanlabCore Github repository on your path, with subfolders, to run these scripts.')
    %
    %     else
    %         disp('Found CANlabCore toolbox');
    %         g = genpath(toolboxdir);
    %         addpath(g);
    %     end

    disp('Finding repo: Running command:')
    disp(['find ' myrepodir ' -name "' key_files_to_find{i} '"']);

    [status,result] = system(['find ' myrepodir ' -name "' key_files_to_find{i} '"']);

    if ~status && ~isempty(result)

        was_missing(i) = false;

        % Get enclosing folder
        % mydir = fileparts(fileparts(result(1, :)));

        % Add with subfolders
        g = genpath(myrepodir);
        addpath(g)

        fprintf('Found and added to path with subfolders: %s\n', myrepodir);

    else

        fprintf('Did not find: %s\n', key_files_to_find{i});

    end

end % for

fprintf('\n%s\n', dashes);

% ------------------------------------------------------------------------
% Possibly get repositories that are missing
% ------------------------------------------------------------------------

if ~any(was_missing)
    fprintf('\nAll CANlab repositories previously installed. Saving path\n');
    savepath
else
    fprintf('\nMissing repositories:\n');
    disp(char(repo_descrip{was_missing}));

    doinstall = input(sprintf('\nTry to download from Github and install in %s?\nPress 1 for yes, 0 for no: ', main_repository_dir));

    if doinstall

        for i = 1:n_repos

            if was_missing(i)

                fprintf('\n%s\n', dashes);

                % Make sure we are in the main directory above CANlab Core
                % canlab_clone_github_repository also adds to Path
                cd(main_repository_dir)

                try

                canlab_clone_github_repository('repo', repo_names{i}, 'noverbose');

                catch

                    warning('DID NOT CLONE REPO: %s\n', repo_names{i})

                end
            end

        end

    end % doinstall

    % remove .git dirs
    strip_git_dirs;

    disp('Finished install. Saving path.');
    savepath

end % any missing

end % function
%%

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function test_ft_componentanalysis_methods

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_componentanalysis ft_rejectcomponent ft_checkdata runica fastica sobi jader dss parafac

% note 3) binica does not run on all computers, which is why it is disabled
% note 9) parafac depends on an external toolbox, which I don't have access to at the moment

% construct a simple and minimalistic raw data representation
nchan  = 10;
ntime  = 1000;
ntrial = 5;
data = [];
data.fsample = 1;
for i=1:nchan
  data.label{i} = sprintf('chan%d', i);
end
for i=1:ntrial
  data.trial{i} = sqrt(abs(randn(nchan,ntime)));
  data.time{i} = 1:ntime;
end

% run the decomposition with various methods
cfg = [];
cfg.method = 'fastica';
comp1f = ft_componentanalysis(cfg, data);
cfg.numcomponent = nchan-1;
comp1s = ft_componentanalysis(cfg, data);

cfg = [];
cfg.method = 'runica';
cfg.runica.maxsteps = 50;
comp2f = ft_componentanalysis(cfg, data);
cfg.numcomponent = nchan-1;
comp2s = ft_componentanalysis(cfg, data);

% cfg = [];
% cfg.method = 'binica';
% comp3f = ft_componentanalysis(cfg, data);
% cfg.numcomponent = nchan-1;
% comp3s = ft_componentanalysis(cfg, data);

cfg = [];
cfg.method = 'jader';
comp4f = ft_componentanalysis(cfg, data);
cfg.numcomponent = nchan-1;
comp4s = ft_componentanalysis(cfg, data);

cfg = [];
cfg.method = 'varimax';
comp5f = ft_componentanalysis(cfg, data);
cfg.numcomponent = nchan-1;
comp5s = ft_componentanalysis(cfg, data);

cfg = [];
cfg.method = 'cca';
comp6f = ft_componentanalysis(cfg, data);
cfg.numcomponent = nchan-1;
comp6s = ft_componentanalysis(cfg, data);

cfg = [];
cfg.method = 'pca';
comp7f = ft_componentanalysis(cfg, data);
cfg.numcomponent = nchan-1;
comp7s = ft_componentanalysis(cfg, data);

cfg = [];
cfg.method = 'svd';
comp8f = ft_componentanalysis(cfg, data);
cfg.numcomponent = nchan-1;
comp8s = ft_componentanalysis(cfg, data);

% cfg = [];
% cfg.method = 'parafac';
% comp9f = ft_componentanalysis(cfg, data);
% cfg.numcomponent = nchan-1;
% comp9s = ft_componentanalysis(cfg, data);

cfg = [];
cfg.method = 'dss';
comp10f = ft_componentanalysis(cfg, data);
cfg.numcomponent = nchan-1;
comp10s = ft_componentanalysis(cfg, data);

cfg = [];
cfg.method = 'sobi';
comp11f = ft_componentanalysis(cfg, data);
cfg.numcomponent = nchan-1;
comp11s = ft_componentanalysis(cfg, data);

% reproject the data to the channel level
cfg = [];
[dataout] = ft_rejectcomponent(cfg, comp1f);
[dataout] = ft_rejectcomponent(cfg, comp2f);
% [dataout] = ft_rejectcomponent(cfg, comp3f);
[dataout] = ft_rejectcomponent(cfg, comp4f);
[dataout] = ft_rejectcomponent(cfg, comp5f);
[dataout] = ft_rejectcomponent(cfg, comp6f);
[dataout] = ft_rejectcomponent(cfg, comp7f);
[dataout] = ft_rejectcomponent(cfg, comp8f);
% [dataout] = ft_rejectcomponent(cfg, comp9f);
[dataout] = ft_rejectcomponent(cfg, comp10f);
[dataout] = ft_rejectcomponent(cfg, comp11f);

cfg = [];
[dataout] = ft_rejectcomponent(cfg, comp1s);
[dataout] = ft_rejectcomponent(cfg, comp2s);
% [dataout] = ft_rejectcomponent(cfg, comp3s);
[dataout] = ft_rejectcomponent(cfg, comp4s);
[dataout] = ft_rejectcomponent(cfg, comp5s);
[dataout] = ft_rejectcomponent(cfg, comp6s);
[dataout] = ft_rejectcomponent(cfg, comp7s);
[dataout] = ft_rejectcomponent(cfg, comp8s);
% [dataout] = ft_rejectcomponent(cfg, comp9s);
[dataout] = ft_rejectcomponent(cfg, comp10s);
[dataout] = ft_rejectcomponent(cfg, comp11s);

cfg = [];
cfg.component = [1 3];
[dataout] = ft_rejectcomponent(cfg, comp1f);
[dataout] = ft_rejectcomponent(cfg, comp2f);
% [dataout] = ft_rejectcomponent(cfg, comp3f);
[dataout] = ft_rejectcomponent(cfg, comp4f);
[dataout] = ft_rejectcomponent(cfg, comp5f);
[dataout] = ft_rejectcomponent(cfg, comp6f);
[dataout] = ft_rejectcomponent(cfg, comp7f);
[dataout] = ft_rejectcomponent(cfg, comp8f);
% [dataout] = ft_rejectcomponent(cfg, comp9f);
[dataout] = ft_rejectcomponent(cfg, comp10f);
[dataout] = ft_rejectcomponent(cfg, comp11f);

cfg = [];
cfg.component = [1 3];
[dataout] = ft_rejectcomponent(cfg, comp1s);
[dataout] = ft_rejectcomponent(cfg, comp2s);
% [dataout] = ft_rejectcomponent(cfg, comp3s);
[dataout] = ft_rejectcomponent(cfg, comp4s);
[dataout] = ft_rejectcomponent(cfg, comp5s);
[dataout] = ft_rejectcomponent(cfg, comp6s);
[dataout] = ft_rejectcomponent(cfg, comp7s);
[dataout] = ft_rejectcomponent(cfg, comp8s);
% [dataout] = ft_rejectcomponent(cfg, comp9s);
[dataout] = ft_rejectcomponent(cfg, comp10s);
[dataout] = ft_rejectcomponent(cfg, comp11s);

% perform some checks on the decomposed data properties
assert(length(comp1f.label)==nchan);
assert(length(comp2f.label)==nchan);
% assert(length(comp3f.label)==nchan);
assert(length(comp4f.label)==nchan);
assert(length(comp5f.label)==nchan);
assert(length(comp6f.label)==nchan);
assert(length(comp7f.label)==nchan);
assert(length(comp8f.label)==nchan);
assert(length(comp8f.label)==nchan);
% assert(length(comp9f.label)==nchan);
assert(length(comp10f.label)==nchan);
assert(length(comp11f.label)==nchan);

% see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=903
assert(length(comp1s.label)==nchan-1);
assert(length(comp2s.label)==nchan-1);
% assert(length(comp3s.label)==nchan);
assert(length(comp4s.label)==nchan-1);
assert(length(comp5s.label)==nchan-1);
assert(length(comp6s.label)==nchan-1);
assert(length(comp7s.label)==nchan-1);
assert(length(comp8s.label)==nchan-1);
% assert(length(comp9s.label)==nchan-1);
assert(length(comp10s.label)==nchan-1);
assert(length(comp11s.label)==nchan-1);


text:
%% Copyright (C) 2014-2017, 2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @deftypemethod  @@sym {@var{c} =} coeffs (@var{p}, @var{x})
%% @deftypemethodx @@sym {@var{c} =} coeffs (@var{p})
%% @deftypemethodx @@sym {@var{c} =} coeffs (@dots{}, 'all')
%% @deftypemethodx @@sym {[@var{c}, @var{t}] =} coeffs (@var{p}, @var{x})
%% @deftypemethodx @@sym {[@var{c}, @var{t}] =} coeffs (@var{p})
%% @deftypemethodx @@sym {[@var{c}, @var{t}] =} coeffs (@dots{}, 'all')
%% Return non-zero (or all) coefficients of symbolic polynomial.
%%
%% @var{c} contains the coefficients and @var{t} the corresponding
%% terms.
%%
%% Example:
%% @example
%% @group
%% syms x
%% [c, t] = coeffs (x^6 + 3*x - 4)
%%   @result{} c = (sym) [1  3  -4]  (1×3 matrix)
%%   @result{} t = (sym 1×3 matrix)
%%        ⎡ 6      ⎤
%%        ⎣x   x  1⎦
%% @end group
%% @end example
%%
%% The polynomial can be multivariate:
%% @example
%% @group
%% syms x y
%% [c, t] = coeffs (x^2 + y*x)
%%   @result{} c = (sym) [1  1]  (1×2 matrix)
%%   @result{} t = (sym 1×2 matrix)
%%       ⎡ 2     ⎤
%%       ⎣x   x⋅y⎦
%% @end group
%%
%% @group
%% [c, t] = coeffs (x^2 + y*x, [x y])   % same
%%   @result{} c = (sym) [1  1]  (1×2 matrix)
%%   @result{} t = (sym 1×2 matrix)
%%       ⎡ 2     ⎤
%%       ⎣x   x⋅y⎦
%%
%% [c, t] = coeffs (x^2 + y*x, @{x y@})   % same
%%   @result{} c = (sym) [1  1]  (1×2 matrix)
%%   @result{} t = (sym 1×2 matrix)
%%       ⎡ 2     ⎤
%%       ⎣x   x⋅y⎦
%% @end group
%% @end example
%%
%% You can use the second argument to specify a vector or list of
%% variables:
%% @example
%% @group
%% [c, t] = coeffs (x^2 + y*x, x)
%%   @result{} c = (sym) [1  y]  (1×2 matrix)
%%   @result{} t = (sym 1×2 matrix)
%%       ⎡ 2   ⎤
%%       ⎣x   x⎦
%% @end group
%% @end example
%%
%% Omitting the second output is not recommended, especially for non-interactive
%% code, because it gives only the non-zero coefficients, and additionally
%% the output is in the ``wrong order'' compared to other polynomial-related
%% commands:
%% @example
%% @group
%% c = coeffs (x^6 + 3*x - 4)
%%   @result{} c = (sym) [-4  3  1]  (1×3 matrix)
%% @end group
%% @end example
%% @strong{Warning:} Again, note the order is reversed from the two-output
%% case; this is for compatibility with Matlab's Symbolic Math Toolbox.
%%
%% If the optional input keyword @qcode{'all'} is passed, the zero
%% coefficients are returned as well, and in the familiar order.
%% @example
%% @group
%% c = coeffs (x^6 + 3*x - 4, 'all')
%%   @result{} c = (sym) [1  0  0  0  0  3  -4]  (1×7 matrix)
%% @end group
%% @end example
%% @strong{Note:} The @qcode{'all'} feature does not yet work with
%% multivariate polynomials (https://github.com/cbm755/octsympy/issues/720).
%%
%% @seealso{@@sym/sym2poly}
%% @end deftypemethod

function [c, t] = coeffs(p, x, all)

  if (nargin == 1)
    x = [];
    all = false;
  elseif (nargin == 2)
    if (ischar (x))
      assert (strcmpi (x, 'all'), ...
              'coeffs: invalid 2nd input: if string, should be "all"')
      x = [];
      all = true;
    else
      all = false;
    end
  elseif (nargin == 3)
    assert (strcmpi (all, 'all'), ...
            'coeffs: invalid 3rd input: should be string "all"')
    all = true;
  elseif (nargin > 3)
    print_usage ();
  end

  assert (isscalar (p), 'coeffs: works for scalar input only')

  p = sym(p);

  if (isempty (x))
    x = symvar (p);
    if (isempty (x))
      x = sym('x');  % any symbol
    end
  end

  x = sym(x);

  cmd = { '(f, xx, all) = _ins'
          'if not xx.is_Matrix:'
          '    xx = sp.Matrix([xx])'
          'xx = list(xx)'
          'p = Poly.from_expr(f, *xx)'
          'if all:'
          '    terms = p.all_terms()'
          'else:'
          '    terms = p.terms()'
          'cc = [q[1] for q in terms]'
          'tt = [1]*len(terms)'
          'for i, x in enumerate(p.gens):'
          '    tt = [t*x**q[0][i] for (t, q) in zip(tt, terms)]'
          'return (Matrix([cc]), Matrix([tt]))' };

  [c, t] = pycall_sympy__ (cmd, p, x, all);

  %% SMT compat:
  % reverse the order if t is not output.
  if (nargout <= 1) && (all == false)
    c = fliplr(c);
  end

  % if nargout == 1, its simplier to use 'p.coeffs()'
end

%!error coeffs (sym(1), 2, 3, 4)
%!error <should be> coeffs (sym(1), 2, 'al')
%!error <should be> coeffs (sym(1), 'al')

%!test
%! % simple
%! syms x
%! [c, t] = coeffs(6*x*x + 27);
%! assert (isequal (c, [6 27]))
%! assert (isequal (t, [x*x 1]))

%!test
%! % specify a variable
%! syms x
%! [c, t] = coeffs(6*x*x + 27, x);
%! assert (isequal (c, [6 27]))
%! assert (isequal (t, [x*x 1]))

%!test
%! % specify another variable
%! syms x y
%! [c, t] = coeffs(6*x + 27, y);
%! assert (isequal (c, 6*x + 27))
%! assert (isequal (t, 1))

%!test
%! % weird SMT order
%! syms x
%! a1 = [27 6];
%! a2 = [6 27];
%! c = coeffs(6*x*x + 27);
%! assert (isequal (c, a1))
%! coeffs(6*x*x + 27);
%! assert (isequal (ans, a1))
%! [c, t] = coeffs(6*x*x + 27);
%! assert (isequal (c, a2))

%!test
%! % no weird order with "all"
%! syms x
%! c = coeffs(6*x*x + 27, 'all');
%! assert (isequal (c, [6 0 27]))

%!test
%! % "all"
%! syms x
%! [c, t] = coeffs(6*x*x + 27, 'all');
%! assert (isequal (c, [6 0 27]))
%! assert (isequal (t, [x^2 x 1]))

%!test
%! % "All"
%! syms x
%! [c, t] = coeffs(6*x, 'All');
%! assert (isequal (c, [6 0]))
%! assert (isequal (t, [x 1]))

%!test
%! % multivariable array
%! syms x y
%! [c, t] = coeffs(6*x*x + 27*y*x  + 36, [x y]);
%! a = [6  27  36];
%! s = [x^2  x*y  1];
%! assert (isequal (c, a))
%! assert (isequal (t, s))
%! % with list
%! [c, t] = coeffs(6*x*x + 27*y*x  + 36, {x y});
%! assert (isequal (c, a))
%! assert (isequal (t, s))

%!test
%! % other symbols treated as part of coeffs
%! syms x y
%! [c, t] = coeffs(6*x*x + 27*y*x  + 36, x);
%! a = [6  27*y  36];
%! s = [x^2  x  1];
%! assert (isequal (c, a))
%! assert (isequal (t, s))

%!error <multivariate polynomials not supported>
%! % TODO: multivariate all not working (https://github.com/cbm755/octsympy/issues/720)
%! syms x y
%! [c, t] = coeffs(6*x^2 + 7*y + 19, [x y], 'all');

%!test
%! % empty same as not specifying; maybe not SMT compatible:
%! % https://github.com/cbm755/octsympy/pull/708#discussion_r94292831
%! syms x y
%! [c, t] = coeffs(6*x*x + 27*y*x  + 36, {});
%! a = [6  27  36];
%! assert (isequal (c, a))
%! [c, t] = coeffs(6*x*x + 27*y*x  + 36);
%! assert (isequal (c, a))

%!test
%! % no input defaults to all symbols (not symvar to get x)
%! syms x y
%! [c, t] = coeffs(6*x*x + 27*y*x  + 36);
%! assert (isequal (c, [6 27 36]))

%!test
%! % non sym input
%! syms x
%! assert (isequal (coeffs(6, x), sym(6)))

%!test
%! % constant input without x
%! assert (isequal (coeffs(sym(6)), sym(6)))

%!test
%! % constant input without x
%! assert (isequal (coeffs (sym(6), {}), sym(6)))

%! % irrational coefficients
%! syms x
%! f = x^2 + sqrt(sym(2))*x;
%! [c1, t1] = coeffs (f);
%! [c2, t2] = coeffs (f, x);
%! assert (isequal (c1, c2))
%! assert (isequal (t1, t2))

repo: octsympy


text:
function y = sdpcone(varargin)
%SDPCONE Low-level operator to define several SDP constraints in vectorized
%form
%
% Input
%    X       : Linear SDPVAR object of size n^2 x N
%
% Example
%
% The typical use is when we want to define a very large number of LMI
% constraints with reduced overhead (no analysis or check for symmetry etc)
%
% This operator is very specialized and low-level, and not normally used...
%
% X = sdpvar(5);Y = sdpvar(5);
% F = sdpcone([X(:) Y(:)]); % Equivalent to [X>=0, Y>=0]
%
% See also  @SDPVAR/CONE, @SDPVAR/RCONE

y.typeflag = 57;
y=lmi(y);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Transforms a covariance matrix into a correlation matrix

function C = cov2corr(C)

S = diag(diag(C).^(-0.5));
C = S*C*S;

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function out = isequal(f, g)
%ISEQUAL   Test if two CLASSICFUN objects are equal.
%   ISEQUAL(F, G) returns logical 1 (TRUE) if the CLASSICFUN objects F and G have the
%   same domain and .onefun fields which are equal, as determined by their
%   respective class methods.

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

% Check whether .domains and .onefun fields agree:
out = (all(f.domain == g.domain) && isequal(f.onefun, g.onefun));

end


text:
function w = hermite_cubic_spline_quad_rule ( nn, xn )

%*****************************************************************************80
%
%% HERMITE_CUBIC_SPLINE_QUAD_RULE: Hermite cubic spline quadrature rule.
%
%  Discussion:
%
%    The integral is taken over the definition interval [X(1),X(NN)].
%
%    Note that if the intervals are equal in size, then the derivative
%    information in DN has no effect on the integral value,
%    except for the first and last entries.
%
%    The quadrature rule is
%
%      Integral ( XN(1) <= X <= XN(NN) ) F(X) dX is approximately
%
%      Sum ( 1 <= I <= NN ) W(1,I) * F(X(I)) + W(2,I) * F'(X(I))
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 March 2011
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Fred Fritsch, Ralph Carlson,
%    Monotone Piecewise Cubic Interpolation,
%    SIAM Journal on Numerical Analysis,
%    Volume 17, Number 2, April 1980, pages 238-246.
%
%  Parameters:
%
%    Input, integer NN, the number of data points.
%
%    Input, real XN(NN), the coordinates of the data points.
%    The entries in XN must be in strictly ascending order.
%
%    Output, real W(2,NN), the quadrature weights for F(1:NN)
%    and DN(1:NN).
%
  w = zeros ( 2, nn );

  w(1,1)      = 0.5 * ( xn(1,2)    - xn(1,1)      );
  w(1,2:nn-1) = 0.5 * ( xn(1,3:nn) - xn(1,1:nn-2) );
  w(1,nn)     = 0.5 * ( xn(1,nn)   - xn(1,nn-1)   );

  w(2,1)      =   ( xn(1,2) - xn(1,1) ).^2 / 12.0;
  w(2,2:nn-1) =   ( xn(1,3:nn) - xn(1,1:nn-2) ) ...
                .* ( xn(1,3:nn) - 2.0 * xn(1,2:nn-1) + xn(1,1:nn-2) ) / 12.0;
  w(2,nn)     = - ( xn(1,nn-1) - xn(1,nn) ).^2 / 12.0;

  return
end


text:
% gray image denoising

% @inproceedings{zhang2017learning,
%   title={Learning Deep CNN Denoiser Prior for Image Restoration},
%   author={Zhang, Kai and Zuo, Wangmeng and Gu, Shuhang and Zhang, Lei},
%   booktitle={IEEE Conference on Computer Vision and Pattern Recognition},
%   year={2017}
% }

% If you have any question, please feel free to contact with me.
% Kai Zhang (e-mail: cskaizhang@gmail.com)

% clear; clc;
addpath('utilities');
imageSets     = {'BSD68','Set12'}; %%% testing dataset
folderTest    = 'testsets';
folderModel   = 'models';
folderResult  = 'results';
taskTestCur   = 'Denoising';
if ~exist(folderResult,'file')
    mkdir(folderResult);
end

setTestCur    = imageSets{1};
imageSigmaS   = [15,25,50];
modelSigmaS   = [15,25,50];
showResult    = 1;
saveResult    = 0;
useGPU        = 1;
pauseTime     = 1;

%%% folder to store results
folderResultCur = fullfile(folderResult, [taskTestCur,'_',setTestCur]);
if ~exist(folderResultCur,'file')
    mkdir(folderResultCur);
end

%%% read images
ext         =  {'*.jpg','*.png','*.bmp'};
filePaths   =  [];
folderTestCur = fullfile(folderTest,setTestCur);
for i = 1 : length(ext)
    filePaths = cat(1,filePaths, dir(fullfile(folderTestCur,ext{i})));
end

%%% PSNR and SSIM
PSNRs = zeros(length(modelSigmaS),length(filePaths));
SSIMs = zeros(length(modelSigmaS),length(filePaths));

load(fullfile(folderModel,'modelgray.mat'));

for i = 1:length(modelSigmaS)

    disp([i,length(modelSigmaS)]);
    net = loadmodel(modelSigmaS(i),CNNdenoiser);
    net = vl_simplenn_tidy(net);
    % for i = 1:size(net.layers,2)
    %     net.layers{i}.precious = 1;
    % end
    %%% move to gpu
    if useGPU
        net = vl_simplenn_move(net, 'gpu');
    end

    for j = 1:length(filePaths)

        %%% read images
        label = imread(fullfile(folderTestCur,filePaths(j).name));
        [~,imageName,extCur] = fileparts(filePaths(j).name);
        label = im2double(label);
        randn('seed',0);
        input = single(label + imageSigmaS(i)/255*randn(size(label)));

        %%% convert to GPU
        if useGPU
            input = gpuArray(input);
        end
        res    = vl_simplenn(net,input,[],[],'conserveMemory',true,'mode','test');
        output = input - res(end).x;

        %%% convert to CPU
        if useGPU
            output = gather(output);
            input  = gather(input);
        end

        %%% calculate PSNR and SSIM
        [PSNRCur, SSIMCur] = Cal_PSNRSSIM(im2uint8(label),im2uint8(output),0,0);
        if showResult
            imshow(cat(2,im2uint8(label),im2uint8(input),im2uint8(output)));
            title([filePaths(j).name,'    ',num2str(PSNRCur,'%2.2f'),'dB','    ',num2str(SSIMCur,'%2.4f')])
            drawnow;
            if saveResult
                imwrite(im2uint8(output),fullfile(folderResultCur,[imageName,'_',num2str(imageSigmaS(i)),'_',num2str(modelSigmaS(i)),'_',num2str(PSNRCur,'%2.2f'),'.png']));
            end
            pause(pauseTime)
        end
        PSNRs(i,j) = PSNRCur;
        SSIMs(i,j) = SSIMCur;

    end
end

%%% save PSNR and SSIM metrics
save(fullfile(folderResultCur,['PSNR_',taskTestCur,'_',setTestCur,'.mat']),'PSNRs')
save(fullfile(folderResultCur,['SSIM_',taskTestCur,'_',setTestCur,'.mat']),'SSIMs')

disp([mean(PSNRs,2),mean(SSIMs,2)]);

author: cszn

repo: IRCNN

save_path: github-repos/MATLAB/cszn-IRCNN


text:
function [eta,etaf]=finishat(frac,tol,fmt)
%FINISHAT print estimated finish time of a long computation (FRAC,TOL,FMT)
% Usage:  (1)  for i=1:many
%                  finishat((i-1)/many);  % initializes on first pass when i=1
%                  ... computation ...
%              end
%
%         (2)  for i=1:many
%                  finishat([i-1 many]);  % alternative argument format
%                  ... computation ...
%              end
%
%         (3)  finishat(0);               % explicit initialization before loop
%              for i=1:many
%                  ... computation ...
%                  finishat(i/many);      % calculate fraction completed
%              end
%
%         (4)  for i=1:NI
%                  for j=1:NJ
%                      for k=1:NK
%                          finishat([i NI; j NJ; k-1 NK]); % one row per nested loop
%                          ... computation ...
%                      end
%                  end
%              end
%
% Inputs: FRAC = fraction of total comutation that has been completed
%                Alternatively at start of inner loop: [i NI; j NJ; k-1 NK ...] where i, j, k are
%                loop indices and NI, NJ, NK their limits. Use k instead of k-1 if placed at
%                the end of the inner loop. As a special case, FRAC=0 initializes the routine.
%         TOL  = Tolerance in minutes. If the estimated time has changed by less
%                than this, then nothing will be printed. [default 10% of remaining time]
%         FMT  = Format string which should include %s for estimated finish time, %d for remaining minutes and %f for fraction complete
%
% Output: ETA  = string containing the expected finish time
%                specifying this will suppress printing message to std err (fid=2)
%         ETAF = expected finish time as a daynumber
%
% Example:       finishat(0);
%                for i=1:many
%                    long computation;
%                    finishat(i/many);
%                end

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

persistent oldt oldnw
if nargin<3
    fmt='Estimated finish at %s (%.2f done, %d min remaining)\n';

end
nf=size(frac,1);
if all(frac(:,1)<=[ones(nf-1,1); 0]) % initialize if fraction done is <=0
    oldt=0;
    eta='Unknown';
    tic;
else
    if size(frac,2)==2
        fp=cumprod(frac(:,2));
        frac=sum((frac(:,1)-1)./fp)+1/fp(end);
    end
    nw=now;                             % current time as serial number
    sectogo=(1/frac-1)*toc;      % seconds to go
    newt=nw+sectogo/86400;       % add estimated time in days
    if nargin<2 || ~numel(tol)
        tol=max(0.1*(newt-nw)*1440,1);
    end
    if ~exist('oldt','var') || oldt==0 || (abs(newt-oldt)>tol/1440 && (nw-oldnw)>10/86400) || (nw-oldnw)>10/1440 || nargout>0
        oldt=newt;
        if floor(oldt)==floor(nw)
            df='HH:MM';
        else
            df='HH:MM dd-mmm-yyyy';
        end
        eta=datestr(oldt,df);
        if ~nargout
            ix=find(fmt=='%',1);
            while ~isempty(ix)
                fprintf(2,fmt(1:ix-1));
                fmt=fmt(ix:end);
                ix=find(fmt>='a' & fmt<='z',1); % find letter
                switch fmt(ix)
                    case 's'
                        fprintf(2,fmt(1:ix),eta);
                    case 'd'
                        fprintf(2,fmt(1:ix),round(sectogo/60));
                    case 'f'
                        fprintf(2,fmt(1:ix),frac);
                end
                fmt=fmt(ix+1:end);
                ix=find(fmt=='%',1);
            end
            fprintf(2,fmt);
        end
        oldnw=nw;                           %
    end
end
etaf=oldt;

author: jtkim-kaist

repo: Speech-enhancement


text:
function [X,cost,test,stats] = completion_rankincrease_adaptive_new2( method, A_Omega, Omega, A_Omega_C, Omega_C, A_Gamma, Gamma, X0, opts )

    if ~isfield( opts, 'maxrank');      opts.maxrank = 4  ;         end
    if ~isfield( opts, 'cg');           opts.cg = true;             end
    if ~isfield( opts, 'tol');          opts.tol = 1e-6;            end
    if ~isfield( opts, 'reltol');       opts.reltol = 1e-8;         end
    if ~isfield( opts, 'reltol_final'); opts.reltol_final = eps;    end
    if ~isfield( opts, 'maxiter');      opts.maxiter = 10;          end
    if ~isfield( opts, 'maxiter_final');opts.maxiter_final = 20;    end
    if ~isfield( opts, 'locked_tol');   opts.locked_tol = 1;        end
    if ~isfield( opts, 'epsilon');      opts.epsilon = 1e-8;        end
    if ~isfield( opts, 'verbose');      opts.verbose = false;       end

    if strcmpi( method, 'GeomCG' )
        completion = @( A_Omega, Omega, A_Gamma, Gamma, X0, opts ) ...
                            completion_orth( A_Omega, Omega, A_Gamma, Gamma, X0, opts )
    elseif strcmpi( method, 'ALS' )
        completion = @( A_Omega, Omega, A_Gamma, Gamma, X0, opts ) ...
                            completion_als( A_Omega, Omega, A_Gamma, Gamma, X0, opts )
    end
    d = X0.order;

    test = [];
    control_old = inf;

    % ===========================================
    disp('____________________________________________________________________');
    disp(['Completion with with starting rank r = [ ' num2str(X0.rank) ' ] ...']);
    [X,cost,control,stats] = completion( A_Omega, Omega, A_Gamma, Gamma, X0, opts);

    stats.rankidx = [length(cost)];

    disp('____________________________________________________________________');
    disp(['Increasing rank ... ']);

    locked = zeros(1,d+1);

    for k = 2:opts.maxrank
        for i = 2:d

            disp(['Locked cores:' num2str(locked) ])
            if locked(i)
                disp(['Rank r(' num2str(i) ') is locked. Skipping.']);
            else
                r = X.rank;
                disp(['Trying to increase rank r(' num2str(i) ') from ' num2str(r(i)) ' to ' num2str(r(i)+1) ':']);
                Xnew = increaseRank(X, 1, i, opts.epsilon);
                Xnew = orthogonalize(Xnew, d);
                if i==d && k == opts.maxrank
                    opts.maxiter = opts.maxiter_final;
                end
                [Xnew,cost_tmp,control_tmp,stats_tmp] = completion( A_Omega, Omega, A_Omega_C, Omega_C, Xnew, opts);
                stats.rankidx = [stats.rankidx, length(cost_tmp)];
                disp( ['Current cost function:            ', num2str(cost_tmp(end)) ]);

                progress = (control_tmp(end) - control_old )/control_old;
                disp( ['Current rel. progress on control: ' num2str(progress)]);

                if  progress > opts.locked_tol
                    disp(['     ... failed. Reverting.']);
                else
                   disp(['     ... accepted.']);
                   X = Xnew;
                   control_old = control_tmp(end)
                   test_current = norm(X(Gamma) - A_Gamma)/ norm(A_Gamma)
                   disp( ['Current error on test set Gamma:  ', num2str(test_current) ]);
                   test = [test, test_current];
                end

                if ~isempty(stats.time)
                    stats_tmp.time = stats_tmp.time + stats.time(end);
                end

                cost = [cost; cost_tmp];
                control = [control; control_tmp];
                stats.time = [stats.time, stats_tmp.time];

            end

        end
    end

end

author: NicolasBoumal

repo: manopt


text:
function [pindx, pval] = peakdetect3(dat, threshold, mindist)

% PEAKDETECT3 detects peaks above a certain threshold in single-channel data
%
% Use as
%   [pindx, pval] = peakdetect3(dat, threshold, mindist)
%
% See also PEAKDETECT, PEAKDETECT2

% Copyright (C) 2000-2005, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% threshold the data
tr = dat>threshold;

% the derivative of the data changes its sign at the peak
td = diff(dat);
td = td(1:(end-1))>0 & td(2:end)<0;
td = [0 td 0];

pindx = find(td & tr);

if nargin>2 && length(pindx)>0
  % find the peaks that are too close to each other
  pd = [inf diff(pindx)];
  pindx = pindx(pd>mindist);
end

if nargout>1
  pval = dat(pindx);
end

author: fieldtrip

repo: fieldtrip


text:
function c = cond_2x2(A)
%COND   Condition number with respect to inversion.
%   COND2d(X) returns the 2-norm condition number (the ratio of the
%   largest singular value of X to the smallest).  Large condition
%   numbers indicate a nearly singular matrix.
%
%   COND(X,P) returns the condition number of X in 2-norm:
%
%      NORM(X,2) * NORM(INV(X),2).
%
%   Class support for input X:
%      float: double, single
%
%   See also RCOND, CONDEST, CONDEIG, NORM, NORMEST.

%   Copied from MATLAB's cond function.

assert(~issparse(A));
[d1,d2,n] = size(A);
assert(d1==2 && d2==2);
A = reshape(A,[4,n]);

% A2 = A'*A
tmp = A(1,:).*A(3,:) + A(2,:).*A(4,:);
A2 = cat(1,A(1,:).^2 + A(2,:).^2,...
  tmp,tmp,...
  A(3,:).^2 + A(4,:).^2);

s = sqrt(eigs_2x2(A2));
c = zeros([1,n],class(A));

issing = any(s==0,1);
c(issing) = inf;
c(~issing) = max(s(:,~issing),[],1)./min(s(:,~issing),[],1);

author: kristinbranson

repo: JAABA


text:
function [pv01, fx]=pitchestm(data, fs, nfr10, pv01)
% function [pv01, pvblk, pvblkb, pv]=pitchestm(data, fs, nfr10, pv01)

[fx,tt,pv,fv]=fxpefac(data, fs); % should plus 3
npv=length(pv);
pv01=zeros(nfr10,1); sign_pv=0;
for i=1:npv
    if pv(i)>0.25
        pv01(i+3) =1;
        if sign_pv==0
            sign_pv=1;
            nstart=i;
        end
    else
        if sign_pv==1
            sign_pv=0;
            nstop=i-1;
            if nstop-nstart<3
                pv01(nstart+3:nstop+3)=0;  % Remove 2 frames only pitch
            end
        end
    end
end
pv01(1:3)=pv01(4);
fxtmp(1:3)=fx(4); fxtmp(4:npv+3)=fx(1:npv);
if (npv+3) < nfr10
    pv01(npv+4:nfr10)=pv01(npv+3);
    fxtmp(npv+4:nfr10)=fx(npv);
else
    pv01=pv01(1:nfr10);
    fxtmp=fxtmp(1:nfr10);
end
fx=fxtmp;

author: zhenghuatan

repo: rVAD

save_path: github-repos/MATLAB/zhenghuatan-rVAD


text:
function [cnt1, cnt2] = elec1020_follow(pnt, dhk, v1, v2, v3, feedback)

% ELEC1020_FOLLOW

% Copyright (C) 2003, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.

if nargin<6
  feedback = false;
end

% compute the distribution of edge lengths
edge = [
  pnt(dhk(:,1),:) - pnt(dhk(:,2),:)
  pnt(dhk(:,2),:) - pnt(dhk(:,3),:)
  pnt(dhk(:,3),:) - pnt(dhk(:,1),:)
  ];
edgelength = sqrt(sum(edge.^2,2));

% the tolerance depends on the edge length
tolerance       = min(edgelength)*1e-6;
tolerance_limit = min(edgelength)*1e-3;

npnt = size(pnt,1);
ndhk = size(dhk,1);

pside = nan(1,npnt);
for i=1:npnt
  % determine on which side of the plane each vertex lies
  pside(i) = ptriside(v1, v2, v3, pnt(i,:), tolerance);
end

% dcut = zeros(ndhk,1);
% for i=1:ndhk
%   % find the triangles that are intersected by the plane
%   if sum(pside(dhk(i,:))==0)==2
%     dcut(i) = 1;
%   elseif sum(pside(dhk(i,:))==0)==1 & sum(pside(dhk(i,:))==1)==1 & sum(pside(dhk(i,:))==-1)==1
%     dcut(i) = 1;
%   elseif sum(pside(dhk(i,:))==1)==2 & sum(pside(dhk(i,:))==-1)==1
%     dcut(i) = 1;
%   elseif sum(pside(dhk(i,:))==1)==1 & sum(pside(dhk(i,:))==-1)==2
%     dcut(i) = 1;
%   end
% end
tmp  = pside(dhk);
dcut = true(ndhk,1);
dcut(all(tmp== 1,2)) = false;
dcut(all(tmp==-1,2)) = false;

% continue working with only the intersecting triangles
dhk = dhk(dcut,:);
ndhk = size(dhk,1);

% for each triangle determine the neighbouring triangles
neighb = zeros(ndhk,ndhk);
for i=1:ndhk
  for j=(i+1):ndhk
    if dhk(i,1)==dhk(j,1)
      neighb(i,j) = 1;
    elseif dhk(i,1)==dhk(j,2)
      neighb(i,j) = 1;
    elseif dhk(i,1)==dhk(j,3)
      neighb(i,j) = 1;
    elseif dhk(i,2)==dhk(j,1)
      neighb(i,j) = 1;
    elseif dhk(i,2)==dhk(j,2)
      neighb(i,j) = 1;
    elseif dhk(i,2)==dhk(j,3)
      neighb(i,j) = 1;
    elseif dhk(i,3)==dhk(j,1)
      neighb(i,j) = 1;
    elseif dhk(i,3)==dhk(j,2)
      neighb(i,j) = 1;
    elseif dhk(i,3)==dhk(j,3)
      neighb(i,j) = 1;
    else
      neighb(i,j) = 0;
    end
    neighb(j,i) = neighb(i,j);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% make a contour from v1 via v2 to v3

% find the nearest triangle on which point v1 projects
v1_dist = inf;
for i=1:ndhk
  % shift a fraction towards v2 to avoid starting in the vertex of a wrong triangle
  [proj, dist] = ptriproj(pnt(dhk(i,1),:), pnt(dhk(i,2),:), pnt(dhk(i,3),:), v1+tolerance*(v2-v1)/norm(v2-v1), 1);
  if dist<v1_dist
    v1_dist = dist;
    v1_indx = i;
    v1_proj = ptriproj(pnt(dhk(i,1),:), pnt(dhk(i,2),:), pnt(dhk(i,3),:), v1, 1);
  end
end

% find the nearest triangle on which point v3 projects
v3_dist = inf;
for i=1:ndhk
  [proj, dist] = ptriproj(pnt(dhk(i,1),:), pnt(dhk(i,2),:), pnt(dhk(i,3),:), v3, 1);
  if dist<v3_dist
    v3_dist = dist;
    v3_indx = i;
    v3_proj = proj;
  end
end

% intersect the triangle containing the projection of v1 with the plane
[l1, l2] = tritrisect(v1, v2, v3, pnt(dhk(v1_indx,1),:), pnt(dhk(v1_indx,2),:), pnt(dhk(v1_indx,3),:));
if pntdist(l1,v2) < pntdist(l2,v2)
  % remember the projection point and l1 as the first line segment
  cnt1 = [v1_proj];
  cnt2 = [l1];
  prev_proj = l1;
  prev_indx = v1_indx;
else
  % remember the projection point and l2 as the first line segment
  cnt1 = [v1_proj];
  cnt2 = [l2];
  prev_proj = l2;
  prev_indx = v1_indx;
end

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

while(1)

  for i=find(neighb(prev_indx,:))

    ncnt = size(cnt1,1);
    c1 = pnt(dhk(i,1),:);
    c2 = pnt(dhk(i,2),:);
    c3 = pnt(dhk(i,3),:);
%
%     if feedback
%       ft_plot_mesh(struct('pos', pnt, 'tri', dhk(i,:)), 'edgecolor', 'b', 'facecolor', 'none')
%       ft_plot_mesh(c1, 'vertexsize', 30, 'vertexcolor', 'r')
%       ft_plot_mesh(c2, 'vertexsize', 30, 'vertexcolor', 'r')
%       ft_plot_mesh(c3, 'vertexsize', 30, 'vertexcolor', 'r')
%       ft_plot_mesh(prev_proj, 'vertexsize', 30, 'vertexcolor', 'g')
%       drawnow
%     end

    [proj, dist] = ptriproj(c1, c2, c3, prev_proj, 1);

    if dist<tolerance
      [l1, l2] = tritrisect(v1, v2, v3, c1, c2, c3);

      if pntdist(l1, cnt1(ncnt,:))<tolerance && pntdist(l2, cnt2(ncnt,:))<tolerance
        continue
      elseif pntdist(l1, cnt2(ncnt,:))<tolerance && pntdist(l2, cnt1(ncnt,:))<tolerance
        continue
      end

      if pntdist(l1, prev_proj) < pntdist(l2, prev_proj)
        cnt1 = [cnt1; l1];
        cnt2 = [cnt2; l2];
        prev_proj = l2;
        prev_indx = i;
        break
      else
        cnt1 = [cnt1; l2];
        cnt2 = [cnt2; l1];
        prev_proj = l1;
        prev_indx = i;
        break
      end
    end

  end

  % stop if no new segment was added
  if ~(size(cnt1,1)>ncnt)
    tolerance = 2*tolerance;
    if tolerance>=tolerance_limit
      ft_warning('premature end of contour')
      break
    else
      ft_warning('increasing tolerance');
    end
  end

  % stop if we arrive on the triangle with the endpoint
  if pntdist(prev_proj, v3_proj)<tolerance
    % replace the current endpoint with the projection of v3
    cnt2(size(cnt2,1),:) = v3_proj;
    break
  end

  % stop if we arrive on the triangle with the endpoint
  if prev_indx==v3_indx
    % replace the current endpoint with the projection of v3
    cnt2(size(cnt2,1),:) = v3_proj;
    break
  end

end % while

if feedback
  X = [cnt1(:,1) cnt2(:,1)];
  Y = [cnt1(:,2) cnt2(:,2)];
  Z = [cnt1(:,3) cnt2(:,3)];
  line(X, Y, Z, 'color', 'k')
  drawnow
end

end % function