text:
function [expdirlistfile,expdirs,expdirsCreated,expdirs_linux] = setUpDirUsingCluster(rootdirs,varargin)

if nargin < 1,
  rootdirs = {};
end

expdirs = {};
expdirsCreated = {};
expdirs_linux = {};

% ignore features arguments -- these will be computed on the cluster
[expdirlistfile,doforce,linuxroot,args] = myparse_nocheck(varargin,'expdirlistfile','','doforce',false,...
  'linuxroot','/groups/hantman/hantmanlab');

% file to output directories to process to
if isempty(expdirlistfile),
  expdirlistfile = sprintf('JAABADirsForCluster%s.txt',datestr(now,'yyyymmddTHHMMSS'));
end

argisremove = find(strcmpi(args(1:2:end-1),'features'));
if ~isempty(argisremove),
  args([argisremove,argisremove+1]) = [];
end
args = [args,'features',false,'doforce',doforce];

% select root directories if not given
if isempty(rootdirs),
  p = fileparts(mfilename('fullpath'));
  rcfilename = fullfile(p,'.setUpDir.mat');
  startdir = '.';
  if exist(rcfilename,'file'),
    try %#ok<TRYNC>
      rc = load(rcfilename);
      if exist(rc.startdir,'dir'),
        startdir = rc.startdir;
      end
    end
  end
  rootdirs = uipickfiles('FilterSpec',startdir,'DirsOnly',true,'Output','cell');
  if isnumeric(rootdirs) || isempty(rootdirs),
    fprintf('No root directories selected, returning\n');
    return;
  end
end

% if only one root directory is given
if ~iscell(rootdirs),
  rootdirs = {rootdirs};
end

if ~isempty(rootdirs),
  p = fileparts(mfilename('fullpath'));
  rcfilename = fullfile(p,'.setUpDir.mat');
  p = fileparts(rootdirs{end});
  if exist(p,'dir'),
    startdir = p; %#ok<NASGU>
    try %#ok<TRYNC>
      save(rcfilename,'startdir');
    end
  end
end

% run the annotation and directory creating steps
expdirsCreated = {};
for i = 1:numel(rootdirs),
  expdirs1 = setUpDir(rootdirs{i},args{:});
  if isempty(expdirs1),
    fprintf('%d: No directories created within root directory %s.\n',i,rootdirs{i});
  else
    fprintf('%d: Created the following directories within root directory %s:\n',i,rootdirs{i});
    fprintf('  %s\n',expdirs1{:});
    expdirsCreated = [expdirsCreated,expdirs1(:)']; %#ok<AGROW>
  end
end

% get a list of experiment directories to compute features for
expdirs = {};
for i = 1:numel(rootdirs),
  expdirs1 = genAllFeaturesExpDirList(rootdirs{i},'doforce',doforce,args{:});
  if isempty(expdirs1),
    fprintf('%d: No directories within root directory %s to compute features for.\n',i,rootdirs{i});
  else
    fprintf('%d: The following directories within root directory %s need processing:\n',i,rootdirs{i});
    fprintf('  %s\n',expdirs1{:});
    expdirs = [expdirs,expdirs1]; %#ok<AGROW>
  end
end

if isempty(expdirs),
  fprintf('No directories to process.\n');
  return;
end

expdirs_linux = {};
if ispc,
  for i = 1:numel(expdirs),
    % strip off the <driveletter>:\
    m = regexp(expdirs{i},'^[A-Z]:\\(.*)$','tokens','once');
    if isempty(m),
      fprintf('JAABA directory %s does not start with <driveletter>:\ -- not sure how to make this compatible with the cluster, skipping.\n',expdirs{i});
      continue;
    end
    % change file separator
    m = strrep(m{1},'\','/');
    expdirs_linux{end+1} = [linuxroot,'/',m]; %#ok<AGROW>
  end
elseif ismac,
  % file share is in the Mac path
  [~,fileshare] = fileparts(linuxroot);
  for i = 1:numel(expdirs),
    % strip off the /Volumes part
    m = regexp(expdirs{i},'^/Volumes/([^/]*)/(.*)','tokens','once');
    if isempty(m),
      fprintf('JAABA directory %s could not be parsed as /Volumes/(something)/(rest) -- not sure how to make this compatible with the cluster, skipping.\n',expdirs{i});
      continue;
    end
    if ~strcmp(fileshare,m{1}),
      fprintf('JAABA directory %s does not match specified file share %s -- not sure how to make this compatible with the cluster, skipping.\n',expdirs{i},fileshare);
      continue;
    end
    expdirs_linux{end+1} = fullfile(linuxroot,m{2}); %#ok<AGROW>
  end
else
  expdirs_linux = cell(size(expdirs));
  for i = 1:numel(expdirs),
    expdirs_linux{i} = absolutifyFileName(expdirs{i},pwd);
  end
end

fprintf('Outputting list of JAABA directories to process on the cluster to file %s...\n',expdirlistfile);
fid = fopen(expdirlistfile,'w');
if fid < 0,
  error('Could not write to file %s',expdirlistfile);
end

fprintf(fid,'%s\n',expdirs_linux{:});
fclose(fid);

fprintf('\nOutput %d JAABA experiments to file %s to be processed on the cluster.\n\n',numel(expdirs_linux),expdirlistfile);

if isunix && ~ismac,
  expdirlistfile_abs = absolutifyFileName(expdirlistfile,pwd);
  fprintf('To process on the cluster, do the following:\n');
  fprintf('1. Log into login2 (this is done with the terminal command "ssh login2")\n\n');
  fprintf('2. On login2, run the command:\n');
  fprintf('    /groups/branson/bransonlab/share/qsub_genAllFeatures.pl %s\n\n',expdirlistfile_abs);
  fprintf('To check on the progress of your jobs, you can type "qstat" at the command prompt to get a list of jobs running.\n');
else
  fprintf('To process on the cluster, do the following:\n');
  fprintf('1. Put the file created by this function, %s, on the Janelia file system, e.g. in the hantmanlab file share.\n',expdirlistfile);
  fprintf('2. Log into login2 (this is done with the terminal command "ssh login2")\n\n');
  fprintf('3. On login2, run the command:\n');
  fprintf('    /groups/branson/bransonlab/share/qsub_genAllFeatures.pl <JAABADirectoryListFile>\n');
  [~,n,e] = fileparts(expdirlistfile);
  n = [n,e];
  fprintf('where <JAABADirectoryListFile> is the *full Linux path* to the location of the file %s on the Janelia file server,\n',n);
  fprintf('Example: /groups/branson/bransonlab/share/qsub_GenAllFeatures.pl /groups/hantman/hantmanlab/JAABADirLists/%s\n\n',n);
  fprintf('To check on the progress of your jobs, you can type "qstat" at the command prompt to get a list of jobs running.\n');
end

author: kristinbranson

repo: JAABA


text:
function dfdupfigure(dffig)
%DFDUPFIGURE Make a duplicate, editable copy of the distribution fitting figure

%   $Revision: 1.1.6.5 $  $Date: 2004/03/26 13:30:58 $
%   Copyright 2003-2004 The MathWorks, Inc.

% Copy the regular axes, not the legend axes
f = figure;
ax = findall(dffig,'Type','axes','Tag','main');
copyobj(ax,f);
newax = findall(f,'Type','axes','Tag','main');

% Adjust layout in new figure, but don't add axis controls
dfadjustlayout(f,'off');

for i=1:length(ax)
   % Remove any context menus and callbacks associated with the old figure
   set(findall(newax(i),'Type','line'),...
       'DeleteFcn','','UIContextMenu',[],'ButtonDownFcn','');

   % Make a new legend based on the original, if any
   [legh,unused,h0,txt] = legend(ax(i));
   if length(h0)>0
      c0 = get(ax(i),'Child');
      c1 = get(newax(i),'Child');
      h1 = h0;
      for j=length(h0):-1:1
         k = find(c0==h0(j));
         if isempty(k)
            h1(j) = [];
            txt(j) = [];
         else
            % Convert to lineseries
            h1(j) = hgline2lineseries(c1(k(1)));
         end
      end

      % It's hard to match the old legend position adequately because the
      % two figure sizes may be very different.  Since the original legend
      % positions were TL or TR, we'll just pick whichever of these is
      % closest to the current position.  If the user moved the legend in
      % the original figure, he or she may need to do the same here.
      oldpos = get(legh,'Position');
      if oldpos(1)<.4
         newpos = 'NW';
      else
         newpos = 'NE';
      end
      legend(newax(i),h1,txt,'Location', newpos);
   end
end

author: zouchuhang

repo: LayoutNet


text:
function [data, labels, imind] = mcmcGenerateRandomSegments3(imsegs, imdir, ...
            adjlist, spdata, edata, vclassifierSP, hclassifierSP, eclassifier)
% Generate random segments for training good (1) vs. bad (-1) segments
% segment, where good segments consist entirely of one label

nimages = numel(imsegs);
nfeatures = 36;

nsegments = 60000;

disp(num2str(nimages))

data = zeros(nsegments, nfeatures);

labels = zeros(nsegments, 1);

imind = zeros(nsegments, 1);

count = 0;

for f = 1:nimages

    disp([num2str(f) ': ' imsegs(f).imname])

    im = im2double(imread([imdir '/' imsegs(f).imname]));
    %imdata = mcmcComputeImageData(im, imsegs(f)) ;

    segimage = imsegs(f).segimage;

    %features{f} = mcmcGetSegmentFeatures(imsegs(f), spfeatures{f}, imdata,
    %gtmaps{f}, [1:max(gtmaps{f})]);

    [pvSP, phSP, pE, smap] = mcmcInitialize(spdata{f}, edata{f}, ...
        adjlist{f}, imsegs(f), vclassifierSP, hclassifierSP, eclassifier, 'labels');

    origmap = smap;

    while count < round(nsegments*(f/nimages))

        %figure(1), hold off, imshow(label2rgb(origmap(segimage)))
        [smap, newseg, origseg, segadjmat] =  ...
            mcmcGenerateProposals(pE, adjlist{f}, smap, im, segimage);

        % display segments
        %figure(2), hold off, imagesc(label2rgb(smap(segimage)));
        %figure(1), displaySegmentGraph(im, smap(segimage), segadjmat);

        neighborseg = [newseg find(segadjmat(newseg, :))];

        % keep all potential segments as data points
        sind = find(smap==newseg);
        for ni = neighborseg
            smap2 = smap;
            smap2(sind) = ni;
            if ni~=newseg
                sind2 = find(smap2 > newseg);
                smap2(sind2) = smap2(sind2) - 1;
                if ni > newseg
                    ni = ni - 1;
                end
            end
            count = count + 1;

            data(count, :) = mcmcGetSegmentationFeatures(pvSP, phSP, pE, adjlist{f}, imsegs(f).npixels, smap2, ni);
            %data(count, :) = mcmcGetSegmentFeatures(imsegs(f), spdata{f}, imdata, smap2, ni);

            labels(count) = getMixUniLabel(imsegs(f), smap2, ni); % 0 for mix, 1 for uni

            imind(count) = f;

            if count == nsegments
                break;
            end

        end

        % remove possibility of returning to previous state
        if numel(neighborseg)>1
            neighborseg(find(neighborseg==origseg)) = [];
        end

        % randomly select neighbor
        nn = numel(neighborseg);
        goodind = find(labels((count-nn+1):count)); % good neighbors
        w = 0.25*ones(nn, 1);
        w(goodind) = 0.5; % so that a uniform segment is twice as likely to be picked
        w = cumsum(w / sum(w));
        ri = find(rand(1) < w);

        ni = neighborseg(ri(1)); % get neighbor segment

        % attach segment newseg to segment ni (unless newseg==ni)
        smap(sind) = ni;
        if ni~=newseg
            sind2 = find(smap > newseg);
            smap(sind2) = smap(sind2) - 1;
        end
        %figure(3), hist(smap, [1:max(smap)])

%         mask = zeros([size(im, 1) size(im, 2)]);
%         mask(find(smap(segimage)==ni)) = 1;
%         disp(num2str([max(smap) labels(count-nn+ri(1))]))
%         figure(2), hold off, imagesc(im .* repmat(mask, [1 1 3])), axis image
    end

    disp(num2str(mean(labels(1:count)==1)))

end

% remove segments that are not clearly good or bad
ind = find(labels==0);
labels(ind) = [];
data(ind, :) = [];
imind(ind) = [];
% end mcmcGenerateRandomSegments3

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

function label = getMixUniLabel(imsegs, smap, si)
% returns -1 for mixed; 0 for neither; 1 for uniform

lcount = zeros(numel(imsegs.label_names), 1);
sind = find(smap==si);

for k = sind'
    slab = imsegs.labels(k);
    if slab > 0
        lcount(slab) = lcount(slab) + imsegs.npixels(k);
    end
end
npix = sum(lcount);

if npix > 0
    lcount = lcount / npix;
else
    label = 0;
    return;
end

label = 0;
if max(lcount) < 0.95 || ((1-max(lcount))*npix > 500)
    label = -1;
elseif max(lcount) > 0.99
    label = 1;
end


text:
%mg_diff   GMG preconditioner for diffusion problem
%IFISS scriptfile: AR, HCE, DJS; 15 April 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

nc=log2(length(y)-1);
%
% compute new MG data or reload existing data?
compute_mg = default('compute / load MG data? 1/2 (default 1)',1);
if compute_mg==2
   load mgdata_diff
else
   h=2^(1-nc);
   fprintf('Setting up MG data ...')
   % top level
   mgdata(nc).matrix=Agal;
   if domain==1 % square
      mgdata(nc).prolong=mg_prolong(2^nc,2^nc,x,y);
   elseif domain==2 % L-shaped
      x=x';y=y';
      mgdata(nc).prolong=mg_prolong_ell(nc,x,y);
   elseif domain==3 % step
      x=x';y=y';
      mgdata(nc).prolong=mg_prolong_step(nc,x,y);
   else
      error('Multigrid solver not implemented for this domain')
   end
   xn=x;yn=y;
   % loop over remaining levels
   for level = nc-1:-1:2;
      xn=xn(1:2:end);yn=yn(1:2:end);
      if domain==1 % square
         mgdata(level).matrix=mg_diff_setup(xn,yn);
         mgdata(level).prolong=mg_prolong(2^level,2^level,xn,yn);
      elseif domain==2
         mgdata(level).matrix=mg_diff_setup_ell(xn,yn);
         mgdata(level).prolong=mg_prolong_ell(level,xn,yn);
      else
         mgdata(level).matrix=mg_diff_setup_step(xn,yn);
         mgdata(level).prolong=mg_prolong_step(level,xn,yn);
      end
   end
   fprintf('done\n')
   gohome, cd datafiles, save mgdata_diff mgdata
   if domain~=1, x=x';y=y';, end
end
%
% MG parameters
smooth = default('Jacobi / ILU smoother? 1/3 (default is ILU)',3);
if smooth==3
   % point ILU
   sweeps=1;stype=1;
elseif smooth==2
   % Gauss-Seidel
   stype = default('point / line Gauss-Seidel? 1/2 (default is line)',2);
   if stype==2
      sweeps = default('number of Gauss-Seidel sweeps? 1/2/3/4 (default is 2)',2);
   else
      sweeps=1;
   end
else
   % point Jacobi
   sweeps=1;stype=1;
end
npre = default('number of pre-smoothing steps? (default is 1)',1);
npost = default('number of post-smoothing steps? (default is 1)',1);
%
% construct smoother
if domain==1
   smooth_data = mg_smooth(mgdata,nc,sweeps,smooth,stype);
elseif domain==2
   smooth_data = mg_smooth_ell(mgdata,nc,sweeps,smooth,stype);
else
   smooth_data = mg_smooth_step(mgdata,nc,sweeps,smooth,stype);
end

author: johannesgerer

repo: jburkardt-m


function texture = faceTexture( FV,R,t,s,im)
%FACETEXTURE Summary of this function goes here
%   Detailed explanation goes here

rotpts = R*FV.vertices';
pts2D =[s.*(rotpts(1,:)+t(1)); s.*(rotpts(2,:)+t(2))];
[rows,cols]=meshgrid(1:size(im,2),1:size(im,1));

FV.facevertexcdata = zeros(0);

for col=1:3%size(im,1)+1-
FV.facevertexcdata(:,col) = interp2(rows,cols,im(:,:,col),pts2D(1,:),size(im,1)+1-pts2D(2,:));
end

%non visible texture set to zero
nonVisibleVertices = setdiff(1:length(FV.vertices), visiblevertices(FV, R));
for col=1:3
FV.facevertexcdata(nonVisibleVertices,col) = NaN;
end

texture = FV.facevertexcdata;
end


text:
function spm_dartel_dotprods(job)
% Generate a kernel from dot-products of images
% FORMAT spm_dartel_dotprods(job)
% job.images  - Images to use
% job.dotprod - Part of filename for results
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_dartel_dotprods.m 4492 2011-09-16 12:11:09Z guillaume $

P      = strvcat(job.images);
[pth,nam,ext] = fileparts(job.dotprod);
ofname = fullfile(pwd,['dp_' nam '.mat']);

N = nifti(P);
n = numel(N);
dm= size(N(1).dat);
dat=cell(1,numel(N));
for i=1:numel(N),
    dat{i} = reshape(N(i).dat,[prod(dm),1]);
end
Phi = zeros(n,n);

if isfield(job,'weight') && ~isempty(job.weight),
    Pmsk = strvcat(job.weight);
    Nmsk = nifti(Pmsk);
    msk  = Nmsk.dat;
    dmsk = size(msk);
    if any(dmsk(1:3) ~= dm(1:3)),
        error('Wrong sized weighting image.');
    end
    msk = reshape(msk,[prod(dmsk),1]);
    if numel(dmsk)==3,
        msk1 = msk;
        for i=2:prod(dm(4:end)),
            msk = [msk;msk1];
        end
    end
end

mem = 32*1024*1024;  % Mbytes of RAM to use
bs  = ceil(mem/8/n); % Block size
nd  = prod(dm);
nblock = ceil(prod(dm)/bs);
spm_progress_bar('Init',nblock,...
                 'Generating kernel','Blocks complete');
for k=1:nblock,
    o = bs*(k-1)+(1:bs);
    o = o(o<nd);
    if exist('msk','var'),
        wt  = msk(o);
        tmp = wt>0;
        o   = o(tmp);
        wt  = wt(tmp);
    end
    if ~isempty(o),
        X = zeros(numel(o),numel(dat));
        for i=1:n,
            tmp    = dat{i}(o);
            tmp(~isfinite(tmp)) = 0;
            if exist('wt','var'), tmp = tmp.*wt; end
            X(:,i) = tmp;
        end
        Phi    = Phi + X'*X;
        clear X
    end
    spm_progress_bar('Set',k);
end
spm_progress_bar('Clear');
input = job;
typ   = 'images';
save(ofname,'Phi','input','typ', spm_get_defaults('mat.format'));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function geometry_test049 ( )

%*****************************************************************************80
%
%% TEST049 tests PARALLELOGRAM_CONTAINS_POINT_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  ntest = 5;
%
%  In, Out, Out, Out, Out
%
  ptest(1:3,1:ntest) = [ ...
    1.0,  1.0,  0.5; ...
    3.0,  3.0,  0.0; ...
    0.5,  0.5, -0.1; ...
    0.1,  0.1,  0.5; ...
    1.5,  1.6,  0.5 ]';

  p1(1:3,1) = [ 0.0; 0.0; 0.0 ];
  p2(1:3,1) = [ 2.0; 2.0; 0.0 ];
  p3(1:3,1) = [ 1.0; 1.0; 1.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST049\n' );
  fprintf ( 1, '  PARALLELOGRAM_CONTAINS_POINT_3D determines if a point\n' );
  fprintf ( 1, '    is within a parallelogram in 3D.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '           P           Inside?\n' );
  fprintf ( 1, '\n' );

  for j = 1 : ntest

    p(1:3,1) = ptest(1:3,j);

    inside = parallelogram_contains_point_3d ( p1, p2, p3, p );

    fprintf ( 1, '  %12f  %12f  %12f  %d\n', p(1:3,1), inside );

  end

  return
end


text:
function z = trace_inv( Y )

% TRACE_INV   Trace of the inverse of a PSD matrix.
%     For square matrix X, TRACE_INV(X) is TRACE(INV(X)) if X is Hermitian
%     or symmetric and positive definite; and +Inf otherwise.
%
%     An error results if X is not a square matrix.
%
%     Disciplined convex programming information:
%         TRACE_INV is convex and nonmonotonic (at least with respect to
%         elementwise comparison), so its argument must be affine.

error( nargchk( 1, 1, nargin ) );
if ndims( Y ) > 2 || size( Y, 1 ) ~= size( Y, 2 ),
    error( 'Input must be a square matrix.' );
end
err = Y - Y';
Y   = 0.5 * ( Y + Y' );
if norm( err, 'fro' )  > 8 * eps * norm( Y, 'fro' ),
    z = Inf;
else
    z = eig( full( Y ) );
    if any( z <= 0 ),
        z = Inf;
    else
        z = sum(1.0./z);
    end
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function [] = Example_ABOD()
    fprintf('In this demo, we will show you how to use ABOD and fastABOD.\n');
    fprintf('\n');
    fprintf('Behold, red "+" expresses anomaly points and blue "+" expresses normal points.\n');
    fprintf('\n');
    fprintf('ABOD with baseline dataset demo will start.\n');
    fprintf('\n');
    fprintf('Baseline.mat contains 6 points.\n');
    fprintf('\n');
    fprintf(' 1  0\n');
    fprintf(' 0  1\n');
    fprintf(' 0  0\n');
    fprintf('-1  0\n');
    fprintf(' 0 -1\n');
    fprintf(' 5  5\n');
    fprintf('\n');
    fprintf('First 5 points are very close to each other and aroind the origin,\n');
    fprintf('and 6th point is far away from origin.\n');
    fprintf('\n');
    fprintf('So obviously 6th will be treat as anomaly.\n');
    fprintf('\n');

    baseline = [1 0;0 1;0 0;-1 0;0 -1;5 5];
    [suspicious_index abof]=ABOD(baseline);
    anomaly=suspicious_index(1);
    normal=suspicious_index(2:end);
    plot(baseline(anomaly,1),baseline(anomaly,2),'r +',baseline(normal,1),baseline(normal,2),'b +');
    axis([-2, 6, -2, 6]);

    fprintf('\n');
    fprintf('\n');
    fprintf('\n');
    fprintf('Press Enter before continuing\n');
    pause;

    fprintf('Loading Iris dataset from UCI\n');
    fprintf('\n');

    load('iris.mat', '-mat');
    fprintf('Processing...\n');
    fprintf('\n');
    fprintf('\n');
    fprintf('\n');

    [suspicious_index abof]=ABOD(iris(:,[1;2;5]));
    anomaly=suspicious_index(1:15,1);
    normal=suspicious_index(16:end,1);
    plot3(iris(anomaly(:,1),1),iris(anomaly(:,1),2),iris(anomaly(:,1),5),'r +',iris(normal(:,1),1),iris(normal(:,1),2),iris(normal(:,1),5),'b +');

    fprintf('ABOD with iris dataset demo(1) will start.\n');
    fprintf('\n');
    fprintf('In iris demo(1), only sepal length, sepal width, class will be considered.\n');
    fprintf('\n');
    fprintf('Other attributes will be used in demo(2), so that we can plot points in 3d graph for you.\n');
    fprintf('\n');
    fprintf('Because this dataset contains 3 class, you will see 3 groups in graph.\n');
    fprintf('\n');
    fprintf('And each group has its anomaly.\n');
    fprintf('\n');

    fprintf('\n');
    fprintf('\n');
    fprintf('\n');
    fprintf('Press Enter before ABOD with iris dataset demo(2).\n');
    pause;

    fprintf('Processing...\n');
    fprintf('\n');
    fprintf('\n');
    fprintf('\n');

    [suspicious_index abof]=ABOD(iris(:,[3;4;5]));
    anomaly=suspicious_index(1:15,1);
    mormal=suspicious_index(16:end,1);
    plot3(iris(anomaly(:,1),3),iris(anomaly(:,1),4),iris(anomaly(:,1),5),'r +',iris(normal(:,1),1),iris(normal(:,1),3),iris(normal(:,1),3),'b +');

    fprintf('In iris demo(2), only petal length, petal width, class will be considered.\n');
    fprintf('\n');

    fprintf('Now you already know how to use ABOD.\n');
    fprintf('\n');
    fprintf('fastABOD is a approximation version of ABOD, so basicly they are same thing.\n');
    fprintf('\n');
    fprintf('You can find all detail information in ABOD.m and fastABOD.m\n');
    fprintf('\n');
    fprintf('Enjoy! :)\n');
    fprintf('\n');

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox


text:
function [tri, bc] = gridDataFast2D(x, y, xi, yi)
% Copyright (c) 2012, Chao Huang
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% enforce x, y, xi, yi to be column vectors
x = x(:);
y = y(:);
xi = xi(:);
yi = yi(:);

% check if the MATLAB version is new enough to have DelaunayTri
if exist('DelaunayTri') %#ok<EXIST>
    % triangulize the data
    tri = DelaunayTri(x, y);
else
    % call the old version of gridDataFast using tsearch
    [tri, bc] = gridDataFast2D_tsearch(x, y, xi, yi);
    return
end

% catch trinagulation error
if isempty(tri)
    error('Data cannot be triangulated.');
end

% find the nearest triangle and the corresponding Barycentric coordinates
[t, bc] = pointLocation(tri,[xi yi]);

% check points are valid
if any(isnan(t))
    error('Cartesian points must lie within the k-space grid defined by kgrid');
end

tri = tri(t,:);

repo: Aero-matlab


text:
%h = guidata(6);
h = guidata(gcf);
fg = h.fiberGroups(h.curFiberGroup);
% interp method can be 'nearest'|'trilin'|'spline'. If you are going to
% just average a bunch of these values, it's probably best to use nearest.
interp = 'nearest';
md = dtiGetValFromFibers(h.dt6, fg, inv(dtiGet(h,'dt6xform')),'md',interp);
fa = dtiGetValFromFibers(h.dt6, fg, inv(dtiGet(h,'dt6xform')),'fa',interp);
ln = dtiGetValFromFibers(h.dt6, fg, inv(dtiGet(h,'dt6xform')),'linearity',interp);

%alignSlice = [nan -55 nan]; % the alignment position (in mm from AC)
alignSlice = [nan nan 28]; % the alignment position (in mm from AC)
ax = find(~isnan(alignSlice));
sl = alignSlice(ax);
if(length(ax)~=1) error('specify just one alignment axis!'); end
if(ax==1) axName='LR'; elseif(ax==2) axName='AP'; else axName='SI'; end
nPts = 20;
nFibers = length(fg.fibers);
fiberFa = zeros(nFibers, nPts*2+1)*nan;
fiberMd = fiberFa; fiberLn = fiberFa;
fiberCoords = zeros(nFibers, nPts*2+1,3)*nan;
for(ii=1:length(fg.fibers))
    d = abs(fg.fibers{ii}(ax,:)-alignSlice(ax));
    nearInd = find(d==min(d));
    if(nearInd-nPts<=1 | nearInd+nPts>=length(fg.fibers{ii}))
        continue;
    end
    fDir = fg.fibers{ii}(ax,[nearInd-1,nearInd+1]);
    inds = [nearInd-nPts:nearInd+nPts];
    % We skip the fiber end-points, since the FA there can be misleading
    if(fDir(2)<fDir(1)) inds = fliplr(inds); end
    goodInds = inds>1 & inds<length(d);
    fiberFa(ii,goodInds) = fa{ii}(inds(goodInds));
    fiberMd(ii,goodInds) = md{ii}(inds(goodInds));
    fiberLn(ii,goodInds) = ln{ii}(inds(goodInds));
    fiberCoord(ii,goodInds,:) = fg.fibers{ii}(:,inds(goodInds))';
end

xlab = [axName ' step (mm)'];
x = [sl-nPts:sl+nPts];
%x = mean(fiberCoord(:,:,ax),1);

mnFa = nanmean(fiberFa);
sdFa = nanstd(fiberFa);
figure;errorbar(x,mnFa,sdFa,sdFa);
hold on;
plot(x,min(fiberFa),'c');
hold off;
y = get(gca,'YLim');
line([sl,sl],y,'color','r');
xlabel(xlab); ylabel('Mean FA');

mnMd = nanmean(fiberMd);
sdMd = nanstd(fiberMd);
figure;errorbar(x,mnMd,sdMd,sdMd);
y = get(gca,'YLim');
line([sl,sl],y,'color','r');
xlabel(xlab); ylabel('Mean Diffusivity (mm^2/s*10^-^6)');

mnLn = nanmean(fiberLn);
sdLn = nanstd(fiberLn);
figure;errorbar(x,mnLn,sdLn,sdLn);
y = get(gca,'YLim');
line([sl,sl],y,'color','r');
xlabel(xlab); ylabel('Linearity Index');

error('stop here');

tdir = '/biac2/wandell2/data/reading_longitude/templates/child/SIRL53warp3';
h = guidata(1);
fg = h.fiberGroups(h.curFiberGroup);
alignSlice = [nan nan 28]; % the alignment position (in mm from AC)
ax = find(~isnan(alignSlice));
sl = alignSlice(ax);
if(length(ax)~=1) error('specify just one alignment axis!'); end
if(ax==1) axName='LR'; elseif(ax==2) axName='AP'; else axName='SI'; end
nPts = 20;
nFibers = length(fg.fibers);
% interp method can be 'nearest'|'trilin'|'spline'. If you are going to
% just average a bunch of these values, it's probably best to use nearest.
interp = 'nearest';

snFiles = findSubjects(tdir, '*_sn*',{});
N = length(snFiles);
gFiberFa = zeros(nFibers, nPts*2+1, N)*nan;
gFiberMd = fiberFa; fiberLn = fiberFa;
gFiberCoords = zeros(nFibers, nPts*2+1,3)*nan;
for(subNum=1:N)
    disp(['Loading ' snFiles{subNum} '...']);
    dt = load(snFiles{subNum});
    dt.dt6(isnan(dt.dt6)) = 0;
    md = dtiGetValFromFibers(dt.dt6, fg, inv(dt.xformToAcPc), 'md', interp);
    fa = dtiGetValFromFibers(dt.dt6, fg, inv(dt.xformToAcPc), 'fa', interp);
    ln = dtiGetValFromFibers(dt.dt6, fg, inv(dt.xformToAcPc), 'linearity', interp);

    for(ii=1:length(fg.fibers))
        d = abs(fg.fibers{ii}(ax,:)-alignSlice(ax));
        nearInd = find(d==min(d));
        fDir = fg.fibers{ii}(ax,[nearInd-1,nearInd+1]);
        inds = [nearInd-nPts:nearInd+nPts];
        % We skip the fiber end-points, since the FA there can be misleading
        if(fDir(2)<fDir(1)) inds = fliplr(inds); end
        goodInds = inds>1 & inds<length(d);
        gFiberFa(ii,goodInds,subNum) = fa{ii}(inds(goodInds));
        gFiberMd(ii,goodInds,subNum) = md{ii}(inds(goodInds));
        gFiberLn(ii,goodInds,subNum) = ln{ii}(inds(goodInds));
        gFiberCoord(ii,goodInds,:) = fg.fibers{ii}(:,inds(goodInds))';
    end
end

ssMnFa = nanmean(fiberFa, 1);
%sdFa = nanstd(fiberFa);
gMnFa = squeeze(nanmean(gFiberFa, 1));
mnGMnFa = nanmean(gMnFa, 2);
sdGMnFa = nanstd(gMnFa, 0, 2);
semGMnFa = sdGMnFa/sqrt(N-1);
figure;errorbar(x, mnGMnFa, semGMnFa, semGMnFa);
hold on;
plot(x,ssMnFa,'r');
hold off;
y = get(gca,'YLim');
line([sl,sl],y,'color','r');
xlabel(xlab); ylabel('Mean FA');

ssMnMd = nanmean(fiberMd, 1);
%sdFa = nanstd(fiberFa);
gMnMd = squeeze(nanmean(gFiberMd, 1));
mnGMnMd = nanmean(gMnMd, 2);
sdGMnMd = nanstd(gMnMd, 0, 2);
semGMnMd = sdGMnMd/sqrt(N-1);
figure;errorbar(x, mnGMnMd, semGMnMd, semGMnMd);
hold on;
plot(x,ssMnMd,'r');
hold off;
y = get(gca,'YLim');
line([sl,sl],y,'color','r');
xlabel(xlab); ylabel('Mean MD (mm^2/s*10^-^6)');

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [gradnew] = ft_convert_grad(grad, amplitude, distance, scaling)

% FT_CONVERT_GRAD is deprecated, please use ft_datatype_sens

% DEPRECATED by roboos on 1 July 2013
% see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=963 for more details
% this function can be completely removed medio 2014
warning('this function is deprecated, please use ft_datatype_sens');

% at the time of writing this code, this is only supported in the "upcoming" version, not the "latest/2011v2"
gradnew = ft_datatype_sens(grad, 'version', 'upcoming', 'amplitude', amplitude, 'distance', distance, 'scaling', scaling);


text:
function [ gc_sum, gc_avg, gc_std  ] = gcFeature( gc, indicator, coor2D_indi )
%GCFEATURE Summary of this function goes here
%   Detailed explanation goes here
if size(indicator,1) ~= size(coor2D_indi,1)
    fprintf('Warning: indicator and coor mismatch!\n');
end
% vector_num = size(indicator,1);
data_num = size(indicator, 2);

gc_sum = zeros( data_num, 11);
gc_avg = zeros( data_num, 11);
gc_std = zeros( data_num, 11);

[sH, sW, ~] = size(gc);
% feature = zeros(data_num, 100);

gc = cat(3, gc, sum(gc(:,:,1:4),3), max(gc(:,:,1:4),[],3), sum(gc(:,:,1:6),3), max(gc(:,:,1:6),[],3));
for did = 1:data_num
    vector_valid = indicator(:,did);
    if ~any(vector_valid)
        continue;
    end
    cood_ind = coor2D_indi(vector_valid);
    loc_orit = zeros(length(cood_ind),11);
    for k = 1:11
        loc_orit(:,k) = gc(cood_ind+(k-1)*sW*sH);
    end

    gc_sum(did,:) = sum(loc_orit, 1);
    gc_avg(did,:) = gc_sum(did,:)./length(cood_ind);
    gc_std(did,:) = std(loc_orit, 1);
end

end


text:
function [W, b] = initializeOneLayerParams(layerSize, inputSize)
%% Initialize Parameters Randomly Based on Layer Sizes
r  = sqrt(6) / sqrt(layerSize+inputSize+1);
W = rand(layerSize, inputSize) * 2 * r - r;
b = zeros(layerSize, 1);
end

author: zellyn


text:
function idxRRtoBeRemoved = FindSpikesInRR(RR, th)
%
% idxRRtoBeRemoved = FindSpikesInRR(RR, th)
%
% OVERVIEW : Code used to clean RR intervals that change more than a given
% threshold (eg., th = 0.2 = 20%) with respect to the median value of the
% previous 5 and next 5 RR intervals (using a forward-backward approach).
%
% INPUTS:
%       RR : a single row of rr interval data in seconds
%       th : threshold percent limit of change from one interval to the next
% OUTPUTS:
%       idxRRtoBeRemoved : a single vector of indexes related to RR
%                          intervals corresponding to a change > th
%
%   DEPENDENCIES & LIBRARIES:
%       PhysioNet Cardiovascular Signal Toolbox
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%
%   REFERENCE:
%       Vest et al. "An Open Source Benchmarked HRV Toolbox for Cardiovascular
%       Waveform and Interval Analysis" Physiological Measurement (In Press), 2018.
%
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%
%   Written by Giulia Da Poian (giulia.dap@gmail.com), 09-13-2017
%
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%

if size(RR,1)>size(RR,2)
    RR = RR';
end

% Forward search
FiveRR_MedianVal = medfilt1(RR,5); % compute as median RR(-i-2: i+2)
% shift of three position to aligne with to corresponding RR
FiveRR_MedianVal = [RR(1:5) FiveRR_MedianVal(3:end-3)];
rr_above_th = find((abs(RR-FiveRR_MedianVal)./FiveRR_MedianVal)>=th);

RR_forward = RR;
RR_forward(rr_above_th) = NaN;

% Backward search
RRfilpped = fliplr(RR);

FiveRR_MedianVal = medfilt1(RRfilpped,5); % compute as median RR(-i-2: i+2)
% shift of three position to aligne with to corresponding RR
FiveRR_MedianVal = [RRfilpped(1:5) FiveRR_MedianVal(3:end-3)];
rr_above_th = find(abs(RRfilpped-FiveRR_MedianVal)./FiveRR_MedianVal>=th);
rr_above_th = sort(length(RR)-rr_above_th+1);

RR_backward = RRfilpped;
RR_backward(rr_above_th) = NaN;

% Combine

idxRRtoBeRemoved = (isnan(RR_forward) & isnan(RR_backward));

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [psth] = ft_spike_psth(cfg, spike)

% FT_SPIKE_PSTH computes the peristimulus histogram of spiketrains.
%
% Use as
%   [psth] = ft_spike_psth(cfg, spike)
%
% The input SPIKE should be organised as either the spike datatype,
% obtained from FT_SPIKE_MAKETRIALS, or the raw datatype, containing binary
% spike trains, obtained from FT_APPENDSPIKE or FT_CHECKDATA. In this case
% the raw datatype is converted to the spike datatype.
%
% Configurations:
%   cfg.binsize          =  [binsize] in sec or string.
%                          If 'scott', we estimate the optimal bin width
%                          using Scott's formula (1979). If 'sqrt', we take
%                          the number of bins as the square root of the
%                          number of observations. The optimal bin width is
%                          derived over all neurons; thus, this procedure
%                          works best if the input contains only one neuron
%                          at a time.
%   cfg.outputunit       = 'rate' (default) or 'spikecount' or
%                          'proportion'. If 'rate', we
%                          convert the output per trial to firing rates
%                          (spikes/sec). If 'spikecount', we count the
%                          number spikes per trial. If 'proportion', we
%                          normalize the area under the PSTH to 1.
%   cfg.spikechannel     = See FT_CHANNELSELECTION for details. cfg.trials
%                          is vector of indices (e.g., 1:2:10)
%                          logical selection of trials (e.g., [1010101010])
%                          'all' (default), selects all trials
%   cfg.vartriallen      = 'yes' (default)
%                          Accept variable trial lengths and use all
%                          available trials and the samples in every trial.
%                          Missing values will be ignored in the
%                          computation of the average and the variance and
%                          stored as NaNs in the output psth.trial. 'no'
%                          Only select those trials that fully cover the
%                          window as specified by cfg.latency and discard
%                          those trials that do not.
%   cfg.latency          = [begin end] in seconds
%                          'maxperiod' (default), i.e., maximum period
%                          available 'minperiod', i.e., the minimal period
%                          all trials share, 'prestim' (all t<=0) 'poststim'
%                          (all t>=0).
%   cfg.keeptrials       = 'yes' or 'no' (default)
%   cfg.trials           =  numeric or logical selection of trials (default = 'all')
%
% Outputs:
%   Psth is a timelock datatype (see FT_DATATYPE_TIMELOCK)
%     Psth.time        = center histogram bin points
%	    Psth.fsample     = 1/binsize;
%     Psth.avg         = contains average PSTH per unit
%     Psth.trial       = contains PSTH per unit per trial
%     Psth.var         = contains variance of PSTH per unit across trials
%
% For subsequent processing you can use
%   FT_SPIKE_PLOT_PSTH    : plot only the PSTH, for a single neuron
%   FT_TIMELOCKSTATISTICS : compute statistics on the PSTH
%   FT_SPIKE_PLOT_RASTER  : plot PSTH with raster for one or more neurons
%   FT_SPIKE_JPSTH        : compute the JPSTH

%  Copyright (C) 2010-2013, Martin Vinck
%
% 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$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble provenance spike

% control input spike structure and convert to spike if raw structure
spike = ft_checkdata(spike,'datatype', 'spike', 'feedback', 'yes');

% get the default options
cfg.outputunit   = ft_getopt(cfg, 'outputunit','rate');
cfg.binsize      = ft_getopt(cfg, 'binsize','scott');
cfg.spikechannel = ft_getopt(cfg, 'spikechannel', 'all');
cfg.trials       = ft_getopt(cfg, 'trials', 'all');
cfg.latency      = ft_getopt(cfg, 'latency','maxperiod');
cfg.vartriallen  = ft_getopt(cfg, 'vartriallen', 'yes');
cfg.keeptrials   = ft_getopt(cfg, 'keeptrials', 'no');

% ensure that the options are valid
cfg = ft_checkopt(cfg,'outputunit',   'char',  {'rate', 'spikecount'});
cfg = ft_checkopt(cfg,'binsize',     {'char',  'doublescalar'});
cfg = ft_checkopt(cfg,'spikechannel',{'cell',  'char', 'double'});
cfg = ft_checkopt(cfg,'latency',     {'char',  'ascendingdoublebivector'});
cfg = ft_checkopt(cfg,'trials',      {'char',  'doublevector', 'logical'});
cfg = ft_checkopt(cfg,'vartriallen' , 'char', {'yes', 'no'});
cfg = ft_checkopt(cfg,'keeptrials'  , 'char', {'yes', 'no'});

cfg = ft_checkconfig(cfg, 'allowed', {'outputunit', 'binsize', 'spikechannel', 'trials', 'latency', 'vartriallen', 'keeptrials'});

% get the number of trials or convert to indices
cfg        = trialselection(cfg,spike);

% select the unit - this should be done with channelselection function
cfg.spikechannel = ft_channelselection(cfg.spikechannel, spike.label);
spikesel         = match_str(spike.label, cfg.spikechannel);
nUnits           = length(spikesel);
if nUnits==0, error('no spikechannel selected by means of cfg.spikechannel'); end

% determine the duration of each trial
begTrialLatency = spike.trialtime(cfg.trials,1); % remember: already selected on trial here
endTrialLatency = spike.trialtime(cfg.trials,2);
trialDur 	    	= endTrialLatency - begTrialLatency;

% select the latencies, use the same modular function in all the scripts
cfg = latencyselection(cfg,begTrialLatency,endTrialLatency);

% compute the optimal bin width if desired
if ischar(cfg.binsize)
  h = zeros(1,nUnits);
  for iUnit = 1:nUnits
    unitIndx       = spikesel(iUnit); % select the unit
    spikesInWin    = spike.time{unitIndx}>=cfg.latency(1) & spike.time{unitIndx}<=cfg.latency(2); % get spikes in trial

    % automatically determine an 'optimal' binwidth
    N              = sum(spikesInWin);
    if strcmp(cfg.binsize,'scott')
      sd           = nanstd(spike.time{unitIndx}(spikesInWin));
      h(iUnit)     = 3.49*sd./(N^(1/3));
    elseif strcmp(cfg.binsize,'sqrt')
      k            = ceil(sqrt(N));
      h(iUnit)     = (cfg.latency(2)-cfg.latency(1))/k;
    else
      error('unsupported option for cfg.binsize');
    end
  end
  cfg.binsize = nanmean(h);
end

% do some error checking on the binsize
if cfg.binsize<=0 || cfg.binsize>(cfg.latency(2)-cfg.latency(1))
  error('cfg.binsize should be greater than zero and not exceed the trialduration');
end

% end of trial should be late enough, beginning should be early enough
fullDur       = trialDur>=cfg.binsize; % trials which have full duration
overlaps      = endTrialLatency>=(cfg.latency(1)+cfg.binsize) & begTrialLatency<=(cfg.latency(2)-cfg.binsize);
if strcmp(cfg.vartriallen,'no') % only select trials that fully cover our latency window
  startsLater    = single(begTrialLatency) > (single(cfg.latency(1)) + 0.0001); % last factor just for rounding errors
  endsEarlier    = single(endTrialLatency) < (single(cfg.latency(2)) - 0.0001);
  hasWindow      = ~(startsLater | endsEarlier); % it should not start later or end earlier
else
  hasWindow     = ones(1,length(cfg.trials));
end
trialSel        = fullDur(:) & overlaps(:) & hasWindow(:);
cfg.trials      = cfg.trials(trialSel); % note that endTrialLatency was of length cfg.trials
if isempty(cfg.trials), warning('No trials were selected after latency selection'); end
nTrials         = length(cfg.trials);
begTrialLatency = begTrialLatency(trialSel); % note that begTrialLatency was of length cfg.trials here
endTrialLatency = endTrialLatency(trialSel);

% create the bins, we start arbitrarily at cfg.latency(1)
bins  = cfg.latency(1) : cfg.binsize : cfg.latency(end);
nBins = length(bins);

% preallocate before computing the psth, otherwise, compute stuff 'on the run'
if strcmp(cfg.keeptrials,'yes'), singleTrials = NaN(nTrials,nUnits,nBins-1); end
[s,ss]   = deal(zeros(nUnits, nBins-1)); %sum, and sum of squared in the loop, to get mean & var

% preallocate and compute degrees of freedom
allStartEarlier =  all(begTrialLatency<=cfg.latency(1));
allEndLater     =  all(endTrialLatency>=cfg.latency(2));
if  ~ (allStartEarlier && allEndLater),
  dof = zeros(1,nBins-1);
else
  dof = ones(1,nBins-1)*nTrials; % if all exceed latency, then this is dof
end

for iTrial = 1:nTrials
  origTrial = cfg.trials(iTrial);
  if  ~ (allStartEarlier && allEndLater) % select bins and count dof + 1
    binSel = begTrialLatency(iTrial)<=bins(1:end-1) & endTrialLatency(iTrial)>=bins(2:end);
    dof(binSel)      = dof(binSel) + 1;
  else
    binSel           = 1:(nBins-1); % always deselect the last bin
  end

  for iUnit = 1:nUnits
    unitIndx      = spikesel(iUnit); % select the unit
    spikesInTrial = (spike.trial{unitIndx}==origTrial); % get spikes in trial
    spikeTimes    = spike.time{unitIndx}(spikesInTrial);

    % compute the psth
    trialPsth   = histc(spikeTimes(:), bins); % we deselect the last bin per default
    trialPsth   = trialPsth(:)'; % force into row vec
    if isempty(trialPsth), trialPsth = zeros(1,length(bins)); end

    % convert to firing rates if requested, with spikecount do nothing
    if strcmp(cfg.outputunit,'rate'),
      trialPsth = trialPsth/cfg.binsize;
    elseif strcmp(cfg.outputunit,'proportion'),
      trialPsth = trialPsth./nansum(trialPsth);
    end

    % compute the sum and the sum of squares for the var and the mean on the fly
    s(iUnit,binSel)  = s(iUnit,binSel)   + trialPsth(binSel); % last bin is single value
    ss(iUnit,binSel) = ss(iUnit,binSel)  + trialPsth(binSel).^2;

    if strcmp(cfg.keeptrials,'yes'),
      singleTrials(iTrial,iUnit,binSel) = trialPsth(binSel);
    end
  end
end

% get the results
dof            = dof(ones(nUnits,1),:);
psth.avg       = s ./ dof;
psth.var       = (ss - s.^2./dof)./(dof-1); % since sumPsth.^2 ./ dof = dof .* (sumPsth/dof).^2
psth.dof       = dof; % combined with psth.var we can get SEM
psth.fsample   = 1/(cfg.binsize);   % might be more compatible with feeding psth in other funcs
psth.time      = bins(1:end-1) + 0.5*cfg.binsize;
psth.label     = spike.label(spikesel);
if (strcmp(cfg.keeptrials,'yes'))
  psth.trial  = singleTrials;
  psth.dimord = 'rpt_chan_time';
else
  psth.dimord = 'chan_time';
end
if isfield(spike,'sampleinfo'), psth.sampleinfo = spike.sampleinfo(cfg.trials,:); end
if isfield(spike,'trialinfo'),  psth.trialinfo  = spike.trialinfo(cfg.trials,:);  end

% do the general cleanup and bookkeeping at the end of the function

ft_postamble previous   spike
ft_postamble provenance psth
ft_postamble history    psth

%%%%%%%%% SUB FUNCTIONS %%%%%%%%%
function [cfg] = latencyselection(cfg,begTrialLatency,endTrialLatency)

if strcmp(cfg.latency,'minperiod')
  cfg.latency = [max(begTrialLatency) min(endTrialLatency)];
elseif strcmp(cfg.latency,'maxperiod')
  cfg.latency = [min(begTrialLatency) max(endTrialLatency)];
elseif strcmp(cfg.latency,'prestim')
  cfg.latency = [min(begTrialLatency) 0];
elseif strcmp(cfg.latency,'poststim')
  cfg.latency = [0 max(endTrialLatency)];
end
% check whether the time window fits with the data
if (cfg.latency(1) < min(begTrialLatency)), cfg.latency(1) = min(begTrialLatency);
  warning('Correcting begin latency because it is before all trial beginnings');
end
if (cfg.latency(2) > max(endTrialLatency)), cfg.latency(2) = max(endTrialLatency);
  warning('Correcting end latency because it is after all trial ends');
end

function [cfg] = trialselection(cfg,spike)

nTrials = size(spike.trialtime,1);
if  strcmp(cfg.trials,'all')
  cfg.trials = 1:nTrials;
elseif islogical(cfg.trials) || all(cfg.trials==0 | cfg.trials==1)
  cfg.trials = find(cfg.trials);
end
cfg.trials = sort(cfg.trials(:));
if max(cfg.trials)>nTrials,
  error('maximum trial number in cfg.trials should not exceed number of rows of spike.trialtime');
end
if isempty(cfg.trials), error('No trials were selected by you, rien ne va plus'); end


text:
function [t,x] = prz(bits, bitrate)
% PRZ Encode bit string using polar RZ code.
%   [T, X] = PRZ(BITS, BITRATE) encodes BITS array using polar RZ
%   code with given BITRATE. Outputs are time T and encoded signal
%   values X.

% Copyright (c) 2013 Yuriy Skalko <yuriy.skalko@gmail.com>

T = length(bits)/bitrate; % full time of bit sequence
n = 200;
N = n*length(bits);
dt = T/N;
t = 0:dt:T;
x = zeros(1,length(t)); % output signal

for i = 0:length(bits)-1
  if bits(i+1) == 1
    x(i*n+1:(i+0.5)*n) = 1;
    x((i+0.5)*n+1:(i+1)*n) = 0;
  else
    x(i*n+1:(i+0.5)*n) = -1;
    x((i+0.5)*n+1:(i+1)*n) = 0;
  end
end

author: Sable

repo: mcbench-benchmarks


text:
if(exist('E:\datasets\FERA_2015\semaine/SEMAINE-Sessions/', 'file'))
    SEMAINE_dir = 'E:\datasets\FERA_2015\semaine/SEMAINE-Sessions/';
elseif(exist('I:\datasets\FERA_2015\Semaine\SEMAINE-Sessions/', 'file'))
    SEMAINE_dir = 'I:\datasets\FERA_2015\Semaine\SEMAINE-Sessions/';
elseif(exist('C:\tadas\face_datasets\fera_2015\semaine/SEMAINE-Sessions/', 'file'))
    SEMAINE_dir = 'C:\tadas\face_datasets\fera_2015\semaine/SEMAINE-Sessions/';
elseif(exist('D:\datasets\face_datasets\fera_2015\semaine\SEMAINE-Sessions/', 'file'))
    SEMAINE_dir = 'D:\datasets\face_datasets\fera_2015\semaine\SEMAINE-Sessions/';
elseif(exist('D:\Datasets\FERA_2015\semaine\SEMAINE-Sessions/', 'file'))
    SEMAINE_dir = 'D:\Datasets\FERA_2015\semaine\SEMAINE-Sessions/';
elseif(exist('D:/fera_2015/semaine/SEMAINE-Sessions/', 'file'))
    SEMAINE_dir = 'D:/fera_2015/semaine/SEMAINE-Sessions/';
else
    fprintf('DISFA location not found (or not defined)\n');
end

if(exist('SEMAINE_dir', 'var'))
    hog_data_dir = 'E:\datasets\face_datasets_processed\semaine/';
end

train_recs = {'rec1', 'rec12', 'rec14', 'rec19', 'rec23', 'rec25', 'rec37', 'rec39', 'rec43', 'rec45', 'rec48', 'rec50', 'rec52', 'rec54', 'rec56', 'rec60'};
devel_recs = {'rec9', 'rec13', 'rec15', 'rec20', 'rec24', 'rec26', 'rec38', 'rec42', 'rec44', 'rec46', 'rec49', 'rec51', 'rec53', 'rec55', 'rec58'};

author: TadasBaltrusaitis

repo: OpenFace


text:
function Gimg = gaborConvF(img,G,gWinLen)
%GIMG = GABORCONV(IMG,G,GWINLEN) Computes Gabor transform using FFT
%   G is the generated FFTed Gabor kernels using GENGABORKERNELF, GIMG will
%   be a cell the same size w/ G. GWINLEN is the radius of the origin kernel
%	Yan Ke @ THUEE, xjed09@gmail.com

[scale_num angle_num] = size(G);
Gimg = cell(scale_num,angle_num);
[fftM fftN] = size(G{1});
[imgM imgN] = size(img);

f = padarray(img,[gWinLen gWinLen],'rep');
fimg = fft2(f,fftM,fftN);
for r = 1:scale_num
	for s = 1:angle_num
		filtered = ifft2(fimg.*G{r,s}); % conv in freq domain
		cropped = filtered(gWinLen*2+(1:imgM),gWinLen*2+(1:imgN));
		Gimg{r,s} = abs(cropped); % use magnitude
	end
end

end


text:
function [hdr] = read_wdq_header(filename)

% READ_WDQ_HEADER reads header information from wdq files
%
% Use as
%  [hdr] = read_wdq_header(filename)
%

% Copyright (C) 2011, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% 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: read_wdq_header.m $

% information about how to interpret the file are taken from the document
% 'CODAS data storage format'

fid = fopen(filename, 'r');
%tmp = fread(fid, 110); % the first 110 bytes are always interpretable equally
hdr = [];
hdr.nchan         = fread(fid, 1, 'uint16');
hdr.nchan         = hdr.nchan - bitand(hdr.nchan, 32);
i2                = fread(fid, 1, 'uint16');
hdr.offset        = fread(fid, 1, 'uint8');
hdr.nbyteschanhdr = fread(fid, 1, 'uint8');
hdr.nbyteshdr     = fread(fid, 1, 'int16');
hdr.nbytesdat     = fread(fid, 1, 'uint32'); %FIXME only for unpacked files
hdr.nbytestrailer1 = fread(fid, 1, 'uint32');
hdr.nbytestrailer2 = fread(fid, 1, 'uint16');
i9  = fread(fid, 1, 'uint16'); % not used
i10 = fread(fid, 1, 'uint16'); % not used
i11 = fread(fid, 1, 'int16');
i12 = fread(fid, 4, 'int8');
i13 = fread(fid, 1, 'double');

% FIXME there seems to be a discrepancy here, according to the
% documentation it should be the next (commented) line. according to the
% file I have at hand it should be the uncommented line.
%hdr.fsample        = double(hdr.nchan)./i13;
hdr.fsample = 1./i13;
i14 = fread(fid, 1, 'uint32');
i15 = fread(fid, 1, 'uint32');
i16 = fread(fid, 1, 'uint32');
i17 = fread(fid, 1, 'uint32');
i18 = fread(fid, 1, 'uint32');
i19 = fread(fid, 2, 'uint16');
i20 = fread(fid, 1, 'int16');
i21 = fread(fid, 1, 'int16');
i22 = fread(fid, 1, 'int8');
i23 = fread(fid, 1, 'int8');
i24 = fread(fid, 1, 'int8');
i25 = fread(fid, 1, 'int8');
i26 = fread(fid, 32, 'int8');
i27 = fread(fid, 1, 'uint16');
i28 = fread(fid, 1, 'uint8');
empty = fread(fid, 1, 'uint8');
i29 = fread(fid, 1, 'int8');
i30 = fread(fid, 1, 'int8');
i31 = fread(fid, 1, 'int8');
empty = fread(fid, 1, 'uint8');
i32 = fread(fid, 1, 'int8');
i33 = fread(fid, 1, 'int8');

% interpret relevant stuff from i27
% other stuff may or may not be relevant
hdr.ispacked  = bitand(i27, 2^15)>0;
hdr.lowchan   = bitand(i27, 2^9) >0;

% get channel info
for k = 1:(hdr.nbyteshdr-hdr.offset-2)/hdr.nbyteschanhdr
  fseek(fid, hdr.offset + (k-1)*hdr.nbyteschanhdr, 'bof');
  chanhdr(k).scalingslope     = fread(fid, 1, 'float');
  chanhdr(k).scalingintercept = fread(fid, 1, 'float');
  chanhdr(k).scale            = fread(fid, 1, 'double');
  chanhdr(k).intercept        = fread(fid, 1, 'double');
  chanhdr(k).unit             = fread(fid, 1, 'uint8=>char');
  fread(fid, 1, 'uint8');
  fread(fid, 1, 'uint8'); % may be relevant for packed files
  chanhdr(k).channr = fread(fid, 1, 'uint8');
  chanhdr(k).label  = num2str(k);
end
hdr.chanhdr = chanhdr;

fseek(fid, hdr.nbyteshdr + hdr.nbytesdat, 'bof');
x = fread(fid, hdr.nbytestrailer1, 'uint8');
hdr.x = x; % what shall we do with this?

fseek(fid, hdr.nbyteshdr + hdr.nbytesdat + hdr.nbytestrailer1, 'bof');
annot = fread(fid, hdr.nbytestrailer2, 'uint8');
sel   = find(annot==0);
indx  = find(annot>0, 1, 'first');
sel   = [0;sel];
for k = 1:numel(sel)-1
  hdr.chanhdr(k).annot = char(annot(sel(k)+1:sel(k+1)-1))';
end

fclose(fid);

author: PatternRecognition

repo: OpenBMI


text:
function classplot2random(X,y,flag,colorstr,MarkerSize)
%+++ flag: 0 No marker;
%          1 using clorstr
%          2 using class number;

if nargin<5;MarkerSize=5;end
if nargin<4;colorstr={'g.';'b.';'k*';'gp';'cs';'g>';'c.';};end
if nargin<3;flag=0;end

nclass=unique(y);
n=length(y);
nperm=randperm(n);
hold on;
for i=1:n
  permi=nperm(i);
  if y(permi)==nclass(1);colortemp=colorstr{1};else; colortemp=colorstr{2};end
  plot(X(permi,1),X(permi,2),colortemp,'MarkerSize',MarkerSize);
end

author: viggin

repo: domain-adaptation-toolbox


text:
function [satclk0,varc0,stat]=pephclk(time,sat,nav,satclk0,varc0)

global glc
stat=1;

if nav.nc<2||timediff(time,nav.peph(1).time)<-900||...
        timediff(time,nav.peph(end).time)>900
    return;
end

%binary search
i=0; j=nav.nc-1;
while i<j
    k=fix((i+j)/2);
    if timediff(nav.pclk(k+1).time,time)<0
        i=k+1;
    else
        j=k;
    end
end

if i<=0,idx=0; else,idx=i-1; end
idx=idx+1;

%linear interpolation for clock
t0=timediff(time,nav.pclk(idx).time);
t1=timediff(time,nav.pclk(idx+1).time);
c0=nav.pclk(idx).clk(sat);
c1=nav.pclk(idx+1).clk(sat);

if t0<=0
    satclk=c0;
    if satclk==0,stat=0;return;end
    varc=(nav.pclk(idx).std(sat)*glc.CLIGHT-1E-3*t0)^2;
elseif t1>=0
    satclk=c1;
    if satclk==0,stat=0;return;end
    varc=(nav.pclk(idx+1).std(sat)*glc.CLIGHT+1E-3*t1)^2;
elseif c0~=0 && c1~=0
    satclk=(c1*t0-c0*t1)/(t0-t1);
    if abs(t0)-abs(t1)>0
        varc=(nav.pclk(idx+1).std(sat)*glc.CLIGHT+1E-3*abs(t1))^2;
    else
        varc=(nav.pclk(idx).std(sat)*glc.CLIGHT+1E-3*abs(t0))^2;
    end
else
    stat=0; return;
end

if satclk~=0
    satclk0=satclk;
    varc0=varc;
end

return

repo: GINav


text:
function c = tapas_unitsq_sgm_mu3_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the unit square sigmoid observation model for binary responses
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The unit square sigmoid (ussgm) is the function
%
% f(x) = x^zeta/(x^zeta + (1-x)^zeta) = ustapas_sgm(x; zeta),
%
% where x is in the unit interval, and zeta > 0 is a parameter that determines the shape of the
% sigmoid. Since both its argument and value are always in the unit interval, its graph is
% restricted to the unit square, hence the name unit square sigmoid.
%
% In the application here, the ussgm is the probability of observing a decision y=1 (rather than
% the only alternative y=0) given the current probability mu1hat (or value) of input u=1:
%
% p(y=1|mu1hat) = ustapas_sgm(mu1hat; zeta)
%
% The parameter zeta regulates the steepness of the sigmoid such that it forms the diagonal of
% the unit square for zeta=1 and approaches a step function at 0.5 as zeta approaches infinity.
% Values of 0 < zeta < 1 lead to sigmoids with reverse concavity than usual, but they still
% represent valid observation models.
%
% Zeta can be interpreted as inverse decision noise. To have a shrinkage prior on this, choose a
% high value. It is estimated log-space since it has a natural lower bound at zero.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'tapas_unitsq_sgm_mu3';

% This model has no free parameters

% Gather prior settings in vectors
c.priormus = [
         ];

c.priorsas = [
         ];

% Model filehandle
c.obs_fun = @tapas_unitsq_sgm_mu3;

% Handle to function that transforms observation parameters to their native space
% from the space they are estimated in
c.transp_obs_fun = @tapas_unitsq_sgm_mu3_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
function [out] = ratioScore(x,Prob)

ration = Prob.user.ration;
ratiod = Prob.user.ratiod;
% x should be in alpha coordinates.
% expdata = Prob.user.expdata;
% model = Prob.user.model;
% N = model.N;

out = (ration'*x)/(ratiod'*x);

return;

author: opencobra

repo: cobratoolbox


text:
function matRad_progress(currentIndex, totalNumberOfEvaluations)
% matRad progress bar
%
% call
%   matRad_progress(currentIndex, totalNumberOfEvaluations)
%
% input
%   currentIndex:               current iteration index
%   totalNumberOfEvaluations:   maximum iteration index
%
% output
%   graphical display of progess. make sure there is no other output
%   written during the loop to prevent confusion
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% If it's not the first step, erase the stuff printed before
if (currentIndex == 1 || nargin > 2)
    fprintf('Progress: ');
end

if (currentIndex > 1 && nargin < 3)
  Length = numel(sprintf('%3.2f %%',(currentIndex-1)/totalNumberOfEvaluations*100));
  fprintf(repmat('\b',1,Length));
end

% Print the progress tool
fprintf('%3.2f %%',currentIndex/totalNumberOfEvaluations*100);

% After the last iteration print a newline command
if (currentIndex == totalNumberOfEvaluations)
    fprintf('\n');
end

end

repo: matRad


text:
function enableM(plant)
    plant.isMEnabled = 1;
end

author: deng-haoyang

repo: ParNMPC


text:
function test_bug3369

% WALLTIME 00:20:00
% MEM 2gb
% DEPENDENCY read_neuromag_maxfilterlog

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/bug3369'));

file0 = 'tactile_stim_raw.fif_tsss_mc.log';
file1 = 'tactile_stim_raw-1.fif_tsss_mc.log';
file2 = 'tactile_stim_raw-2.fif_tsss_mc.log';

%%

cfg = [];
cfg.dataset = file0;
cfg.viewmode = 'vertical';
ft_databrowser(cfg)

%%

cfg = [];
cfg.dataset = file0;
data0 = ft_preprocessing(cfg);
cfg.dataset = file1;
data1 = ft_preprocessing(cfg);
cfg.dataset = file2;
data2 = ft_preprocessing(cfg);

cfg = [];
cfg.keepsampleinfo = 'no'; % the sampleinfo of the three does not match
data = ft_appenddata(cfg, data0, data1, data2);

%%

cfg = [];
cfg.viewmode = 'vertical';
cfg.preproc.demean = 'yes';
ft_databrowser(cfg, data);


text:
%% Copyright (C) 2022 Alex Vong
%%
%% 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 <https://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @defun  python_env_is_cygwin_like (pyexec)
%% Check if Python @var{pyexec} is running in a Cygwin-like POSIX environment,
%% such as Cygwin or MSYS2.  The result is memoized to speed up subsequent
%% calls.
%%
%% @seealso{cygpath}
%% @end defun

function r = python_env_is_cygwin_like (pyexec)
  persistent python_env_is_cygwin_like_memo

  if ~isempty (python_env_is_cygwin_like_memo)
    r = python_env_is_cygwin_like_memo;
    return
  end

  if ispc ()
    if system ('where /q cygpath') ~= 0
      r = false;
      python_env_is_cygwin_like_memo = r;
      return
    end

    [status, out] = system ([pyexec ' -c "import os; print(os.name)"']);
    if status ~= 0
      error ('python_env_is_cygwin_like: %s exited with status %d', ...
             pyexec, status);
    end
    r = ~isempty (regexp (out, 'posix', 'match'));
    python_env_is_cygwin_like_memo = r;
  else
    r = false;
    python_env_is_cygwin_like_memo = r;
  end
end

repo: octsympy


text:
function [mm] = cm2mm(cm)
% Convert length from centimeters to millimeters.
mm = cm*10;


text:
function state = evolutioncomplete(options,state,flag)
% Plays a notification when genetic algorithm finishes. Requires a
% Starcraft installation, in its default directory, typically:
% C:\Program Files\Starcraft
%
% Works with Windows, not tested on other platforms.

if ispc && strcmp(flag,'done') && ~strcmp(options.Display,'off')
    sounddir = [getenv('ProgramFiles') '\Starcraft\Sound\Zerg\Advisor'] ;
    if isdir(sounddir) && exist([sounddir '\ZAdUpd02.wav'],'file')
        [y, Fs, nbits] = audioread([sounddir '\ZAdUpd02.wav']) ;
        obj = audioplayer(y, Fs, nbits);
        playblocking(obj)
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function uvHCurrCand = sc_get_prop_cand(uvHCurr, dx, dy)

% Get propagation neighbor

uvHCurrCand = zeros(size(uvHCurr));

Ht = [1, 0, dx; 0, 1, dy; 0 0 1];
numUvPixels = size(uvHCurr, 2);

% for i = 1: numUvPixels
%     H = reshape(uvHCurr(:,i), 3, 3);
%     Hc = H*Ht;
%     uvHCurrCand(:, i) = Hc(:);
% end
%
% uvHCurrCand = uvHCurrCand./repmat(uvHCurrCand(9,:), 9, 1);

end

repo: StructCompletion


text:
% This function sort the utterance ID, and also rearrange the corresponding
% data cell array.
%
function [ID_sort, data_sort] = groupSort(ID, data)

[ID_sort, idx] = sort(ID);
data_sort = data(idx);

author: singaxiong

repo: SignalGraph


text:
function DoEditValue(handles,parent_handle,Attributes,AttrhNumber,dims)

% Create Editbox for adjusting attributes at the frame given by parent_handle
% special mark
% '$value' ----popupmenu
% '^value' ----checkbox
% '@' ----static text

TextW=dims(1);
EditW=dims(2);
SpaceW=dims(3);
TextH=dims(4);
EditH=dims(5);
SpaceH=dims(6);

%--old schema
% switch get(parent_handle,'Type')
%     case 'uitab'
%         Position=get(get(get(parent_handle,'Parent'),'Parent'),'Position');
%     case 'uipanel'
%         Position=get(parent_handle,'Position');
% end
% Hpair=floor(Position(4)/(TextH+SpaceH));
% Wpair=floor(Position(3)/(TextW+SpaceW));
%--end

Hpair=floor(1/(TextH+SpaceH));
Wpair=floor(1/(TextW+SpaceW));
ind=1;
for j=1:Wpair
    for i=1:Hpair
        if ind<=length(Attributes)
            Attrh.([Attributes(ind).Name 'Text'])=uicontrol(parent_handle,'Style', 'text', 'String', Attributes(ind).Name,'FontWeight','bold','Units','normalized',...
                                                                    'Position', [(j-1)*(TextW+SpaceW+EditW) (i-1)*(TextH+SpaceH) TextW TextH],...
                                                                    'TooltipString', DoTellMeInfo(Attributes(ind).Name));
            eval(['handles.Attrh' num2str(AttrhNumber) '.' Attributes(ind).Name 'Text=Attrh.' Attributes(ind).Name 'Text;']);
            if ~isempty(Attributes(ind).Value)
                switch Attributes(ind).Value(1)
                    case '$'
                        eval(['AttributeOpt={' Attributes(ind).Value(3:end) '};']);
                        Attrh.(Attributes(ind).Name)=uicontrol(parent_handle,'Style', 'popupmenu', 'String', AttributeOpt,'BackgroundColor',[1 1 1],...
                                                                        'Value',str2double(Attributes(ind).Value(2)),'Units','normalized',...
                                                                        'Position', [(j-1)*(TextW+SpaceW+EditW)+TextW (i-1)*(TextH+SpaceH) EditW EditH],...
                                                                        'TooltipString',['Default: ' Attributes(ind).Value(2:end)]);
                        eval(['handles.Attrh' num2str(AttrhNumber) '.' Attributes(ind).Name '=Attrh.' Attributes(ind).Name ';']);
                    case '^'
                        Attrh.(Attributes(ind).Name)=uicontrol(parent_handle,'Style', 'checkbox', 'String',[],'BackgroundColor',[1 1 1],...
                                                                        'Value',str2double(Attributes(ind).Value(2)),'Units','normalized',...
                                                                        'Position', [(j-1)*(TextW+SpaceW+EditW)+TextW (i-1)*(TextH+SpaceH) EditW EditH],...
                                                                        'TooltipString',['Default: ' Attributes(ind).Value(2:end)]);
                        eval(['handles.Attrh' num2str(AttrhNumber) '.' Attributes(ind).Name '=Attrh.' Attributes(ind).Name ';']);
                    case '@'
                        eval(['AttributeOpt=' Attributes(ind).Value(2:end) ';']);
                        Attrh.(Attributes(ind).Name)=uicontrol(parent_handle,'Style', 'text', 'String',AttributeOpt,'Units','normalized',...
                                                                        'Position', [(j-1)*(TextW+SpaceW+EditW)+TextW (i-1)*(TextH+SpaceH) EditW EditH],...
                                                                        'TooltipString',['Default: ' AttributeOpt]);
                        eval(['handles.Attrh' num2str(AttrhNumber) '.' Attributes(ind).Name '=Attrh.' Attributes(ind).Name ';']);
                    otherwise
                        Attrh.(Attributes(ind).Name)=uicontrol(parent_handle,'Style', 'edit', 'String', Attributes(ind).Value,'Units','normalized','BackgroundColor',[1 1 1],...
                                                                    'Position', [(j-1)*(TextW+SpaceW+EditW)+TextW (i-1)*(TextH+SpaceH) EditW EditH],...
                                                                    'TooltipString',['Default: ' Attributes(ind).Value]);
                        eval(['handles.Attrh' num2str(AttrhNumber) '.' Attributes(ind).Name '=Attrh.' Attributes(ind).Name ';']);
                end
            else
                Attrh.(Attributes(ind).Name)=uicontrol(parent_handle,'Style', 'edit', 'String', Attributes(ind).Value,'Units','normalized','BackgroundColor',[1 1 1],...
                                                                    'Position', [(j-1)*(TextW+SpaceW+EditW)+TextW (i-1)*(TextH+SpaceH) EditW EditH],...
                                                                    'TooltipString',['Default: ' Attributes(ind).Value]);
                eval(['handles.Attrh' num2str(AttrhNumber) '.' Attributes(ind).Name '=Attrh.' Attributes(ind).Name ';']);
            end
        end
        ind=ind+1;
    end
end
guidata(parent_handle,handles);
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function ok = t_ok(cond, msg)
%T_OK  Tests whether a condition is true.
%   T_OK(EXPR, MSG) increments the global test count and if the EXPR
%   is true it increments the passed tests count, otherwise increments
%   the failed tests count. Prints 'ok' or 'not ok' followed by the
%   MSG, unless T_BEGIN was called with input QUIET equal true. Intended
%   to be called between calls to T_BEGIN and T_END.
%
%   Optionally returns a true or false value indicating whether or
%   not the test succeeded.
%
%   Example:
%       quiet = 0;
%       t_begin(5, quiet);
%       t_ok(pi > 3, 'size of pi');
%       t_skip(3, 'not yet written');
%       t_is(2+2, 4, 12, '2+2 still equals 4');
%       t_end;
%
%   See also T_IS, T_SKIP, T_BEGIN, T_END, T_RUN_TESTS.

%   MP-Test
%   Copyright (c) 2004-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Test.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mptest for more info.

global t_quiet;
global t_counter;
global t_ok_cnt;
global t_not_ok_cnt;

if nargin < 2 || strcmp(msg, '')
    msg = '';
else
    msg = [' - ', msg];
end
if cond
    t_ok_cnt = t_ok_cnt + 1;
else
    t_not_ok_cnt = t_not_ok_cnt + 1;
    if ~t_quiet
        fprintf('not ');
    end
end
if ~t_quiet
    fprintf('ok %d%s\n', t_counter, msg);
end
t_counter = t_counter + 1;
if nargout
    ok = cond;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
samplen=1000;%Total number of initial states
InitialMat=zeros(4,samplen);%Store initialized x, y, theta, v into each row by sequence
for i=1:samplen
rr=0;
while rr<=1
xini1=(-10+rand(1)*(20));%-10~10
yini1=(-10+rand(1)*(20));%-10~10
rr = (xini1)^2+(yini1)^2;%Initialize location states which are inside safe region
end
thetaini1=(-10+rand(1)*(20));%-10~10
vini1=(-10+rand(1)*(20));%-10~10
InitialMat(1,i)=xini1;
InitialMat(2,i)=yini1;
InitialMat(3,i)=thetaini1;
InitialMat(4,i)=vini1;
end
save('InitialStateData','InitialMat');

author: HybridRobotics

repo: NMPC-DCLF-DCBF


text:
function varargout = plotlattice(varargin)
%PLOTLATTICE  Plots an integer lattice
%
% p = plotlattice(C,which,c,size,options)
%
% Note that only convex sets C in R^2 are supported.
%
% C    :  Constraint object
% which:  'inner' or 'outer'
% color:  color [double] ([r g b] format) or char from 'rymcgbk'
% size :  Size of marker
% options: options structure from sdpsettings

varargin{1} = lmi(varargin{1});
if nargout == 0
     plotlattice(varargin{:});
else
    varargout{1} = plotlattice(varargin{:});
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% op_CSIFOurierTransform.m
% Brenden Kadota, SunnyBrook Hosptial 2021.
%
% Does spectral and spatial fourier transform on MRSI structure. The fast
% fourier transform is done on the spectral dimension. If the spatial k
% space is cartesian, the fast fourier transform is applied.
% USAGE:
% out = op_CSIFourierTransform(out, spatialFT, spectralFT, isCartesian)
%
%
% input:    in         = Twix object of CSI data
%           spatialFT  = Boolean flag (1 or 0) to compute fourier
%                        transfrom along spatial dimension
%           spectralFT = Same as above but for the spectral dimension
%           isCartesian= boolean flag for fft or slow fourier tranform alon
%                        the spatial dimension
%
% output:   in         = Twix object with new out.specs field of the
%                         fourier transformed data
function MRSIStruct = op_CSIFourierTransform(MRSIStruct, k_file, fourierTransform)
    arguments
        MRSIStruct (1,1) struct
        k_file (1,:) char {mustBeFileorDefault} = ""
        fourierTransform.spatial (1,1) logical {mustHaveSpatial(fourierTransform.spatial, MRSIStruct)}
        fourierTransform.spectral (1,1)  logical {mustHaveSpectral(fourierTransform.spectral, MRSIStruct)}
    end
    % set default values for spatial and spectral fourier transform
    fourierTransform = setDefaultFlags(fourierTransform, MRSIStruct);

    % spatial dimension fourier transform
    if (fourierTransform.spatial == 1)
        disp('Calculating spatial dimension');
        if(k_file == "")
            %applying the fast fourier transform if k space is cartesian
            MRSIStruct = applyFastFourierTransformSpatial(MRSIStruct);
        else
            [kTrajectory, numSpatial, numSpectral] = readKSpaceFile(k_file, MRSIStruct);
            MRSIStruct = slowFourierTransfrom(MRSIStruct, kTrajectory, numSpatial, numSpectral);
            MRSIStruct = calculateSpectralValues(MRSIStruct, numSpatial, numSpectral);
        end

        MRSIStruct = setFlags(MRSIStruct, 'spatialFT', true);
        %waterRemoved = op_CSIRemoveWater(MRSIStruct);
    end

    % spectral fourier transform
    if (fourierTransform.spectral)
        % fast fourier transform along the time dimension
        disp('Calculating spectral dimension');
        MRSIStruct = fastFourierTransformTime(MRSIStruct);
    end

end

function sft2_Oper = sft2_Operator(InTraj, OutTraj, Ift_flag)
    %
    % sft2 Perform 2-dimensional slow Fourier transform
    %
    % This function was written by Bernhard Strasser, April 2018.
    %
    %
    % The function performs a slow Fourier transform in two spatial dimensions by using the definition of the Fourier Transform
    % FT(f)(a_j) = sum(k=-N/2;N/2-1) f(x_k)exp(-2pi i x_k a_j)
    % (For even N. For odd, change sum slightly)
    %
    %
    % sft2Operator = sft2(A,Ifft_flag)
    %
    % Input:
    % -         Ift_flag                    ...     Flag of 0 or 1 determining if an inverse Fourier transform or normal should be performed.
    % -         InputSize                   ...     Size of Input data.
    %
    % Output:
    % -         sft2Operator                ...     Fourier transformation matrix which can be applied to the data by Out = sft2Operator*In
    %
    %
    % Feel free to change/reuse/copy the function.
    % If you want to create new versions, don't degrade the options of the function, unless you think the kicked out option is totally useless.
    % Easier ways to achieve the same result & improvement of the program or the programming style are always welcome!
    % File dependancy: ?
    % 0. Preparations

    % Define Exponent of exp
    if(~Ift_flag)
        Expy = -2*pi*1i;
    else
        Expy = 2*pi*1i;
    end

    % Define Output Size
    NOut = size(OutTraj,1);

    % 1. FT

    sft2_Oper = zeros([size(OutTraj,1) size(InTraj,1)]);
    for j=1:NOut
        sft2_Oper(j,:) = exp(Expy*(OutTraj(j,1)*InTraj(:,1) ...
            + OutTraj(j,2)*InTraj(:,2)));
    end

    if(Ift_flag)
        sft2_Oper = sft2_Oper / (size(InTraj,1));
    end

end

function mustHaveSpatial(a, in)
    if isfield(in, 'spatialFT') && (a == true && in.spatialFT == 1)
        eidType = 'spatialTransform:noSpatialTransform';
        msgType = 'Spaitial Fourier Transform already done!';
        throwAsCaller(MException(eidType,msgType))
    end
end

function mustHaveSpectral(a, in)
    if isfield(in, 'spectral') && (a == true && in.spectral == 1)
        eidType = 'spectralTransform:noSpectralTransform';
        msgType = 'Spectral Fourier Transform already done!';
        throwAsCaller(MException(eidType,msgType))
    end
end

function mustBeFileorDefault(file)
    if ~isfile(file) && ~strcmp(file, "")
        eidType = 'mustBeFile:notAFile';
        msgType = 'Invalid value fo k_file, must be a file';
        throwAsCaller(MException(eidType,msgType))
    end
end

function fourierTransform = setDefaultFlags(fourierTransform, in)
    if(~isfield(fourierTransform, 'spatial'))
        if(in.flags.spatialFT)
            fourierTransform.spatial = 0;
        else
            fourierTransform.spatial = 1;
        end
    end
    if(~isfield(fourierTransform, 'spectral'))
        if(in.flags.spectralFT)
            fourierTransform.spectral = 0;
        else
            fourierTransform.spectral = 1;
        end
    end
end

function MRSIStruct = applyFastFourierTransformSpatial(MRSIStruct)
    disp('Applying fast fourier transform');
    %apply half pixel shift
    MRSIStruct = halfPixelShift(MRSIStruct);
    %get data
    data = getData(MRSIStruct);
    xDimension = getDimension(MRSIStruct, 'kx');
    if(mod(getSizeFromDimensions(MRSIStruct, {'kx'}), 2) == 1)
        data = circshift(data, 1, xDimension);
    end
    data = fftshift(fft(fftshift(data, xDimension), [], xDimension), xDimension);

    yDimension = getDimension(MRSIStruct, 'ky');
    if(mod(getSizeFromDimensions(MRSIStruct,{'ky'}), 2) == 1)
        data = circshift(data, 1, yDimension);
    end
    data = fftshift(fft(fftshift(data, yDimension), [], yDimension), yDimension);
    MRSIStruct = setData(MRSIStruct, data);

    MRSIStruct = setDimension(MRSIStruct, 'x', getDimension(MRSIStruct, 'kx'));
    MRSIStruct = setDimension(MRSIStruct, 'y', getDimension(MRSIStruct, 'ky'));
    MRSIStruct = setDimension(MRSIStruct, 'z', getDimension(MRSIStruct, 'kz'));
    MRSIStruct = setDimension(MRSIStruct, 'kx', 0);
    MRSIStruct = setDimension(MRSIStruct, 'ky', 0);
    MRSIStruct = setDimension(MRSIStruct, 'kz', 0);
end

function MRSIStruct = slowFourierTransfrom(MRSIStruct, kTrajectory, numSpatial, numSpectral)

    [xCoordinates, yCoordinates, imageTrajectory] = getImageTrajectory(MRSIStruct);

    %creating fourier transform operator for spatial domain
    sftOperator = sft2_Operator(kTrajectory, imageTrajectory, 0);

    %permute so first 3 dimensions are x, y and t
    [MRSIStruct, prevPermute, prevSize] = reshapeDimensions(MRSIStruct, {'t', 'ky'});
    data = getData(MRSIStruct);

    % apply slow fourier transform matrix to data
    image = applySlowFourierTranformMatrix(MRSIStruct, sftOperator, data, numSpectral, numSpatial);
    MRSIStruct = setData(MRSIStruct, image);

    kyDimension = getDimension(MRSIStruct, 'ky');
    prevPermute = removeDimPrevPermute(prevPermute, kyDimension);
    prevPermute = addDimPrevPermute(prevPermute, 'y', kyDimension);
    prevPermute = addDimPrevPermute(prevPermute, 'x', kyDimension + 1);

    prevSize(1) = numSpectral;
    prevSize(2) = length(yCoordinates);
    prevSize = [prevSize(1:2), length(xCoordinates), prevSize(3:end)];
    MRSIStruct = reshapeBack(MRSIStruct, prevPermute, prevSize);
end

function image = applySlowFourierTranformMatrix(MRSIStruct, sftOperator, data, numSpectral, numSpatial)
    yLength = length(getCoordinates(MRSIStruct, 'y'));
    xLength = length(getCoordinates(MRSIStruct, 'x'));

    %image dimensions after fourier tranform is time, y, x, and extras.
    imageDimensions = [numSpectral, yLength, xLength, ...
                                        getSizeFromDimensions(MRSIStruct, {'extras'})];
    image = zeros(imageDimensions);

    for iPoint = 1:numSpectral
        startingPoint = (iPoint - 1)*numSpatial + 1;
        endingPoint = iPoint * numSpatial;
        kSpaceSlice = data(startingPoint:endingPoint, :, :);
        vectorizedSlice = reshape(kSpaceSlice, [], size(kSpaceSlice,3));
        ftVectorizedSlice = sftOperator*vectorizedSlice;

        imageSlice = reshape(ftVectorizedSlice, [yLength, xLength, size(ftVectorizedSlice,2)]);
        image(iPoint, :, :, :) = imageSlice;
    end
end

function [xCoordinates, yCoordinates, imageTrajectory] = getImageTrajectory(MRSIStruct)
    xCoordinates = getCoordinates(MRSIStruct, 'x');
    yCoordinates = getCoordinates(MRSIStruct, 'y');

    %applying the slow fourier transform if the k space is non cartesian
    [x, y] = meshgrid(xCoordinates, yCoordinates);
    imageTrajectory = [x(:), y(:)];
end

function MRSIStruct = fastFourierTransformTime(MRSIStruct)
    data = getData(MRSIStruct);
    timeDimension = getDimension(MRSIStruct, 't');
    %fourier transform in the spectral domain
    data = fftshift(fft(data, [], timeDimension), timeDimension);

    MRSIStruct = setData(MRSIStruct, data);
    ppm = calculatePPM(MRSIStruct);

    if strcmp(MRSIStruct.nucleus,'1H')
        ppm = ppm + 4.65;
    end

    MRSIStruct = setPPM(MRSIStruct, ppm);
    %flip ppm
    MRSIStruct = setFlags(MRSIStruct, 'spectralFT', true);
end

function ppm = calculatePPM(MRSIStruct)
    %get gamma and nucleus
    gamma = MRSIStruct.gamma;

    %lower bounds of frequency
    spectralWidth = getSpectralWidth(MRSIStruct);
    timeSize = getSizeFromDimensions(MRSIStruct, {'t'});

    step = spectralWidth/timeSize;
    lowerBound = -spectralWidth/2 + step/2;
    %upper bounds of frequency
    upperBound = spectralWidth/2 - step/2;

    %calculating the frequency
    frequency=lowerBound:step:upperBound;

    %calculating the ppm
    ppm=-frequency/(MRSIStruct.Bo*gamma);
end

function MRSIStruct = calculateSpectralValues(MRSIStruct, numSpatial, numSpectral)
    spectralDwellTime = calculateSpectralDwellTime(MRSIStruct, numSpatial);
    spectralWidth = 1/spectralDwellTime;
    spectralTime = calculateSpectralTime(spectralDwellTime, numSpectral);

    MRSIStruct = setSpectralWidth(MRSIStruct, spectralWidth);
    MRSIStruct = setSpectralDwellTime(MRSIStruct, spectralDwellTime);
    MRSIStruct = setSpectralTime(MRSIStruct, spectralTime);
end

function spectralDwellTime = calculateSpectralDwellTime(MRSIStruct, spatialPoints)
    adcDwellTime = getAdcDwellTime(MRSIStruct);
    spectralDwellTime = spatialPoints * adcDwellTime;
end

function spectralTime = calculateSpectralTime(spectralDwellTime, spatialPoints)
    spectralTime = 0:spectralDwellTime:spectralDwellTime*(spatialPoints - 1);
end

function MRSIStruct = halfPixelShift(MRSIStruct)
    kx = getCoordinates(MRSIStruct, 'kx');
    ky = getCoordinates(MRSIStruct, 'ky');
    halfPixelX = getVoxSize(MRSIStruct, 'x')/2;
    halfPixelY = getVoxSize(MRSIStruct, 'y')/2;
    kShift = kx*halfPixelX + ky'*halfPixelY;

    [MRSIStruct, prevPermute, prevSize] = reshapeDimensions(MRSIStruct, {'ky', 'kx'});
    data = getData(MRSIStruct);
    data = data .* exp(-1i*2*pi*kShift);
    MRSIStruct = setData(MRSIStruct, data);
    MRSIStruct = reshapeBack(MRSIStruct, prevPermute, prevSize);
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [ transformedX, transformedY, transformedZ ] = TangentSpaceTransform( x, y, z, meanShape )
%TANGENTSPACETRANSFORM Summary of this function goes here
%   Detailed explanation goes here

    scaling = [ x y z] * [ meanShape(:,1)' meanShape(:,2)' meanShape(:,3)']';
    for i=1:size(x,1)
        x(i,:) = x(i,:) * (1 / scaling(i));
        y(i,:) = y(i,:) * (1 / scaling(i));
        z(i,:) = z(i,:) * (1 / scaling(i));
    end

    transformedX = x * mean(scaling);
    transformedY = y * mean(scaling);
    transformedZ = z * mean(scaling);

end

repo: OpenFace


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

%%
clear; close all; clc;

%% Syntax
% |[varargout]=dicom3Dpar(D);|

%% Description
% This function retrieves the voxel size, position, and orientation
% information from a DICOM file.

%% Examples
%

%%
% Path name for dicom files

defaultFolder = fileparts(fileparts(mfilename('fullpath'))); %Set main folder
pathName=fullfile(defaultFolder,'data','DICOM','0001_human_calf');
fileName=fullfile(pathName,'00001.dcm');

%%
% Get dicom information
dcmInfo=dicominfo(fileName);

%%
% Get 3D/geometry information, e.g. voxel size, position, and orientation
% information.

[v,OR,r,c]=dicom3Dpar(dcmInfo)

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%{
    Copyright (C) Cambridge Electronic Design Limited 2014
    Author: James Thompson
    Web: www.ced.co.uk email: james@ced.co.uk, softhelp@ced.co.uk

    This file is part of CEDS64ML, a MATLAB interface to the SON64 library.

    CEDS64ML 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.

    CEDS64ML 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 CEDS64ML.  If not, see <http://www.gnu.org/licenses/>.
%}

function [ iRead, fVals, i64Time ] = CEDS64ReadWaveF( fhand, iChan, iN, i64From, i64To, maskh )
%CEDS64READWAVE32 Reads wave data from a waveform or realwave channel as singles (32-bit floats).
%   [ iRead, fVals, i64Time ] = CEDS64ReadWaveF( fhand, iChan, iN, i64From {, i64To {, maskh} } )
%   Inputs
%   fhand - An integer handle to an open file
%   iChan - A channel number for a Waveform or Realwave channel
%   iN - The maximum number of data points  to copy
%   i64From - The time in ticks of the earliest time you want to read
%   i64To - (Optional) The time in ticks of the latest time you want to
%   read. If not set or set to -1, read to the end of the channel
%   maskh - (Optional) An integer handle to a marker mask (only used when reading wavemarkers)
%   Outputs
%   iRead - The number of data points read
%   fVals - An array of floats conatining the data points
%   i64Time - The time in ticks of the first data point

if (nargin < 4)
    iRead = -22;
    return;
end

outwavepointer = zeros(iN,1, 'single');
outtimepointer = zeros(1,1,'int64');

if (nargin < 5)
    i64To = -1;
end

if (nargin < 6)
    maskcode = -1;
else
    maskcode = maskh;
end

[iRead, fVals, i64Time] = calllib('ceds64int', 'S64ReadWaveF', fhand, iChan, outwavepointer, iN, i64From, i64To, outtimepointer, maskcode);

if iRead > 0
    fVals(iRead+1:end) = [];
else
    fVals = [];
end
end

author: brainstorm-tools


text:
function [y,nz] = ompdenoise1(params,msgdelta)
%OMPDENOISE1 OMP denoising of 1-D signals.
%  OMPDENOISE1 denoises a 1-dimensional signal using OMP denoising. The
%  function syntax is identical to OMPDENOISE, but it runs significantly
%  faster on 1-D signals. OMPDENOISE1 requires somewhat more memory than
%  OMPDENOISE (approximately the size of the input signal), so if memory is
%  limited, OMPDENOISE can be used instead.
%
%  See also OMPDENOISE.

%  Ron Rubinstein
%  Computer Science Department
%  Technion, Haifa 32000 Israel
%  ronrubin@cs
%
%  August 2009

% parse input arguments %

x = params.x(:);
D = params.dict;
blocksize = params.blocksize;

% maxval %
if (isfield(params,'maxval'))
  maxval = params.maxval;
else
  maxval = 1;
end

% gain %
if (isfield(params,'gain'))
  gain = params.gain;
else
  gain = 1.15;
end

% maxatoms %
if (isfield(params,'maxatoms'))
  maxatoms = params.maxatoms;
else
  maxatoms = floor(blocksize/2);
end

% stepsize %
if (isfield(params,'stepsize'))
  stepsize = params.stepsize;
else
  stepsize = 1;
end
if (any(stepsize<1))
  error('Invalid step size.');
end

% noise mode %
if (isfield(params,'noisemode'))
  switch lower(params.noisemode)
    case 'psnr'
      sigma = maxval / 10^(params.psnr/20);
    case 'sigma'
      sigma = params.sigma;
    otherwise
      error('Invalid noise mode specified');
  end
elseif (isfield(params,'sigma'))
  sigma = params.sigma;
elseif (isfield(params,'psnr'))
  sigma = maxval / 10^(params.psnr/20);
else
  error('Noise strength not specified');
end

% lambda %
if (isfield(params,'lambda'))
  lambda = params.lambda;
else
  lambda = maxval/(10*sigma);
end

% msgdelta %
if (nargin <2)
  msgdelta = 5;
end
if (msgdelta<=0)
  msgdelta = -1;
end

epsilon = sqrt(blocksize) * sigma * gain;   % target error for omp

MEM_LOW = 1;
MEM_NORMAL = 2;
MEM_HIGH = 3;

if (isfield(params,'memusage'))
  switch lower(params.memusage)
    case 'low'
      memusage = MEM_LOW;
    case 'normal'
      memusage = MEM_NORMAL;
    case 'high'
      memusage = MEM_HIGH;
    otherwise
      error('Invalid memory usage mode');
  end
else
  memusage = MEM_NORMAL;
end

% compute G %

G = [];
if (memusage >= MEM_NORMAL)
  G = D'*D;
end

% verify dictionary normalization %

if (isempty(G))
  atomnorms = sum(D.*D);
else
  atomnorms = diag(G);
end
if (any(abs(atomnorms-1) > 1e-2))
  error('Dictionary columns must be normalized to unit length');
end

% denoise the signal %

% process the signal in batches to conserve memory
% choose batchsize so im2colstep returns a matrix of approximately the same
% size as the signal
batchsize = ceil(length(x)*stepsize/blocksize + blocksize);

y = zeros(size(x));
ids = 1:min(batchsize,length(x));
nz = 0;

blocknum = floor((length(x)-blocksize)/stepsize) + 1;
processedblocks = 0;
tid = timerinit('ompdenoise', blocknum);
while (length(ids)>=blocksize)

  % extract the signal blocks
  blocks = im2colstep(x(ids),[blocksize 1],[stepsize 1]);

  % remove DC
  [blocks, dc] = remove_dc(blocks,'columns');

  % denoise the blocks
  if (memusage == MEM_LOW)
    gamma = omp2(D,blocks,[],epsilon,'maxatoms',maxatoms,'checkdict','off');
  else
    gamma = omp2(D'*blocks,sum(blocks.*blocks),G,epsilon,'maxatoms',maxatoms,'checkdict','off');
  end
  nz = nz + nnz(gamma);
  cleanblocks = add_dc(D*gamma, dc, 'columns');

  y(ids) = y(ids) + col2imstep(cleanblocks, [length(ids) 1], [blocksize 1], [stepsize 1]);
  ids = ids + floor((batchsize-blocksize)/stepsize)*stepsize + stepsize;
  if (ids(end)>length(x))
    ids = ids(ids<=length(x));
  end

  % display status
  if (msgdelta>0)
    processedblocks = processedblocks + size(blocks,2);
    timereta(tid, processedblocks, msgdelta);
  end

end

if (msgdelta>0)
  timereta(tid, blocknum);
end

cnt = countcover(size(x),[blocksize 1],[stepsize 1]);
y = (y+lambda*x)./(cnt + lambda);
y = reshape(y,size(params.x));

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function out = EN_SampEn(y,M,r,preProcessHow)
% EN_SampEn     Sample Entropy of a time series
%
% SampEn(m,r), using code from PhysioNet.
%
% Uses a compiled C version of the code if available, otherwise uses a (slower)
% Matlab implementation (which can actually be faster for shorter time series
% due to overheads of reading/writing to disk)
%
% The publicly-available PhysioNet Matlab code, sampenc (renamed here to
% RN_sampenc) is available from:
% http://www.physionet.org/physiotools/sampen/matlab/1.1/sampenc.m
%
% cf. "Physiological time-series analysis using approximate entropy and sample
% entropy", J. S. Richman and J. R. Moorman, Am. J. Physiol. Heart Circ.
% Physiol., 278(6) H2039 (2000)
%
% This function can also calculate the SampEn of successive increments of time
% series, i.e., we using an incremental differencing pre-processing, as
% used in the so-called Control Entropy quantity:
%
% "Control Entropy: A complexity measure for nonstationary signals"
% E. M. Bollt and J. Skufca, Math. Biosci. Eng., 6(1) 1 (2009)
%
%---INPUTS:
% y, the input time series
% M, the embedding dimension
% r, the threshold
% preProcessHow [opt], (i) 'diff1', incremental differencing (as per 'Control Entropy').

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

% Check y is a column vector:
if isempty(y)
    error('Must input a valid data vector');
end
if size(y,1)==1
    y = y';
end

% Embedding dimension:
if nargin < 2
    M = 2;
end

% Tolerance:
if nargin < 3
    r = 0.1*std(y);
end

if nargin < 4
    preProcessHow = ''; % don't apply any preprocessing
end

%-------------------------------------------------------------------------------
% Can specify to first apply an incremental differencing of the time series
% thus yielding the 'Control Entropy':
% "Control Entropy: A complexity measure for nonstationary signals"
% E. M. Bollt and J. Skufca, Math. Biosci. Eng., 6(1) 1 (2009)
if ~isempty(preProcessHow)
    y = BF_PreProcess(y,preProcessHow);
end

% ------------------------------------------------------------------------------
% Use the physionet code to calculate the Sample Entropy using these parameters:
% ------------------------------------------------------------------------------
% Check if a compiled C version exists:
% [a,b] = system('which sampen');

sampEn = sampen_mex(y',M+1,r);
sampEn = sampEn(1:M+1); % always that extra one for the M=0

%-------------------------------------------------------------------------------
% This is dangerous because it could lead to inconsistent implementations across runs
%-------------------------------------------------------------------------------
% warning('No mex file found: using a slower native Matlab implementation instead');
% No mex version available; use (much slower) Matlab implementation
% if isempty(b) || length(y) < 3000 % faster to run within Matlab
% sampEn = PN_sampenc(y,M+1,r,false);
%     fprintf('Using compiled C code~~~\n')
%     % http://www.physionet.org/physiotools/sampen/c/
%     % (use Makefile in Toolboxes/Physionet/ to run make, then make install)
%
%     % Run compiled C code:
%     filePath = BF_WriteTempFile(y);
%     command = sprintf('sampen -m %u -r %f < %s',M,r,filePath);
%     [~,res] = system(command);
%     fprintf(1,'%s\n',res)
%     s = textscan(res,'%[^\n]'); s = s{1};
%     sampEn = zeros(M,1);
%     for i = 1:M
%         [~,params] = regexp(s{i},'\((\S+)\)','tokens','match');
%         params = regexp(params{1}(2:end-1),',','split');
%         [~,result] = regexp(s{i},'= (\S+)','tokens','match');
%         result = str2num(result{1}(3:end));
%         sampEn(i) = result;
%     end
% end
% end

% ------------------------------------------------------------------------------
% Compute outputs from the code
% ------------------------------------------------------------------------------
for i = 1:M+1
    % Sample entropy:
    out.(sprintf('sampen%u',i-1)) = sampEn(i);

    % Quadratic sample entropy (QSE), Lake (2006):
    % (allows better comparison across r values)
    out.(sprintf('quadSampEn%u',i-1)) = sampEn(i) + log(2*r);

    % COSEn (Lake and Moorman, 2011), doesn't really make sense in general;
    % especially for z-scored series!:
    % out.(sprintf('COSEn%u',i)) = sampEn(i) + log(2*r) - log(mean(y));
end

if M > 1
    out.meanchsampen = mean(diff(sampEn));
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function pass = test_nonlinSysDampingBreaks_US(pref)
% Test 2x2 system (sin/cos) where damping is required
%
% Asgeir Birkisson, May 2014.

if ( nargin == 0 )
    pref = cheboppref;
end

tol = 1e-10;

%% Piecewise:
d = [-pi 0 pi];
x = chebfun('x',d);
f = [ 0*x ; 0*x ];

%% ULTRAS:
pref.discretization = @ultraS;

A = chebop(@(x,u,v) [u - diff(v,2); diff(u,2) + cos(v)], d);
A.lbc = @(u,v) [u-1/2; v+1/4];
A.rbc = @(u,v) [v-1/4; u+1/2];

[u1, u2, info] = solvebvp(A, f, pref);

% Want to check BCs as well.
bcFunLeft = chebfun(A.lbc(u1,u2));
bcFunRight = chebfun(A.rbc(u1,u2));

% And check that we're continuous over breakpoint
u1jump = jump(u1, 0);
u2jump = jump(u2, 0);

pass(1) = abs(info.error) < tol;
pass(2) = norm(bcFunLeft(d(1))) < tol && norm(bcFunRight(d(end))) < tol;
pass(3) = norm(u1jump) < tol && norm(u2jump) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text: % This file is here only so that you can put this directory into the dependencies folder of a BCILAB installation.

author: PatternRecognition

repo: OpenBMI