text:
%PRWAITBAR Report PRTools progress by single waitbar
%
%   H = PRWAITBAR(N,M,TEXT)
%   H = PRWAITBAR(N,TEXT,FLAG)
%   S = PRWAITBAR
%
% INPUT
%   N      Integer, total number of steps in loop
%   M      Integer, progress in number of steps in loop
%   TEXT   Text to be displayed in waitbar
%   FLAG   Flag (0/1)
%
% OUTPUT
%   H      Waitbar handle
%   S      Status PRWAITBAR ('on' or 'off')
%
% DESCRIPTION
% This routine may be used to report progress in PRTools experiments.
% It detects and integrates levels of loops. The following calls are
% supported:
%
%   PRWAITBAR(N,TEXT)        initialize loop
%   PRWAITBAR(N,TEXT,FLAG)   initialize loop if FLAG == 1
%   PRWAITBAR(N,M)           update progress
%   PRWAITBAR(N,M,TEXT)      update info
%   PRWAITBAR(0)             closes loop level
%   PRWAITBAR OFF            removes waitbar
%   PRWAITBAR ON             switches waitbar on again
%   PRWAITBAR REPORT         no waitbar, report progress in command window
%   PRWAITBAR                reset prwaitbar
%
% A typical sequence of calls is:
%    .....
%    prwaitbar(nfolds,'cross validation')
%    for j = 1:nfolds
%        prwaitbar(nfolds,j,['cross validation, fold ' int2str(j)])
%        ....
%    end
%    prwaitbar(0)
%    .....
% Calls to PRTWAITBAR may be nested and all progress is merged into a single
% waitbar. In this PRWAITBAR differs from Matlab's WAITBAR.
% A typical example can be visualised by:
%
%    prcrossval({gendatb,gendath},{svc,loglc,fisherc},5,2)
%
% PRWAITBAR may increase the computation time with 25%. For long running
% processes PRWAITBAR REPORT (e.g. in a screen window) is recommended.
%
% Some other, more high-level routines calling PRWAITBAR are
% PRWAITBARINIT, PRWAITBARNEXT, PRWAITBARONCE, PREIG, PRINV, PRPINV, PRRANK,
% PRSVD, PRCOV
%
% SEE ALSO (PRTools Guide)
% PRWAITBARINIT, PRWAITBARNEXT, PRWAITBARONCE

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com

function varargout = prwaitbar(n,m,text)

persistent N                % array, N(DEPTH) is the loop size at level DEPTH
persistent M                % array, M(DEPTH) is the present counter values at level DEPTH
persistent DEPTH            % nesting DEPTH
persistent WHANDLE          % handle of the waitbar
persistent MESS             % cell array with message to be displayed at level DEPTH
persistent STAT             % flag, if STAT = 0, PRWAITBAR is switched off
persistent WPOS             % position of the waitbar on the screen
persistent REP              % flag, if REP = 1, PRWAITBAR reports in the command window
persistent FNAME            % cell array with names of the calling routine at level DEPTH
persistent DOIT             % flag that may switch off tracking progress at deeper levels
persistent NSKIP            % level at which tracking progress has been switched off
persistent OLDPROG
persistent OLDPROG1         % previous progress
persistent OLDPROG2         % previous progress
persistent OLDPROG3         % previous progress for report

% disp(['prwaitbar ' callername])
% disp(nargin)
% if nargin > 0, disp(n),   end
% if nargin > 1, disp(m),   end
% if nargin > 2, disp(text),end
% disp('-------------')
if isempty(DEPTH)           % initialisation of persistent variables
	DEPTH = 0;
	MESS = cell(1,10);
	FNAME = cell(1,10);
	STAT = 1;
	DOIT = 1;
  OLDPROG = 0;
end

decrease_depth = 0;
if nargin == 0                            % ---- call: prwaitbar ----
	if nargout == 0
		%prwaitbar off
		DEPTH = 0;
		STAT = 1;
    REP = 0;
		DOIT = 1;
		MESS = cell(1,10);
		OLDPROG1 = 0;
		OLDPROG2 = 0;
		OLDPROG3 = 0;
		err.message = '##';
		lasterror(err);
% 		if ishandle(WHANDLE)
% 			delete(WHANDLE);
% 		end
	else                                    % ---- call: out = prwaitbar ----
		if nargout == 1
			if STAT == 0, varargout = {'off'};
			else varargout = {'on'};
			end
		else
			varargout = {WHANDLE,N,M,DEPTH,MESS,STAT,REP,WPOS,FNAME,DOIT,NSKIP};
		end
		return
	end

elseif nargin == 1

	if isstr(n) & strcmp(n,'off')           % ---- call: prwaitbar off ----

		if ~isempty(WHANDLE) & ishandle(WHANDLE)
			WPOS = get(WHANDLE,'position');     % save position as user likes it
			close(WHANDLE);
		end
		N = []; M = []; DEPTH = 0;            % set variables at initial values
		WHANDLE = [];
		MESS = cell(1,10);
		STAT = 0;
 		REP = 0;
		return

	elseif isstr(n) & strcmp(n,'on')        % ---- call: prwaitbar on ----
		STAT = 1;                            % restart waitbar
		REP = 0;
    DEPTH = [];
		return

  elseif isstr(n) & strcmp(n,'report')    % ---- call: prwaitbar report ----
		prwaitbar off;
    err.message = '##';
    lasterror(err);
		STAT = 1;
		REP = 1;
		DEPTH = 0;
		if ishandle(WHANDLE)
			delete(WHANDLE);
		end
		WHANDLE = [];
    OLDPROG3 = 0;
		return

	elseif n == 0 & STAT                    % ---- call: prwaitbar(0) ----
		if ~DOIT
			if (NSKIP == DEPTH)                 % we are back at the right level, so
				DOIT = 1;                         % restart tracking progress
			end
			DEPTH = max(DEPTH - 1,1);           % counting problem !!! caused by catch-try loops
      return                              % let us continue with DEPTH = 1
		end
		if DEPTH < 1                        % wrong call, return
			return
		end
		M(DEPTH) = N(DEPTH);                % loop almost ready
		decrease_depth = 1;                 % decrease depth after plotting progress
	else
		return
	end

elseif ~STAT                            % next commands can be skipped if switched off
	return

elseif nargin == 2

	if isstr(m)                           % ---- call: prwaitbar(n,text) ----
			                                  % note: this is a typical call to
			                                  % start a waitbar cycle

		DEPTH = DEPTH + 1;                  % we nest one level deeper
		if ~DOIT                            % skip if we dont track progress
			return
		end

		err = lasterror;
		% we have to find out whether we continue a proper set of prwaitbar calls, or whether
		% we have to recover from some error or possible interrupt and have to re-initialise
		% the admin

    if (isempty(err.message)) | (~ishandle(WHANDLE)) | ...
            (~strcmp(err.message,'##') &  ...
            isempty(strfind(err.identifier,'prtools:')))
        % yes, restart after interrupt
      if REP
        prwaitbar report;
      else
        prwaitbar;
      end
      DEPTH = 1;
%     elseif strcmp(err.message,'')
%       err.message = '##';
%       lasterror(err);
    end
			                                   % find the name of the calling routine
		[ss,ii] = dbstack;
		if length(ss) == 1
			name = '';                         % call from command line!!?? debugging?
		else
			[path,name] = fileparts(ss(2).name);
			if strcmp(name,mfilename) | strncmp(name,'prwait',6) % search deeper in case of internal call
				[path,name] = fileparts(ss(3).name);
			end
		end
		FNAME{DEPTH} = name;                 % store it, and display it in waitbar
		if ~REP
			set(WHANDLE,'name',['                         PRWAITBAR:   ' name]);
		end

		N(DEPTH) = n;                        % complete the admin
		M(DEPTH) = 0;
		MESS{DEPTH} = m;                     % note: m contains text
    if DEPTH == 1
      OLDPROG = 0;
    end

	elseif DOIT                            % ---- call: prwaitbar(n,m) ----
		DEPTH = max(DEPTH,1);                % store admin
		N(DEPTH) = n;
		M(DEPTH) = m-1;                      % store m-1, assuming call is in the start of loop
	else                                   % skip and return if progress should not be reported
		return
	end

elseif nargin == 3

	if isstr(m)                            % ---- call: prwaitbat(n,text,flag) ----
		if ~DOIT                             % there is a flag, but we already skip progress
			DEPTH = DEPTH + 1;                 % just keep track of loop nesting.
		else                                 % we are in reporting mode
			DOIT = (text > 0);                 % check flag
			if ~DOIT                           % stop progress reporting!!!
				DEPTH = DEPTH + 1;               % keep track of loop nesting
				NSKIP = DEPTH;                   % store level
			else                               % flag does not apply,
				prwaitbar(n,m);                  % just print waitbar
			end
		end
		return

	else                                   % ---- call: prwaitbar(n,m,text) ----

		if ~DOIT                             % skip progress report if needed
			return
		end

		if isempty(WHANDLE) & ~REP           % switch waitbar on if needed
			prwaitbar;                         % this only happens if the loop is not properly opened
    end
    if DEPTH > 0
      N(DEPTH) = n;                      % update admin
      M(DEPTH) = m-1;
      MESS{DEPTH} = text;
    end
	end
end

if ~STAT | ~DOIT                         % dont display, (just to be sure)
	return
end

% waitbar creation and display

if isempty(WHANDLE) & ~REP               % if no waitbar, make one

	s = sprintf(' \n \n \n' );             % room for 4 levels of text

	if ~isempty(WPOS)
		h = waitbar(0,'','position',WPOS);
		h = waitbar(0,h,s);
	else
		h = waitbar(0,s);                    % create waitbar
		WPOS = get(h,'position');
	end
	OLDPROG1 = 0;
	OLDPROG2 = 0;
	if DEPTH == 0, fnam = [];
	else, fnam = FNAME{DEPTH}; end         % retrieve name of calling routine
 	set(h,'name',['                         PRWAITBAR:   ' fnam])
	%set(h,'HandleVisibility','on');
	WHANDLE = h;

	err.message = '##';                    % store for proper continuation in next call
	lasterror(err);

end

progress = 0;                            % progress
if DEPTH == 0                            % still on level zero
	s = sprintf(' \n \n \n' );             % no message
else
	for j=DEPTH:-1:1                       % run over all levels for progress
		progress = (progress + M(j))/(N(j)); % compute total progress
	end
	s = wtext(MESS,N,DEPTH,REP);           % construct message
  progress = max(progress,OLDPROG);      % avoid progress to decrease
end

significant1 = (progress - OLDPROG1 >= 0.005 | OLDPROG1 == 0 | decrease_depth);
significant2 = (progress - OLDPROG2 >= 0.0005);
if ishandle(WHANDLE)
	if significant1                          % show significant updates only
		waitbar(progress,WHANDLE,s,'position',WPOS);            % update waitbar
		OLDPROG1 = progress;
		OLDPROG2 = progress;
	elseif significant2                       % update text more often
		set(get(get(WHANDLE,'children'),'title'),'string',s); drawnow
		OLDPROG2 = progress;
	end
else                                      % handle was not proper (anymore)
	if REP
    if progress - OLDPROG3 >= 0.001
      fprintf(['\n%6.3f ' s],progress)
      OLDPROG3 = progress;
    end
  else % reconstruct waitbar
    if isempty(WPOS)
      WHANDLE = waitbar(progress,s);
    else
      WHANDLE = waitbar(progress,s,'position',WPOS);
    end
	end
	OLDPROG1 = progress;
	OLDPROG2 = progress;
end
OLDPROG = progress;

if decrease_depth                         % update admin
	DEPTH = DEPTH - 1;
	if DEPTH == 0                           % take care that new waitbar starts fresh
		if ~REP
			WPOS = get(WHANDLE,'position');       % but in old position
      %delete(WHANDLE);
      %we don't delete but preserve the waitbar for future use
      %so the following cleanup is needed
      p=get(get(WHANDLE,'children'),'children');
      if( length(p) > 1 )
          set(p(2),'erasemode','normal');
      end
      waitbar(0,WHANDLE,sprintf(' \n \n \n' ));
      set(WHANDLE,'name',['                         PRWAITBAR:      ']);
      if( length(p) > 1 )
        set(p(2),'erasemode','none');
      end
      prwaitbar;
      lasterror('reset');
      %now we are ready to leave
    end
	else
		if significant1
	 		set(WHANDLE,'name',['                         PRWAITBAR:   ' FNAME{DEPTH}]);
		end
	end
end

if nargout > 0
	if REP
		varargout = 0;
	else
		varargout = {WHANDLE};
	end
end

return

function s = wtext(MESS,N,DEPTH,REP)

	t = cell(1,DEPTH);
	n = 0;
	for j=1:DEPTH
		if N(j) > 1
			n = n+1;
			t{n} = MESS{j};
		end
	end

	if REP
		if n > 0
			s = sprintf('%s\n',t{1:n});
		else
			s = '';
		end
	elseif n == 0
		s = sprintf([' \n \n \n \n ']);
	elseif n == 1
		s = sprintf([' \n \n' t{1} '\n ']);
	elseif n == 2
		s = sprintf([' \n' t{1} '\n' t{2} '\n ']);
	elseif n == 3
		s = sprintf([' \n' t{1} '\n' t{2} '\n' t{3}]);
	else
		s = sprintf([t{n-3} '\n' t{n-2} '\n' t{n-1} '\n' t{n}]);
	end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function test_suite = test_midwt
initTestSuite;

function test_midwt_1D
       x = makesig('LinChirp',8);
       h = daubcqf(4,'min');
       L = 2;
       [y,L] = mdwt(x,h,L);
       [x_new,L] = midwt(y,h,L);
assertVectorsAlmostEqual(x, x_new,'relative',0.0001);

function test_midwt_2D
       load lena512;
       x = lena512;
       h = daubcqf(6);
       [y,L] = mdwt(x,h);
       [x_new,L] = midwt(y,h);
assertEqual(L,9);
assertVectorsAlmostEqual(x, x_new,'relative',0.0001);

author: ricedsp


text:
%% Upload templates and subject data from AFQ to Archiva
%

rd = RdtClient('vistasoft');
rd.credentialsDialog;

%% Templates
cd(vistaRootPath);
cd(fullfile('..','AFQ'))
baseDir = pwd;
fullDir = fullfile(baseDir,'templates');

rd.crp('/AFQ/templates');

a = rd.publishArtifacts(fullDir);

%% List the artifacts it thinks it sent up
for ii=1:length(a)
    disp(a(ii).artifactId)
end

%% List, but this time ask the repository for the list
rd.crp('/AFQ/templates/');
a = rd.listArtifacts;

%% Now the AFQ data sets

% These are in AFQ/data
% So the base directory starts at AFQ
% We will then move into AFQ/data and the subdirectories
cd(vistaRootPath);
cd(fullfile('..','AFQ'))
baseDir = pwd;

% First the afq.mat file
rd.crp('/AFQ/data');
fullDir = fullfile(baseDir,'data');
rd.publishArtifacts(fullDir);

% rd.openBrowser('fancy',true)

%%  Find the directory tree under AFQ/data
cd(vistaRootPath);
cd(fullfile('..','AFQ','data'))
baseDir = pwd;

pNames = dirwalk(baseDir);

nDirs = length(pNames);
for ii=1:nDirs
    fullDir = pNames{ii};
    remotePath = strrep(pNames{ii},'/Users/wandell/Github','');
    rd.crp(remotePath);
    rd.publishArtifacts(fullDir);
end

%%

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function out=mplog2(precision)

if nargin==0
 mp_defaults
 precision=default_precision;
end

out_rval=mpfr_const_log2(precision);

out=class(struct('rval',out_rval,...
                  'ival','0',...
                  'precision',precision),'mp');

author: opencobra

repo: cobratoolbox


text:
function []=plot_sb_baselines(vargin)
% PLOT_SB_BASELINES replaced by SB_BASELINE_PLOT
% Optional an input argument ix can be specified containing a vector of the
% small baseline interferograms to keep which will be plotted in the baseline plot.
%
% Modifications:
% DB 12/2012    fix error for input arguments and add the some syntax
% DB 01/2016    Suppress the warning. Alternative remove the function.

% warning('name of this function changed to sb_baseline_plot')

if nargin>0                     % [DB] changed to nargin > 0, as ix cannot be specified otherwize
    sb_baseline_plot(vargin)
else
    sb_baseline_plot
end

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
%% Load the data into MATLAB from a binary log file
% Usage: >> [datapoints, numpoints] = readdata('datafile.log')
% Header information format:
%           String "MWLOGV##"
%           Time/Date 4 bytes (time())
%           Number of Signals per record Logged 1 bytes (256 max)
%           Data Type of Signals Logged  1 bytes (1-10)
%           Number of bytes per record 2 (65535 max)
% Plot Data Example: plot([1:numpoints], datapoints(1,:), [1:numpoints], datapoints(2,:))
% MathWorks Pilot Engineering 2015
% Steve Kuznicki
function [datapts, numpts] = readdata(dataFile)
%%
datapts = 0;
numpts = 0;

if nargin == 0
    dataFile = 'data.bin';
end

fid = fopen(dataFile, 'r');
% load the header information
hdrToken = fread(fid, 8, 'char');
if strncmp(char(hdrToken),'MWLOGV',6) == true
    logTime = uint32(fread(fid, 1, 'uint32'));
    numflds = double(fread(fid, 1, 'uint8'));
    typefld = uint8(fread(fid, 1, 'uint8'));
    recSize = uint16(fread(fid, 1, 'uint16'));
    fieldTypeStr = get_elem_type(typefld);
    datapts = fread(fid, double([numflds, Inf]), fieldTypeStr);
    fclose(fid);
    numpts = size(datapts,2);
end

end

%% get the element type string
function [dtypeStr] = get_elem_type(dtype)
    switch(dtype)
        case 1
            dtypeStr = 'double';
        case 2
            dtypeStr = 'single';
        case 3
            dtypeStr = 'int32';
        case 4
            dtypeStr = 'uint32';
        case 5
            dtypeStr = 'int16';
        case 6
            dtypeStr = 'uint16';
        case 7
            dtypeStr = 'int8';
        case 8
            dtypeStr = 'uint8';
        case 9
            dtypeStr = 'logical';
        case 10
            dtypeStr = 'embedded.fi';
    end
end

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function value = cmach ( job )

%*****************************************************************************80
%
%% CMACH computes machine parameters for complex arithmetic.
%
%  Discussion:
%
%    Assume the computer has
%
%      B = base of arithmetic;
%      T = number of base B digits;
%      L = smallest possible exponent;
%      U = largest possible exponent;
%
%    then
%
%      EPS = B**(1-T)
%      TINY = 100.0 * B**(-L+T)
%      HUGE = 0.01 * B**(U-T)
%
%    If complex division is done by
%
%      1 / (X+i*Y) = (X-i*Y) / (X**2+Y**2)
%
%    then
%
%      TINY = sqrt ( TINY )
%      HUGE = sqrt ( HUGE )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 June 2009
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for FORTRAN usage,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, pages 308-323, 1979.
%
%  Parameters:
%
%    Input, integer JOB:
%    1, EPS is desired;
%    2, TINY is desired;
%    3, HUGE is desired.
%
%    Output, real VALUE, the requested value.
%
  epsilon = 1.0;

  while ( 1 )

    epsilon = epsilon / 2.0;
    s = 1.0 + epsilon;
    if ( s <= 1.0 )
      break
    end

  end

  epsilon = 2.0 * epsilon;

  s = 1.0;

  while ( 1 )

    tiny = s;
    s = s / 16.0;

    if ( s * 1.0 == 0.0 )
      break
    end

  end

  tiny = ( tiny / eps ) * 100.0;

  s = ( 1 + i ) / ( tiny + tiny * i );

  if ( s ~= 1.0 / tiny )
    tiny = sqrt ( tiny );
  end

  huge = 1.0 / tiny;

  if ( job == 1 )
    value = epsilon;
  elseif ( job == 2 )
    value = tiny;
  elseif ( job == 3 )
    value = huge;
  else
    value = 0.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function tmp = rmSliceGet(model,slice,id)
% rmSliceGet - extract slice from model struct and place it
% into temporary struct, also convert to single here.
%
% tmp = rmSliceGet(model,slice,id);
%
% 2008/01 SOD: extracted from rmGridFit.

if ~exist('model','var') || isempty(model), error('Need model'); end
if ~exist('slice','var') || isempty(slice), slice = 1;           end
if ~exist('id','var') || isempty(id),       id = 1:numel(model); end

% loop over models
tmp = cell(numel(id),1);
for n=id,
    f = {'x0','y0','s','x02','y02','s2','s_major','s_minor','s_theta','rss','rss2','rawrss','rawrss2', 'exponent'};

    % for all models
    tmp{n}.desc = rmGet(model{n},'desc');
    tmp{n}.df   = single(rmGet(model{n},'dfglm'));
    for fn = 1:numel(f),
        val    = rmGet(model{n},f{fn});
        if ~isempty(val),
            % switch on the number of dimensions, since inplane data has
            % diff dimensionality than other views
            switch length(size(val))
                case 3 % presumably means INPLANE view, where slice is the third dimension
                    temp = single(val(:,:,slice));
                    tmp{n}.(f{fn}) = temp(:)';
                otherwise % otherwise slice is the first dimensions
                    tmp{n}.(f{fn}) = single(val(slice,:));
            end
        end;
    end;

    % put all beta values in one matrix
    val      = rmGet(model{n},'b');
    switch length(size(val)) % switch on the number of dimensions, since inplane data has diff dimesnsionality than other views
        case 4 % 4 dimensions means Inplane model:
               %  3 dimensions of coords, and one dimension of beta values.
               %  Hence the dims are x, y, slice, betas. We want to make a
               %  matrix of betas for one slice. This matrix be 2-D, with
               %  one dimension indexing the x and y values (hence
               %  length x * y), and the other dimension representing the
               %  beta values (hence length size(val, 4))
            nbetas = size(val,4);
            nvoxelsPerSlice = size(val,2)*size(val,1);
            tmp{n}.b = zeros(nbetas,nvoxelsPerSlice,'single');
            for fn = 1:nbetas,
                temp  = single(val(:,:,slice,fn));
                tmp{n}.b(fn,:) = temp(:);
            end;
        otherwise
            tmp{n}.b = zeros(size(val,3),size(val,2),'single');
            for fn = 1:size(val,3),
                tmp{n}.b(fn,:) = single(val(slice,:,fn));
            end;
    end
end;

return;


text:
%-----------------------------------------------------------------------
% Job configuration created by cfg_util (rev $Rev: 7703 $)
%-----------------------------------------------------------------------
matlabbatch{1}.spm.meeg.source.headmodel.D = {'C:\Users\gbarnes\Documents\simdata_peakomatic\simdata_1dip_5mm_0pcamp2_amp10\simdata.mat'};
matlabbatch{1}.spm.meeg.source.headmodel.val = 1;
matlabbatch{1}.spm.meeg.source.headmodel.comment = '';
matlabbatch{1}.spm.meeg.source.headmodel.meshing.meshes.template = 1;
matlabbatch{1}.spm.meeg.source.headmodel.meshing.meshres = 2;
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(1).fidname = 'nas';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(1).specification.select = 'nas';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(2).fidname = 'lpa';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(2).specification.select = 'FIL_CTF_L';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(3).fidname = 'rpa';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(3).specification.select = 'FIL_CTF_R';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.useheadshape = 0;
matlabbatch{1}.spm.meeg.source.headmodel.forward.eeg = 'EEG BEM';
matlabbatch{1}.spm.meeg.source.headmodel.forward.meg = 'Single Shell';
matlabbatch{2}.spm.tools.beamforming.data.dir = {'C:\Users\gbarnes\Documents\batchtest\'};
matlabbatch{2}.spm.tools.beamforming.data.D(1) = cfg_dep('M/EEG head model specification: M/EEG dataset(s) with a forward model', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','D'));
matlabbatch{2}.spm.tools.beamforming.data.val = 1;
matlabbatch{2}.spm.tools.beamforming.data.space = 'MNI-aligned';
matlabbatch{3}.spm.tools.beamforming.sources.BF(1) = cfg_dep('Prepare data: BF.mat file', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{3}.spm.tools.beamforming.sources.plugin.grid.resolution = 5;
matlabbatch{3}.spm.tools.beamforming.sources.plugin.grid.space = 'MNI template';
matlabbatch{4}.spm.tools.beamforming.features.BF(1) = cfg_dep('Define sources: BF.mat file', substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{4}.spm.tools.beamforming.features.whatconditions.all = 1;
matlabbatch{4}.spm.tools.beamforming.features.woi = [-Inf Inf];
matlabbatch{4}.spm.tools.beamforming.features.foi = [0 0];
matlabbatch{4}.spm.tools.beamforming.features.plugin.regcov.lambda = 0;
matlabbatch{5}.spm.tools.beamforming.inverse.BF(1) = cfg_dep('Covariance features: BF.mat file', substruct('.','val', '{}',{4}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{5}.spm.tools.beamforming.inverse.plugin.lcmv = struct([]);
matlabbatch{6}.spm.tools.beamforming.output.BF(1) = cfg_dep('Inverse solution: BF.mat file', substruct('.','val', '{}',{5}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_power.whatconditions.all = 1;
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_power.woi = [-Inf Inf];
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_power.contrast = 1;
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_power.result = 'singleimage';
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_power.modality = 'MEG';
matlabbatch{7}.spm.tools.beamforming.write.BF(1) = cfg_dep('Output: BF.mat file', substruct('.','val', '{}',{6}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{7}.spm.tools.beamforming.write.plugin.nifti.normalise = 'separate';
matlabbatch{7}.spm.tools.beamforming.write.plugin.nifti.space = 'mni';

author: spm


text:
% Codes for CVPR-15 work `Face Alignment by Coarse-to-Fine Shape Searching'
% Any question please contact Shizhan Zhu: zhshzhutah2@gmail.com
% Released on July 25, 2015

function varargout = indexingData(varargin)
%varargout = indexingData(ind,varargin)
%   every element in varargin will be extracted according to ind(index),
%   and the relative output will be sent to varargout.
%   PLEASE NOTE indexing only refer to the 1st dimension!

varargout = cell(1,nargin-1);
ind = varargin{1};
for i = 2:nargin
    k = varargin{i};
    if size(k,1) == 1
        k = k';
    end
    varargout{i-1} = k(ind,:);
end

end

author: zhusz


text:
classdef PrintMicroStructuresForGaussianMean < handle

    properties (Access = private)
       rho
       xi
       phi
       stressProblem
    end

    properties (Access = private)
       rhoV
       xiV
       phiV
       hMesh
       pNorm
       fileName
       outPutPath
    end

    methods (Access = public)

        function obj = PrintMicroStructuresForGaussianMean()
            obj.init();
            for iTest = 3:3%length(obj.rhoV)
                obj.rho = obj.rhoV(iTest);
                obj.xi = obj.xiV(iTest);
                obj.computePhiV();
                for iphi = 1:length(obj.phiV)
                    obj.phi = obj.phiV(iphi);
                    obj.printOptimalMicroStructure(iTest,iphi);
                end
            end
        end

    end

    methods (Access = private)

        function init(obj)
            obj.rhoV = [0.9 0.9 0.5 0.5];
            obj.xiV  = [83.7951 58.0865 39.0219 27.0665]*pi/180;
            obj.pNorm = 16;
            obj.fileName = 'OptimalMicroForStress';
            obj.outPutPath = '/home/alex/git-repos/MicroStructurePaper/';
        end

        function computePhiV(obj)
            obj.phiV = [0,pi/4,pi/2,3*pi/4,pi,obj.xi,pi - obj.xi];
        end

        function printOptimalMicroStructure(obj,iTest,iPhi)
            obj.createStressProblem(iTest,iPhi);
            obj.stressProblem.computeOptimalExponent();
            obj.stressProblem.printOptimalMicroStructure();
        end

        function createStressProblem(obj,iTest,iPhi)
            s.rho = obj.rho;
            s.txi = obj.xi;
            s.fileName = [obj.fileName,'Case',num2str(iTest),'Phi',num2str(iPhi)];
            s.phi = obj.phi;
            s.hMesh = [];
            s.pNorm = obj.pNorm;
            sProblem = OneOptimalExponentComputerAndFunctionVariation(s);
            obj.stressProblem = sProblem;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function Xk = nufft_table_adj(st, X)
%function Xk = nufft_table_adj(st, X)
% adjoint of table-based nufft interpolation.
% in
%	st		structure from nufft_init
%	X [M,?]		DTFT values
% out
%	Xk [*Kd,?]	DFT coefficients
% Copyright 2004-3-30, Jeff Fessler and Yingying Zhang, University of Michigan

dd = length(st.Kd);

% t = omega / gamma
tm = zeros(size(st.om));
for id=1:dd
	gam = 2*pi / st.Kd(id);
	tm(:,id) = st.om(:,id) / gam;
end

if size(X,1) ~= st.M
	error 'X size problem'
end

% adjoint of phase shift
if isvar('st.phase_shift') & ~isempty(st.phase_shift)
	X = X .* conj(st.phase_shift);
end

if isreal(X)
	X = complexify(X);
end

if dd == 1
	Xk = interp1_table_adj_mex(X, st.h{1}, ...
		int32(st.Jd), int32(st.Ld), tm, int32(st.Kd(1)));

elseif dd==2
	Xk = interp2_table_adj_mex(X, st.h{1}, st.h{2}, ...
		int32(st.Jd), int32(st.Ld), tm, int32(st.Kd));

elseif dd==3
	Xk = interp3_table_adj_mex(X, st.h{1}, st.h{2}, st.h{3}, ...
		int32(st.Jd), int32(st.Ld), tm, int32(st.Kd));

else
	error '> 3d not done'
end

author: thomaskuestner


text:
function tapas_ehgf_plotTraj(r)
% Plots the estimated or generated trajectories for the HGF perceptual model
% Usage example:  est = tapas_fitModel(responses, inputs); tapas_hgf_plotTraj(est);
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2020 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/>.

% Optional plotting of standard deviations (true or false)
plotsd = true;

% Set up display
scrsz = get(0,'screenSize');
outerpos = [0.2*scrsz(3),0.2*scrsz(4),0.8*scrsz(3),0.8*scrsz(4)];
figure(...
    'OuterPosition', outerpos,...
    'Name', 'HGF trajectories');

% Time axis
if size(r.u,2) > 1 && ~isempty(find(strcmp(fieldnames(r.c_prc),'irregular_intervals'))) && r.c_prc.irregular_intervals
    t = r.u(:,end)';
else
    t = ones(1,size(r.u,1));
end

ts = cumsum(t);
ts = [0, ts];

% Number of levels
l = length(r.p_prc.p)/5;

% Upper levels
for j = 1:l-1

    % Subplots
    subplot(l,1,j);

    if plotsd == true
        upperprior = r.p_prc.mu_0(l-j+1) +sqrt(r.p_prc.sa_0(l-j+1));
        lowerprior = r.p_prc.mu_0(l-j+1) -sqrt(r.p_prc.sa_0(l-j+1));
        upper = [upperprior; r.traj.mu(:,l-j+1)+sqrt(r.traj.sa(:,l-j+1))];
        lower = [lowerprior; r.traj.mu(:,l-j+1)-sqrt(r.traj.sa(:,l-j+1))];

        plot(0, upperprior, 'ob', 'LineWidth', 1);
        hold all;
        plot(0, lowerprior, 'ob', 'LineWidth', 1);
        fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
             'b', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    end
    plot(ts, [r.p_prc.mu_0(l-j+1); r.traj.mu(:,l-j+1)], 'b', 'LineWidth', 2);
    hold all;
    plot(0, r.p_prc.mu_0(l-j+1), 'ob', 'LineWidth', 2); % prior
    xlim([0 ts(end)]);
    title(['Posterior expectation of x_' num2str(l-j+1)], 'FontWeight', 'bold');
    ylabel(['\mu_', num2str(l-j+1)]);
end

% Input level
subplot(l,1,l);

if plotsd == true
    upperprior = r.p_prc.mu_0(1) +sqrt(r.p_prc.sa_0(1));
    lowerprior = r.p_prc.mu_0(1) -sqrt(r.p_prc.sa_0(1));
    upper = [upperprior; r.traj.mu(:,1)+sqrt(r.traj.sa(:,1))];
    lower = [lowerprior; r.traj.mu(:,1)-sqrt(r.traj.sa(:,1))];

    plot(0, upperprior, 'or', 'LineWidth', 1);
    hold all;
    plot(0, lowerprior, 'or', 'LineWidth', 1);
    fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
         'r', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
end
plot(ts, [r.p_prc.mu_0(1); r.traj.mu(:,1)], 'r', 'LineWidth', 2);
hold all;
plot(0, r.p_prc.mu_0(1), 'or', 'LineWidth', 2); % prior
plot(ts(2:end), r.u(:,1), '.', 'Color', [0 0.6 0]); % inputs
if ~isempty(find(strcmp(fieldnames(r),'y'))) && ~isempty(r.y)
    plot(ts(2:end), r.y(:,1), '.', 'Color', [1 0.7 0]); % responses
    title(['Response y (orange), input u (green), and posterior expectation of x_1 ', ...
           '(red) for \rho=', num2str(r.p_prc.rho), ', \kappa=', ...
           num2str(r.p_prc.ka), ', \omega=', num2str(r.p_prc.om), ...
           ', \pi_u=', num2str(r.p_prc.pi_u)], ...
          'FontWeight', 'bold');
    ylabel('y, u, \mu_1');
else
    title(['Input u (green) and posterior expectation of x_1 ', ...
           '(red) for \rho=', num2str(r.p_prc.rho), ', \kappa=', ...
           num2str(r.p_prc.ka), ', \omega=', num2str(r.p_prc.om), ...
           ', \pi_u=', num2str(r.p_prc.pi_u)], ...
          'FontWeight', 'bold');
    ylabel('u, \mu_1');
end
xlim([0 ts(end)]);
xlabel({'Trial number', ' '}); % A hack to get the relative subplot sizes right
hold off;

author: translationalneuromodeling

repo: tapas


text:
function c = ultra2ultra(c, lam_in, lam_out)
%ULTRA2ULTRA   Convert between Ultraspherical (US) expansions.
%   C_OUT = ULTRA2ULTRA(C_IN, LAM_IN, LAM_OUT) converts the vector C_IN of
%   US-LAM_IN coefficients to a vector of US-LAM_OUT coefficients.
%
%   This code is essentially a wrapper for the JAC2JAC code based on [1].
%
%   References:
%     [1] A. Townsend, M. Webb, and S. Olver, "Fast polynomial transforms
%     based on Toeplitz and Hankel matrices", submitted, 2016.
%
% See also JAC2JAC.

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

n = length(c) - 1;

% Scale input from US to Jacobi basis:
c = c./scl(lam_in, n);
% Call JAC2JAC:
c = jac2jac(c, lam_in-.5, lam_in-.5, lam_out-.5, lam_out-.5);
% Scale output from Jacobi to US:
c = c.*scl(lam_out, n);

end

function s = scl(lam, n)
% Scaling from Jacobi to US polynomials. See DLMF Table 18.3.1.
    if ( lam == 0 )
        nn = (0:n-1).';
        s = [1 ; cumprod((nn+.5)./(nn+1))];
    else
        nn = (0:n).';
        s = ( gamma(2*lam) ./ gamma(lam+.5) ) * ...
            exp( gammaln(lam+.5+nn) - gammaln(2*lam+nn) );
    end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% test_all_ct.m

list = {
'ct_beta_study'
'ct_sys'
%'de_ftab_build'
'de_ftab_curv test'
'de_ftab_fit test'
%'de_ftab_fit_sprad'
'de_ftab_fm test'
%'de_ftab_invert'
%'de_ftab_inv1 test' % 2018-10-01 fails R2018b due to lsqlin solver issue (todo)
%'de_ftab_iwater'
'de_ftab_s_iter test' % slow
%'de_ftab_xform'
'de_ftab test'
%'de_pl_denom'
%'de_pl_obj'
%'de_pl_osps'
'de_poly_eval test'
%'element_density'
%'wls_simplex'
%'xct_poly1_dercurv'
'xray_atten_interp test'
'xray_filters test'
'xray_material_file_name test'
'xray_read_atten test'
'xray_read_dens test'
'xray_read_spectra test'
}

run_mfile_local(list)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%--------------------------------------------------------------------------------------------------------
% The system is created based on the principles described in the following paper
% Jimmy SJ. Ren and Li Xu, "On Vectorization of Deep Convolutional Neural Networks for Vision Tasks",
% The 29th AAAI Conference on Artificial Intelligence (AAAI-15). Austin, Texas, USA, January 25-30, 2015
% email: jimmy.sj.ren@gmail.com
%--------------------------------------------------------------------------------------------------------
addpath applications/image_denoise/
addpath applications/image_denoise/utility/
addpath utils/
addpath cuda/
addpath mem/
addpath layers/
addpath layers_adapters/
addpath pipeline/

clearvars -global config;
clearvars -global mem;
clear;
I = im2double(imread('sample1.png'));
global config;
prepare_net(size(I, 1), size(I, 2), 'w.mat');
final_output = apply_net(config.NEW_MEM(I));
final_output = gather(final_output);
figure, imshow([I final_output]);
drawnow();

author: jimmy-ren


text:
function [y,dy] = div_gmc(dx,x,Mu_x,Mu_dx,Sigma_dxdx,Sigma_xx,Sigma_dxx,Sigma_xdx)
%GMC_DIV Derivative of the conditional Gaussian Mixture model P(dx|x)
%
%    sqrt((2*pi)^nbVar * (abs(det(Sigma))+realmin))
%   check the derivative of  log exp(- 0.5 * (dx - \Mu_{dx|x})'(\Sigma_{dx|x}^-1) (dx - \Mu_{dx|x})  )
%
%   input -----------------------------------------------------------------
%
%       o dx : (P x 1)
%
%       o x  : (Q x 1)
%
%       o Mu_dx : (P x 1), $\mu_{\dot{x}}$
%
%       o Mu_x  : (Q x 1), $\mu_x$
%
%           P + Q = D
%
%       o Sigma_dxdx   : (P x P), covariance between dx and dx
%
%       o Sigma_{dx,x} : (P x Q), covarariance between dx and x
%
%       o Sigma_{x,x}  : (Q x Q), covariance between x and x
%
%       o Sigma_{x,dx} : (Q x P),  covariance between x and dx
%
%

x   = x(:);
dx  = dx(:);
P   = size(Mu_dx,1);

% (P x 1) = (P x 1) + (P x Q)    (Q x Q)    *  (Q x 1) - (Q x 1)
% (P x 1) =           (P x 1)
invSig_xx   =  inv(Sigma_xx);
% (P x Q) = (P x Q) * (Q x Q)
A           = Sigma_dxx * invSig_xx;

% (P x 1)
Mu_dx_x = Mu_dx + A * (x - Mu_x);

% (P x P)
Sigma_dx_x = Sigma_dxdx - A * Sigma_xdx;

%  (P x N)
invSigma_dx_x = inv(Sigma_dx_x);
%norm_fac      = 0.9;%1/( (2*pi^(P/2)) * sqrt(det(Sigma_dx_x)));
%norm_fac      = -P/2 * log(pi) - 0.5*log(det(Sigma_dx_x));
denom         = sqrt((2*pi)^P * (abs(det(Sigma_dx_x))+realmin));
%norm_fac      = log(1/denom);
norm_fac      = -log(denom);
%norm_fac      = -log(sqrt((2*pi)^P) - log(abs(det(Sigma_dx_x))+realmin));

%y             =  norm_fac .* exp(-0.5 .* (dx - Mu_dx_x)' * (invSigma_dx_x *  (dx - Mu_dx_x)));
y             =   norm_fac + (-0.5 .* (dx - Mu_dx_x)' * (invSigma_dx_x *  (dx - Mu_dx_x)));

%
% (dx - Mu_dx_x)' * (invSigma_dx_x *  (dx - Mu_dx_x))
%
% (dx' * invSigma_dx_x - Mu_dx_x' * invSigma_dx_x) * (dx - Mu_dx_x)
% dx' * invSigma_dx_x * dx - Mu_dx_x' * invSigma_dx_x * dx - dx' * invSigma_dx_x * Mu_dx_x +  Mu_dx_x' * invSigma_dx_x * Mu_dx_x
%
% dx' * invSigma_dx_x * dx - 2 * Mu_dx_x' * invSigma_dx_x * dx  +  Mu_dx_x' * invSigma_dx_x * Mu_dx_x
%
% dx' * invSigma_dx_x * dx - 2 * (Mu_dx + A * (x - Mu_x))' * invSigma_dx_x * dx  +  (Mu_dx + A * (x - Mu_x))' * invSigma_dx_x * (Mu_dx + A * (x - Mu_x))
% 'here'
% tmp1  = (x - Mu_x)' * A' * invSigma_dx_x;
%
% part1 = dx' * invSigma_dx_x * dx  - 2 *Mu_dx' *invSigma_dx_x * dx  - 2 * tmp1 * dx;
%
% part1 + (Mu_dx' * invSigma_dx_x + tmp1) * (Mu_dx + A * (x - Mu_x))
%  'max expansion'
% part1 + Mu_dx' * invSigma_dx_x * Mu_dx + tmp1 * Mu_dx + Mu_dx' * invSigma_dx_x * (A * (x - Mu_x)) + tmp1 * (A * (x - Mu_x))
%
% part1 + Mu_dx' * invSigma_dx_x * Mu_dx + tmp1 * Mu_dx + tmp1 * Mu_dx + tmp1 * (A * (x - Mu_x))
%
%
% dx' * invSigma_dx_x * dx  - 2 *Mu_dx' *invSigma_dx_x * dx +  Mu_dx' * invSigma_dx_x * Mu_dx  - 2 * tmp1 * dx  + tmp1 * Mu_dx + tmp1 * Mu_dx + tmp1 * (A * (x - Mu_x))
%
% dx' * invSigma_dx_x * dx  - 2 *Mu_dx' *invSigma_dx_x * dx +  Mu_dx' * invSigma_dx_x * Mu_dx +  tmp1 * (- 2 * dx  + Mu_dx + Mu_dx + (A * (x - Mu_x)));

%fac = (x - Mu_x)' * A' * invSigma_dx_x;
% y2 = dx' * invSigma_dx_x * dx   -2 *Mu_dx' *invSigma_dx_x * dx +  Mu_dx' * invSigma_dx_x * Mu_dx + fac * (-2*dx + 2*Mu_dx + A * (x - Mu_x));
%
% y2 = exp(-0.5 .* y2);

if nargout > 1
    %
    %          -            (1 x P)          (1 x P)                      (1 x P)
%    dy = y .* ((dx' * invSigma_dx_x)' - (Mu_dx' * invSigma_dx_x)' - fac' );
 %   dy = ((dx' * invSigma_dx_x)' - (Mu_dx' * invSigma_dx_x)' - ( (x - Mu_x)' * A' * invSigma_dx_x)' );
    dy = invSigma_dx_x' * ( dx - Mu_dx  -  ((x - Mu_x)' * A')' );

end

end

author: epfl-lasa


text:
function [ y, j, ierror ] = yj_check_hebrew ( y, j )

%*****************************************************************************80
%
%% YJ_CHECK_HEBREW checks a Hebrew YJ date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, J, the YJ date.
%
%    Output, integer IERROR, is 0 if no serious error was found in
%    the date, and 1 otherwise.
%

%
%  Check the year.
%
  [ y, ierror ] = y_check_hebrew ( y );

  if ( ierror ~= 0 )
    return
  end
%
%  Make sure J is not too small or too big.
%
  [ y, j ] = j_borrow_hebrew ( y, j );

  [ y, j ] = j_carry_hebrew ( y, j );

  return
end


text:
% this function solves the odometry calibration problem
% given a measurement matrix Z.
% We assume that the information matrix is the identity
% for each of the measurements
% Every row of the matrix contains
% z_i = [u'x, u'y, u'theta, ux, uy, ytheta]
% Z:	The measurement matrix
% X:	the calibration matrix
% returns the correction matrix X
function X = ls_calibrate_odometry(Z)
  % initial solution (the identity transformation)
  X = eye(3);

  % TODO: initialize H and b of the linear system
  H = zeros(9);
  b = zeros(9,1);
  omega = eye(3);

  % TODO: loop through the measurements and update H and b
  % You may call the functions error_function and jacobian, see below
  % We assume that the information matrix is the identity.
  for i=1:size(Z,1)
	H = H + jacobian(i,Z)'*omega*jacobian(i,Z);
	b = b + jacobian(i,Z)'*omega'*error_function(i,X,Z);
  endfor
  % TODO: solve and update the solution
  delX = -H\b;
  X = X + [delX(1:3)';delX(4:6)';delX(7:9)'];
end

% this function computes the error of the i^th measurement in Z
% given the calibration parameters
% i:	the number of the measurement
% X:	the actual calibration parameters
% Z:	the measurement matrix, each row contains first the scan-match result
%       and then the motion reported by odometry
% e:	the error of the ith measurement
function e = error_function(i, X, Z)
  % TODO compute the error of each measurement
  e = Z(i,1:3)' - X*Z(i,4:6)';
end

% derivative of the error function for the ith measurement in Z
% i:	the measurement number
% Z:	the measurement matrix
% J:	the jacobian of the ith measurement
function J = jacobian(i, Z)
  % TODO compute the Jacobian
  J = -[Z(i,4) Z(i,5) Z(i,6) zeros(1,6); zeros(1,3) Z(i,4) Z(i,5) Z(i,6) zeros(1,3); zeros(1,6) Z(i,4) Z(i,5) Z(i,6)];
end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
function [err] = writeSPMmat (matin, matout)
%
%   [err] = writeSPMmat (matin, [matout]);
%
%Purpose:
%   Write an ascii version of SPM's .mat files
%   containing a 4x4 transformation matrix.
%   The output of this function can be used
%   with 3dWarp to create AFNI-formatted data from
%   SPM datasets.
%
%Input Parameters:
%   matin : Name of .mat file containing transform
%   matout: (optional) Name of output file.
%           If matout is not specified, the _ASCII
%           extension is added to matin
%
%Output Parameters:
%   err : 0 No Problem
%       : 1 Problems
%
%
%More Info :
%   3dWarp -help
%   https://afni.nimh.nih.gov/ssc/ziad/SUMA/SUMA_doc.htm
%   (search for SPM)
%
%     Author : Ziad Saad
%     Date : Thu Aug 7 11:31:54 EDT 2003
%     SSCC/NIMH/ National Institutes of Health, Bethesda Maryland

%Define the function name for easy referencing
FuncName = 'writeSPMmat';

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;

%does the input file exist ?
matin = deblank(matin);
matin = deblank(fliplr(matin));
matin = fliplr(matin);

if (exist(matin) ~= 2),
   fprintf(2,'Error %s:\nFile %s not found.\n', FuncName, matin);
   return;
end

%get its prefix
if (strmatch('tam.', fliplr(matin)) == 1),
   prefixin = matin(1:length(matin)-4);
else
   prefixin = matin;
end

%load it
S = load (matin);
if (~isfield(S,'M')),
   fprintf(2,'Error %s:\nFailed to find M in %s\n', FuncName, matin);
   return;
end

if (size(S.M,2) ~= 4),
   fprintf(2,'Error %s:\nM does not have 4 columns!\nM''s size is %g x %g',...
      FuncName, size(S.M,1), size(S.M,2));
   return;
end

if (size(S.M,1) ~= 3 & size(S.M,1) ~= 4),
   fprintf(2,'Error %s:\nM does not have 3 or 4 rows!\nM''s size is %g x %g',...
      FuncName, size(S.M,1), size(S.M,2));
   return;
end

if (size(S.M,1) == 4),
   if (S.M(4,1) ~= 0 | S.M(4,2) ~= 0 | S.M(4,3) ~= 0 | S.M(4,4) ~= 1),
      beep;
      fprintf(2,'Warning %s:\nAFNI expects the 4th row of M to be [0 0 0 1].\n', FuncName);
      fprintf(2,'This row is currently [%g %g %g %g] and will be ignored by AFNI.\n', ...
         S.M(4,1), S.M(4,2), S.M(4,3), S.M(4,4));
      fprintf(2,'Proceed with caution.\n');
      beep;
   end
end

%write output file
if (nargin == 1),
   matout = [prefixin, '_ASCII.mat'];
end
if (exist(matout) == 2),
   fprintf(2,'Error %s:\nOutput file %s exists already.\nWill not overwrite.\n',...
    FuncName, matout);
   return;
end

fid = fopen(matout,'w');
if (fid < 0),
   fprintf(2,'Error %s:\nFailed in opening %s for writing.\nCheck permissions and disk space.\n',...
      FuncName, matout);
   return;
end
fprintf (fid,'#Conversion of SPM''s %s \n', matin);
fprintf (fid,'# transformation matrix file to ascii format.\n#\n');
fprintf (fid,'#To convert SPM datasets to AFNI''s format, use \n');
fprintf (fid,'# 3dWarp including these options:\n');
fprintf (fid,'# -matvec_in2out %s -matvec_fsl \n', matout);
fprintf (fid,'#See 3dWarp -help for more info.\n#\n');
for (i=1:1:4),
   fprintf (fid,'\t%f\t%f\t%f\t%f\n',...
      S.M(i,1), S.M(i,2), S.M(i,3), S.M(i,4));
end

fclose(fid);
err = 0;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = process_cohere2_time_2021( varargin )
% PROCESS_COHERE2_TIME_2021: Compute the time-resolved coherence between all the pairs of signals, in one file.

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Time-resolved coherence AxB [2021]';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Connectivity';
    sProcess.Index       = 660;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Connectivity';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',     'results',  'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 2;
    sProcess.nMinFiles   = 1;
    sProcess.isPaired    = 1;

    % === CONNECT INPUT
    sProcess = process_corr2('DefineConnectOptions', sProcess);
    % === REMOVE EVOKED REPONSE
    sProcess.options.removeevoked.Comment = 'Remove evoked response from each trial';
    sProcess.options.removeevoked.Type    = 'checkbox';
    sProcess.options.removeevoked.Value   = 0;
    sProcess.options.removeevoked.Group   = 'input';
    % === Time window
    sProcess.options.slide_win.Comment = 'Sliding time window duration:';
    sProcess.options.slide_win.Type    = 'value';
    sProcess.options.slide_win.Value   = {.350, 'ms', []};
    % === WinOverlap for Sliding window (Time)
    sProcess.options.slide_overlap.Comment = 'Sliding window overlap:';
    sProcess.options.slide_overlap.Type    = 'value';
    sProcess.options.slide_overlap.Value   = {50, '%', []};
    % === COHERENCE METHOD
    sProcess.options.cohmeasure.Comment = {...
        ['Magnitude-squared Coherence
' ...
        '|C|^2 = |Gxy|^2/(Gxx*Gyy)'], ...
        ['Imaginary Coherence (2019)
' ...
        'IC    = |imag(C)|'], ...
        ['Lagged Coherence (2019)
' ...
        'LC    = |imag(C)|/sqrt(1-real(C)^2)'], ...
        [' Imaginary Coherence (before 2019)
' ...
        ' IC    = imag(C)^2 / (1-real(C)^2) ']; ...
        'mscohere', 'icohere2019','lcohere2019', 'icohere'};
    sProcess.options.cohmeasure.Type    = 'radio_label';
    sProcess.options.cohmeasure.Value   = 'mscohere';
    % === WINDOW LENGTH
    sProcess.options.win_length.Comment = 'Window length for PSD estimation:';
    sProcess.options.win_length.Type    = 'value';
    sProcess.options.win_length.Value   = {1, 's', []};
    % === OVERLAP
    sProcess.options.overlap.Comment = 'WinOverlap for PSD estimation:' ;
    sProcess.options.overlap.Type    = 'value';
    sProcess.options.overlap.Value   = {50, '%', []};
    % === HIGHEST FREQUENCY OF INTEREST
    sProcess.options.maxfreq.Comment = 'Highest frequency of interest:';
    sProcess.options.maxfreq.Type    = 'value';
    sProcess.options.maxfreq.Value   = {60,'Hz',2};
    % === OUTPUT FILE TAG
    sProcess.options.commenttag.Comment = 'File tag: ';
    sProcess.options.commenttag.Type    = 'text';
    sProcess.options.commenttag.Value   = '';
    sProcess.options.commenttag.Group   = 'output';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputA, sInputB) %#ok<DEFNU>
    % Initialize returned values
    OutputFiles = {};
    % Output mode 2021: Forcing the average cross-spectra of input files (one output file)
    sProcess.options.outputmode.Value = 'avgcoh';
    % Input options
    OPTIONS = process_corr2('GetConnectOptions', sProcess, sInputA, sInputB);
    if isempty(OPTIONS)
        return
    end
    CommentTag = sProcess.options.commenttag.Value;
    % Metric options
    OPTIONS.Method = 'cohere';
    OPTIONS.RemoveEvoked  = sProcess.options.removeevoked.Value;
    OPTIONS.WinLen        = sProcess.options.win_length.Value{1};
    OPTIONS.MaxFreq       = sProcess.options.maxfreq.Value{1};
    OPTIONS.CohOverlap    = 0.50;
    OPTIONS.pThresh       = 0.05;
    OPTIONS.isSave        = 0;
    OPTIONS.CohMeasure    = sProcess.options.cohmeasure.Value;
    % Sliding time windows options
    WinLength  = sProcess.options.slide_win.Value{1};
    WinOverlap = sProcess.options.slide_overlap.Value{1};

    % Read time information
    TimeVectorA  = in_bst(sInputA(1).FileName, 'Time');
    sfreq        = round(1/(TimeVectorA(2) - TimeVectorA(1)));
    % Get time window of first fileA if none specified in parameters
    if isempty(OPTIONS.TimeWindow)
        OPTIONS.TimeWindow = TimeVectorA([1, end]);
    end
    % Select input time window
    TimeVectorA = TimeVectorA((TimeVectorA >= OPTIONS.TimeWindow(1)) & (TimeVectorA <= OPTIONS.TimeWindow(2)));
    nTime       = length(TimeVectorA);

    % Compute sliding windows length
    Lwin  = round(WinLength * sfreq);
    Loverlap = round(Lwin * WinOverlap / 100);
    Nwin = floor((nTime - Loverlap) ./ (Lwin - Loverlap));
    % If window is bigger than the data
    if (Lwin > nTime)
        bst_report('Error', sProcess, sInputA, 'Sliding window for the coherence estimation is too long compared with the epochs in input.');
        return;
    end

    % Check that time is the same for FilesB
    if ~isempty(sInputB)
        TimeVectorB = in_bst(sInputB(1).FileName, 'Time');
        if (length(TimeVectorA) ~= length(TimeVectorB))
            bst_report('Error', sProcess, sInputA, 'Files A and B must share the same time vector.');
            return;
        end
    end

    % Get progress bar position
    posProgress = bst_progress('get');
    % Loop over all the time windows
    for iWin = 1:Nwin
        % Set the progress bar at the same level at every iteration
        bst_progress('set', posProgress);
        % Select time window
        iTimes = (1:Lwin) + (iWin-1)*(Lwin - Loverlap);
        OPTIONS.TimeWindow = TimeVectorA(iTimes([1,end]));
        % Compute metric
        if ~isempty(sInputB)
            ConnectMat = bst_connectivity({sInputA.FileName}, {sInputB.FileName}, OPTIONS);
        else
            ConnectMat = bst_connectivity({sInputA.FileName}, [], OPTIONS);
        end
        % Processing errors
        if isempty(ConnectMat) || ~iscell(ConnectMat) || ~isstruct(ConnectMat{1}) || isempty(ConnectMat{1}.TF)
            bst_report('Error', sProcess, sInputA, 'Coherence for the selected time segment could not be calculated.');
            return;
        end
        % Start a new brainstorm structure
        if (iWin == 1)
            NewMat = ConnectMat{1};
            NewMat.Time = OPTIONS.TimeWindow(1);
            NewMat.TimeBands = [];
        % Add next time point
        else
            NewMat.TF(:,iWin,:) = ConnectMat{1}.TF;
            NewMat.Time(end+1) = OPTIONS.TimeWindow(1);
        end
    end

    % Fix time vector
    if (length(NewMat.Time) == 1)
        bst_report('Warning', sProcess, sInputA, 'Only one sliding time window could be estimated.');
        NewMat.Time = [TimeVectorA(1), TimeVectorA(end)];
    end
    % Add comment tag
    if ~isempty(CommentTag)
        NewMat.Comment = [NewMat.Comment ' | ' CommentTag];
    end
    % File tag
    if (length(NewMat.RefRowNames) == 1)
        fileTag = 'connect1';
    else
        fileTag = 'connectn';
    end
    % Output filename
    sOutputStudy = bst_get('Study', OPTIONS.iOutputStudy);
    OutputFiles{1} = bst_process('GetNewFilename', bst_fileparts(sOutputStudy.FileName), ['timefreq_' fileTag '_cohere_time']);
    % Save file
    bst_save(OutputFiles{1}, NewMat, 'v6');
    % Add file to database structure
    db_add_data(OPTIONS.iOutputStudy, OutputFiles{1}, NewMat);
end

author: brainstorm-tools


text:
function phiFaceAverage = upwindMean(phi, u)
% This function gets the value of the field variable phi defined
% over the MeshStructure and calculates the upwind average on
% the cell faces, for a uniform mesh based on the direction of the velocity
% vector u
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Written by Ali A. Eftekhari
% See the license file

% extract data from the mesh structure

d = phi.domain.dimension;
if (d ==1) || (d==1.5) || (d==1.8)
	phiFaceAverage = upwindMean1D(phi, u);
elseif (d == 2) || (d == 2.5) || (d==2.8)
	phiFaceAverage = upwindMean2D(phi, u);
elseif (d == 3) || (d==3.2)
    phiFaceAverage = upwindMean3D(phi, u);
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
% Read single/multi channel waveforms
% Inputs:
%   files: a cell array of files to be read. Each cell represents one
%     sentence could be a cell and could be array of files itself if it is
%     multi-channel recording. A filename may also contains the starting
%     and ending time to be read, e.g. "abc.wav 0.2 0.5" will read the
%     segment from 0.2s to 0.5s of abc.wav.
%   reader: a structure that may contains following fields
%       fs: sampling frequency, if known in advance.
%       array: set to 1 if the sentence has multiple channels
%       multiArrayFiles: set to 1 if each channel is stored in one file.
%       useChannel: an array of integers specifying which channel(s) to be
%           returned.
%
%   Author: Xiong Xiao, Nanyang Technological University, Singapore
%   Last modified: 27 Jul 2016
%
function [all_wav, fs] = Reader_waveform(files, reader)

for i=1:length(files)
    if isfield(reader, 'array') && reader.array
        % read in multichannel waveforms
        if isfield(reader, 'multiArrayFiles') && reader.multiArrayFiles     % when the different channels are stored in different files
            for j=1:length(files{i})
                [tmp_wav, fs] = Reader_waveform_core(files{i}{j}, reader.fs);
                if j==1
                    wav = zeros(length(files{i}), length(tmp_wav));
                end
                wav(j,:) = tmp_wav;
            end
        else                                                                % when the different channels are stored in the same file
            [wav, fs] = Reader_waveform_core(files{i}, reader.fs);
        end
        % optional selection of channels
        if isfield(reader, 'useChannel')
            wav = wav(reader.useChannel, :);
        end
    else
        % read in single channel waveform
        [wav, fs] = Reader_waveform_core(files{i}, reader);
    end
    if isfield(reader, 'precision')
        switch lower(reader.precision)
            case 'int16'
                wav = StoreWavInt16(wav);
            case 'single'
                wav = single(wav);
            case 'double'
                wav = double(wav);
        end
    end
    all_wav{i} = wav;
end
end

%%
function [wav, fs] = Reader_waveform_core(file, reader)

words = ExtractWordsFromString_v2(file);

selectChannel = 0;
selectTime = 0;

filename = words{1};
if length(words)==4 || length(words)==2
    channelID = str2num(words{end});
    selectChannel = 1;
end
if length(words)>=3
    time1 = str2num(words{2});
    time2 = str2num(words{3});
    selectTime = 1;
end

fs = ReturnFieldWithDefaultValue(reader, 'fs', -1);
big_endian = ReturnFieldWithDefaultValue(reader, 'big_endian', 0);

switch lower(reader.name)
    case 'raw'
        [wav] = read08(filename, big_endian);   % you have to set fs correctly in the reader.
    otherwise
        if selectTime
            [wav, fs] = ReadAudioSeg(filename, time1, time2, fs);
        else
            [wav, fs] = audioread(filename);
        end
end

if selectChannel
    wav = wav(:,channelID);
end

wav = wav';

end

author: singaxiong

repo: SignalGraph


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% calculation of the 'epsgg' matrix for circular holes using
%%% analytical expression; the matrix is symmetric, i.e E'=E
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function epsi=epsgg(r, na, nb, b1, b2, N1, N2)
N = N1*N2;
epsilon=zeros(N1,N2);
epsi=zeros(N,N);
%%% filling factor, ratio between the area occupied by the cylinders and area of the unit cell
f=(2*pi/sqrt(3))*r^2;

for  l=1:N1
      for m=1:N2
          for  n=1:N1
              for p=1:N2
                  GGx=(l-n)*b1(1)+(m-p)*b2(1);
                  GGy=(l-n)*b1(2)+(m-p)*b2(2);
                  GG=sqrt(GGx^2+GGy^2); x=2*pi*GG;
                  %%% GG is a scalar which changes with iteration and goes into Bessel function argument
                  if (GG~=0)
                      epsilon(p,n)=2*f*(na^2-nb^2)*besselj(1,x*r)/(x*r);
                  else
                      epsilon(p,n)=f*na^2+(1-f)*nb^2;
                  end ;
              end;
           end;
           u=(l-1)*N2+m; %%% this is the line index
           epsi(u,:)=reshape(epsilon,1,N);
        end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: Sable

repo: mcbench-benchmarks


text:
function init(fname, nbytes, opts)
% Initialise binary file on disk
% FORMAT init(fname, nbytes[, opts])
% fname   - filename
% nbytes  - data size {bytes}
% opts    - optional structure with fields:
%   .offset   - file offset {bytes} [default: 0]
%   .wipe     - overwrite exisiting values with 0 [default: false]
%   .truncate - truncate file if larger than requested size [default: true]
%
% This function is normally called by file_array/initialise
% _______________________________________________________________________
% Copyright (C) 2013 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: init.m 5458 2013-05-01 14:32:23Z guillaume $

%-This is merely the help file for the compiled routine
error('init.c not compiled - see Makefile');


text:
function besk0_test ( )

%*****************************************************************************80
%
%% BESK0_TEST tests R4_BESK0 and R8_BESK0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BESK0_TEST:\n' );
  fprintf ( 1, '  Test BESK0_VALUES, R4_BESK0, R8_BESK0\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X       BESK0(X)\n' );
  fprintf ( 1, '                  R4_BESK0(X)         Diff\n' );
  fprintf ( 1, '                  R8_BESK0(X)         Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx1 ] = bessel_k0_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_besk0 ( single ( x ) );
    fx3 = r8_besk0 ( x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.6g\n', x, fx1 );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end


text:
name = 'david50kf';
name = 'hand';
name = 'elephant-50kv';
name = 'bunny';
options.name = name;

path(path, '../toolbox_graph_data/off/');

rep = ['results/curvature/' name '/'];
if not(exist(rep))
    mkdir(rep);
end

[vertex,face] = read_mesh(name);
n = size(vertex,2);

options.symmetrize = 0;
options.normalize = 1;
type = 'combinatorial';
type = 'conformal';
if not(exist('L'))
    L = compute_mesh_laplacian(vertex,face,type,options);
end

% initial conditions
sigma = .015;
if strcmp(name, 'bunny')
    sigma = .002;
end
npoints = 10; m =0;
f = zeros(n,1);
while m<npoints
    clf;
    options.face_vertex_color = f;
    plot_mesh(vertex,face,options);
    colormap jet(256);
    disp('Click on mesh and then hit enter.');
    pause; p = select3d;
    if not(isempty(p))
        d = compute_distance_to_points(vertex,p)';
        f = f + (-1)^mod(m,2) *exp( -d/(2*sigma.^2) );
        m = m+1;
    end
end

Tmax = 80;
dt = .5;

rep = 'results/wave-equation/';
if not(exist(rep))
    mkdir(rep);
end

fprev = f;
m = 1;
for i=1:round(Tmax/dt)
    f1 = f;
    f = 2*f - fprev - dt^2 * L*f;
    fprev = f1;
    if mod(i,round(4/dt))==1
        a = f; a(a==max(a)) = max(abs(f));
        a(a==min(a)) = -max(abs(f));
        options.face_vertex_color = a;
        clf;
        plot_mesh(vertex,face,options);
        colormap jet(256);
        saveas(gcf, [rep name '-wave-eq-' num2string_fixeddigit(m,2) '.png'], 'png');
        m = m+1;
    end
end

author: gpeyre

repo: matlab-toolboxes


text:
function clearWarningErrorLabels()
%clearWarningErrorLabels Summary of this function goes here
%   Detailed explanation goes here
    global warn_error_label_handles;

    warnErrorLbls = warn_error_label_handles;
    for(i=1:length(warnErrorLbls)) %#ok<*NO4LP>
        h = warnErrorLbls(i);
        h.Visible = 'off';
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function test_suite=test_naive_bayes_classifier_searchlight
% tests for cosmo_naive_bayes_classifier_searchlight
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_naive_bayes_classifier_searchlight_tiny
    ds=cosmo_synthetic_dataset('ntargets',25,'nchunks',4,'size','small');
    nh=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);

    nh.neighbors={1:size(ds.samples,2)};
    nh.fa=cosmo_slice(nh.fa,1,2,'struct');

    opt=struct();
    opt.partitions=cosmo_nfold_partitioner(ds);
    opt.output='winner_predictions';
    opt.progress=false;

    x=cosmo_naive_bayes_classifier_searchlight(ds,nh,opt);
    assertEqual(x.fa,nh.fa);
    assertEqual(x.a,nh.a);

    nfolds=numel(opt.partitions.train_indices);
    for fold=1:nfolds
        tr=opt.partitions.train_indices{fold};
        te=opt.partitions.test_indices{fold};

        pred=cosmo_classify_naive_bayes(ds.samples(tr,:),...
                                            ds.sa.targets(tr),...
                                            ds.samples(te,:));
        y=x.samples(te,:);
        assertEqual(y,pred);

    end
    % compare with standard searchlight
    assert_same_output_as_classifical_searchlight(ds,nh,opt);

function test_naive_bayes_classifier_searchlight_basics
    ds=cosmo_synthetic_dataset('ntargets',25,'nchunks',4,'size','small');
    nh=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);

    nh.neighbors={1:size(ds.samples,2)};
    nh.fa=cosmo_slice(nh.fa,1,2,'struct');

    opt=struct();
    opt.partitions=cosmo_nfold_partitioner(ds);
    opt.output='winner_predictions';
    opt.progress=false;

    x=cosmo_naive_bayes_classifier_searchlight(ds,nh,opt);
    assertEqual(x.fa,nh.fa);
    assertEqual(x.a,nh.a);

    % compare with standard searchlight
    assert_same_output_as_classifical_searchlight(ds,nh,opt);

    opt.output='accuracy';
    xacc=cosmo_naive_bayes_classifier_searchlight(ds,nh,opt);
    assertEqual(xacc.samples,mean(bsxfun(@eq,x.samples,x.sa.targets)));
    assert_same_output_as_classifical_searchlight(ds,nh,opt);

function test_naive_bayes_classifier_searchlight_multiple_pred
    ds=cosmo_synthetic_dataset('ntargets',25,'nchunks',5,'size','small');
    nh=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);

    opt=struct();
    % multiple predictions per fold
    opt.partitions=cosmo_nchoosek_partitioner(ds,2);
    opt.output='winner_predictions';
    opt.progress=false;

    x=cosmo_naive_bayes_classifier_searchlight(ds,nh,opt);
    assertEqual(x.fa,nh.fa);
    assertEqual(x.a,nh.a);

    % compare with standard searchlight
    assert_same_output_as_classifical_searchlight(ds,nh,opt);

    opt.output='accuracy';
    assert_same_output_as_classifical_searchlight(ds,nh,opt);

function assert_same_output_as_classifical_searchlight(ds,nh,opt)
    opt.progress=false;
    x=cosmo_naive_bayes_classifier_searchlight(ds,nh,opt);

    opt.classifier=@cosmo_classify_naive_bayes;
    y=cosmo_searchlight(ds,nh,@cosmo_crossvalidation_measure,opt);

    sx=x.samples;
    sy=y.samples;
    assertElementsAlmostEqual(sx,sy);

    x=rmfield(x,'samples');
    y=rmfield(y,'samples');
    assertEqual(x,y);

function test_naive_bayes_classifier_searchlight_exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                cosmo_naive_bayes_classifier_searchlight(varargin{:}),'');

    ds=cosmo_synthetic_dataset('size','small','nchunks',4);
    nh=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);

    opt=struct();
    opt.progress=false;
    opt.partitions=cosmo_nchoosek_partitioner(ds,2);

    opt.output='foo';
    aet(ds,nh,opt);

    % missing samples, so illegal partitions
    opt.output='winner_predictions';
    ds=cosmo_slice(ds,ds.sa.chunks<=2);
    opt.partitions=cosmo_nfold_partitioner(ds);
    aet(ds,nh,opt);

    % unsupported output
    opt.output='fold_predictions';
    ds_bad=ds;
    ds_bad=cosmo_slice(ds_bad,ds_bad.sa.chunks<=2);
    opt.partitions=cosmo_nfold_partitioner(ds_bad);
    aet(ds,nh,opt);

function test_naive_bayes_classifier_searchlight_deprecations
    ds=cosmo_synthetic_dataset('size','small','nchunks',4);
    nh=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);

    %
    opt=struct();
    opt.progress=false;
    opt.partitions=cosmo_nchoosek_partitioner(ds,2);
    opt.output='predictions';

    orig_warning_state=cosmo_warning();
    cleaner=onCleanup(@()cosmo_warning(orig_warning_state));
    cosmo_warning('reset');
    cosmo_warning('off');

    % no warnings
    w=cosmo_warning();
    assertEqual(numel(w.shown_warnings),0);

    % output='predictions' is deprecated, so expect a warning
    cosmo_naive_bayes_classifier_searchlight(ds,nh,opt);

    w=cosmo_warning();
    assertEqual(numel(w.shown_warnings),1);

function test_naive_bayes_classifier_searchlight_partial_partitions
    nchunks=4;
    ds=cosmo_synthetic_dataset('ntargets',5,'nchunks',nchunks,...
                        'size','small');
    nh=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);

    partitions=cosmo_nfold_partitioner(ds);
    nsamples=size(ds.samples,1);
    prediction_count=zeros(nsamples,1);
    for k=1:nchunks
        with_missing=partitions.test_indices{k}(2:end);
        partitions.test_indices{k}=with_missing;
        prediction_count(with_missing)=prediction_count(with_missing)+1;
    end

    opt=struct();
    opt.progress=false;
    opt.partitions=partitions;
    opt.output='winner_predictions';

    res=cosmo_naive_bayes_classifier_searchlight(ds,nh,opt);

    ds_sa=rmfield(ds.sa,'chunks');
    assertEqual(res.sa,ds_sa);
    cosmo_check_dataset(res);

    assert_same_output_as_classifical_searchlight(ds,nh,opt);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
% rbm_get_hidden
% Copyright (C) 2011 KyungHyun Cho, Tapani Raiko, Alexander Ilin
%
%This program is free software; you can redistribute it and/or
%modify it under the terms of the GNU General Public License
%as published by the Free Software Foundation; either version 2
%of the License, or (at your option) any later version.
%
%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License
%along with this program; if not, write to the Free Software
%Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
%
function [h] = rbm_get_hidden(x0, R)

n_visible = size(R.W, 1);
n_hidden = size(R.W, 2);

if R.data.binary == 1
    h = sigmoid(bsxfun(@plus,x0 * R.W, R.hbias'));
else
    h = sigmoid(bsxfun(@plus, bsxfun(@rdivide, x0, R.sigmas.^2) * R.W, R.hbias'));
end

author: kyunghyuncho

repo: deepmat


text:
function []=panel3plot(endo,Units,Xi,Yi,N,n,m,p,k,T,beta_median,beta_gibbs,It,Bu,decimaldates1,stringdates1,pref,cband,favar)

% actual vs. fitted

% initiate predictions
Ytilde=[];
% loop over units
for ii=1:N
% obtain a point estimate betatilde of the VAR coefficients
% use the point estimate, which is simply the median
betatilde=beta_median(:,:,ii);
Btilde=reshape(betatilde,k,n);
% compute fitted values for the unit
Ytilde(:,:,ii)=Xi(:,:,ii)*Btilde;
end
% then plot the figure
if pref.plot
actualfitted=figure('Tag','BEARresults');
set(actualfitted,'Color',[0.9 0.9 0.9]);
set(actualfitted,'name','actual vs fitted')
% initiate the count
count=0;
% loop over units
for ii=1:N
   % loop over endogenous variables
   for jj=1:n
   % increment count
   count=count+1;
   % then plot
   subplot(N,n,count)
   hold on
   plot(decimaldates1,Yi(:,jj,ii),'Color',[0 0 0],'LineWidth',2);
   plot(decimaldates1,Ytilde(:,jj,ii),'Color',[1 0 0],'LineWidth',2);
   hold off
   set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
      % top labels
      if count<=n
      title(endo{count,1},'FontWeight','normal');
      end
      % side labels
      if jj==1
      ylabel(Units{ii,1},'FontWeight','normal');
      end
      % legend
      if ii==1 && jj==1
      plotlegend=legend('actual','fitted');
      set(plotlegend,'FontName','Times New Roman');
      end
   end
end
end
% save on Excel
% create the cell that will be saved on excel
afcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},T+3,1);
horzspace=repmat({''},3,4*n);
% loop over units
for ii=1:N
% initiate the cell of results
unitcell={};
   % loop over endogenous variables (horizontal dimension)
   for jj=1:n
   % create a header
   header=[{[Units{ii,1} ': ' endo{jj,1}]} {''} {''};{''} {''} {''};{''} {'sample'} {'fitted'}];
   % complete the cell
   endocell=[[header;stringdates1 num2cell(Yi(:,jj,ii)) num2cell(Ytilde(:,jj,ii))] vertspace];
   % concatenate to the previous parts of unitcell
   unitcell=[unitcell endocell];
   end
% concatenate to the previous parts of afcell
afcell=[afcell;horzspace;unitcell];
end
% trim
afcell=afcell(4:end,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),afcell,'actual fitted','B2');
end

% residuals

% inititate the matrix of residuals
EPStilde=[];
% loop over units
for ii=1:N
% estimate the residuals for this unit
EPStilde(:,:,ii)=Yi(:,:,ii)-Ytilde(:,:,ii);
end
if pref.plot
% using these values, plot
residuals=figure('Tag','BEARresults');
set(residuals,'Color',[0.9 0.9 0.9]);
set(residuals,'name','residuals')
% initiate the count
count=0;
% loop over units
for ii=1:N
   % loop over endogenous variables
   for jj=1:n
   % increment count
   count=count+1;
   % then plot
   subplot(N,n,count)
   plot(decimaldates1,EPStilde(:,jj,ii),'Color',[0 0 0],'LineWidth',2);
   set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
      % top labels
      if count<=n
      title(endo{count,1},'FontWeight','normal');
      end
      % side labels
      if jj==1
      ylabel(Units{ii,1},'FontWeight','normal');
      end
   end
end
end
% save on Excel
% create the cell that will be saved on excel
residcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},T+3,1);
horzspace=repmat({''},3,3*n);
% loop over units
for ii=1:N
% initiate the cell of results
unitcell={};
   % loop over endogenous variables (horizontal dimension)
   for jj=1:n
   % create a header
   header=[{[Units{ii,1} ': ' endo{jj,1}]} {''};{''} {''};{''} {'median'}];
   % complete the cell
   endocell=[[header;stringdates1 num2cell(EPStilde(:,jj,ii))] vertspace];
   % concatenate to the previous parts of unitcell
   unitcell=[unitcell endocell];
   end
% concatenate to the previous parts of afcell
residcell=[residcell;horzspace;unitcell];
end
% trim
residcell=residcell(4:end,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),residcell,'resids','B2');
end

% steady-state

% initiate the cell storing the steady-state draws
ss_record={};
ss_estimates={};
% run the Gibbs sampler for the steady-state (common to all units)
ss_record=bear.ssgibbs(n,m,p,k,Xi(:,:,1),beta_gibbs,It,Bu,favar);
% obtain point estimates and credibility interval
ss_estimates=bear.ssestimates(ss_record,n,T,cband);
if pref.plot
% then plot the figure
sstate=figure('Tag','BEARresults');
set(sstate,'Color',[0.9 0.9 0.9]);
set(sstate,'name','steady-state')
% initiate the count
count=0;
% loop over units
for ii=1:N
   % loop over endogenous variables
   for jj=1:n
   % increment count
   count=count+1;
   % then plot
   subplot(N,n,count)
   hold on
   Xpatch=[decimaldates1' fliplr(decimaldates1')];
   Ypatch=[ss_estimates{jj,1}(1,:) fliplr(ss_estimates{jj,1}(3,:))];
   Fpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
   set(Fpatch,'facealpha',0.6);
   set(Fpatch,'edgecolor','none');
   ss=plot(decimaldates1,ss_estimates{jj,1}(2,:),'Color',[0.4 0.4 1],'LineWidth',2);
   actual=plot(decimaldates1,Yi(:,jj,ii),'Color',[0 0 0],'LineWidth',2);
   plot([decimaldates1(1,1),decimaldates1(end,1)],[0 0],'k--');
   hold off
   minband=min(min(ss_estimates{jj,1}(1,:)),min(Yi(:,jj,ii)));
   maxband=max(max(ss_estimates{jj,1}(3,:)),max(Yi(:,jj,ii)));
   space=maxband-minband;
   Ymin=minband-0.2*space;
   Ymax=maxband+0.2*space;
   set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'YLim',[Ymin,Ymax],'FontName','Times New Roman');
      % top labels
      if count<=n
      title(endo{count,1},'FontWeight','normal');
      end
      % side labels
      if jj==1
      ylabel(Units{ii,1},'FontWeight','normal');
      end
      % legend
      if ii==1 && jj==1
      plotlegend=legend([ss,actual],'steady-state','actual');
      set(plotlegend,'FontName','Times New Roman');
      end
   end
end
end
% save on Excel
% create the cell that will be saved on excel
sscell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},T+3,1);
horzspace=repmat({''},3,6*n);
% loop over units
for ii=1:N
% initiate the cell of results
unitcell={};
   % loop over endogenous variables (horizontal dimension)
   for jj=1:n
   % create a header
   header=[{[Units{ii,1} ': ' endo{jj,1}]} {''} {''} {''} {''};{''} {''} {''} {''} {''};{''} {'actual'} {'lower bound'} {'median'} {'upper bound'}];
   % complete the cell
   endocell=[[header;stringdates1 num2cell(Yi(:,jj,ii)) num2cell(ss_estimates{jj,1}(1,:)') num2cell(ss_estimates{jj,1}(2,:)') num2cell(ss_estimates{jj,1}(3,:)')] vertspace];
   % concatenate to the previous parts of unitcell
   unitcell=[unitcell endocell];
   end
% concatenate to the previous parts of afcell
sscell=[sscell;horzspace;unitcell];
end
% trim
sscell=sscell(4:end,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),sscell,'steady state','B2');
end

author: european-central-bank

repo: BEAR-toolbox


text:
function [grad] = AssignCostGradFrame(grad, nFrOrig, nSeg, mask, Cost_layer)
D = size(grad,1);

words = ExtractWordsFromString_v2(Cost_layer.costFrameSelection);
selectionType = words{1};

switch selectionType
    case 'last'
        if strcmpi(Cost_layer.name, 'MSE')
            return;
        end

        if length(words)>1
            N = str2num(words{2});
        else
            N = 1;
        end

        gradTmp = grad;
        precision = class(gather(grad(1,1,1)));
        if strcmpi(class(grad), 'gpuArray')
            grad = gpuArray.zeros(D,nFrOrig, nSeg, precision);
        else
            grad = zeros(D,nFrOrig, nSeg, precision);
        end
        if numel(mask)>0 && sum(sum(mask))>0   % if the trajectories have variable length
            for i=1:nSeg
                idx = find(mask(:,i)==1);
                if isempty(idx)
                    grad(:,max(1,end-N+1):end,i) = gradTmp(:,:,i);
                else
                    grad(:,max(1,idx(1)-N):max(1,idx(1)-1),i) = gradTmp(:,:,i);
                end
            end
            grad = PadShortTrajectory(grad, mask, -1e10);
        else
            grad(:,max(1,end-N+1):end,:) = gradTmp;
        end
end
end


text:
function id = phWindowIndices(ph,phWindow)
%
% id = phWindowIndices(ph,phWindow)
%
% Returns indices of vector ph within phWindow
%
% if phWindow(1)<phWindow(2) returns phWindow(1) <= ph <= phWindow(2)
% if phWindow(1)>phWindow(2) returns ph >= phWindow(1) or ph <= phWindow(2)
%
% djh, 7/98

if diff(phWindow)>0
  id = find(ph>=phWindow(1) & ph<=phWindow(2));
else
  id = find(ph>=phWindow(1) | ph<=phWindow(2));
end

author: vistalab

repo: vistasoft


text:
function functionals2mrGray(volumeView, coordMask)
%
% functionals2mrGray(volumeView, [coordMask])
%
% AUTHOR:  S. Chial
% DATE:  02.05.99
% PURPOSE:
%   Writes out the functional data currently on display in the
%   volume window into an output file that can be read in by
%   mrGray (and mrUnfold) as an overlay file.
%
% volumeView: structure that contains all the information
% related to the volume functional data and user interface.
%
% sjc, 02/05/99
% rfd, 04/20/99 replaced the call to the now obsolete
%      "getOverlayClip" with a call to "getMapWindow".
% bw/hab 08/09/99
%      Introduced a subtraction of 1 from outLocs.  mrGray
%      locations run from [0 to n-1] and mrLoadRet run from [1 n].
%      This follows the change in readGrayGraph.  We abstracted
%      the conversion from mrLoadRet into mrGray with a new routine,
%      mrLoadRet2mrGray.
% rfd 10/15/99
%      Fixed a minor bug that produced bad cmap indices for some
%      clip ranges.  Basically, we just have to clip cmapIdx after we scale it.
% bw  01.16.01
%       Updated for mrLoadRet-2.5
%
% rfd 2002.08.08 added optional coordMask and updated to mrLoadRet-3
% bw  2003.08.28 added persistent recentDir, handled case of cancel on uiputfile

if(~exist('coordMask','var'))
    coordMask = [];
end
persistent recentDir      % Most recently visited directory

% Figure out which scan is being viewed.
%
scanNum = getCurScan(volumeView);

% Find the data that is above the correlation threshold
% and within the phase window.
%
coThreshold = getCothresh(volumeView);
validCoIdx = find(volumeView.co{scanNum} >= coThreshold);

phaseWindow = getPhWindow(volumeView);
validPhIdx = find((volumeView.ph{scanNum} >= phaseWindow(1)) & ...
    (volumeView.ph{scanNum} <= phaseWindow(2)));

validDataIdx = intersect(validCoIdx,validPhIdx);

% Extract the data being displayed and the corresponding locations
%
cmd = ['validData = volumeView.' volumeView.ui.displayMode '{scanNum}(validDataIdx);'];
eval(cmd);
outLocs = volumeView.coords(:,validDataIdx);
if(~isempty(coordMask))
    [outLocs, keepThese] = intersectCols(outLocs, coordMask);
else
    keepThese = [];
end

% Take mrLoadRet-2.0 locations (here outLocs) and convert them to
% mrGray locations.
%
outLocs = mrLoadRet2mrGray(outLocs);
%
% This is now done in mrLoadRet2mrGray:
%    outLocs = [outLocs(2,:); outLocs(1,:); outLocs(3,:)];

% Get the colormap for the data, convert it from having
% a range of 0 to 1 to a range of 0 to 255.
%
modeStr=['volumeView.ui.',volumeView.ui.displayMode,'Mode'];
mode = eval(modeStr);
cmap = round(255*mode.cmap);
numGrays = mode.numGrays;
numColors = mode.numColors;
clipMode = mode.clipMode;

% Only keep the part of the colormap that is color
%
cmap = cmap([numGrays+1:numGrays+numColors],:);

% Map the data to a color, clipping the colormap appropriately
%
overlayClip = getMapWindow(volumeView);
if ~isempty(validData)
    if strcmp(clipMode,'auto')
        minVal = min(validData);
        maxVal = max(validData);
        overClipMin = min(overlayClip)*(maxVal-minVal) + minVal;
        overClipMax = max(overlayClip)*(maxVal-minVal) + minVal;
    else
        overClipMin = min(clipMode);
        overClipMax = max(clipMode);
    end
    cmapIdx = round(scale_im(validData,1,numColors,overClipMin,overClipMax));
    % RFD: we need to clip, in case there are values above overClipMax or
    % under overClipMin.
    cmapIdx = round(clip(cmapIdx,1,numColors));
end
if(~isempty(keepThese))
    cmapIdx = cmapIdx(keepThese);
end

% Prompt user for a file name to save the data under
%
% This part of the routine should probably be a call to the new routine:
%
%    writeMrGray([outLocs; cmapIdx(:)'],cmap)
%
curDir = pwd;
if ~isempty(recentDir)
    if(exist(recentDir,'dir')) chdir(recentDir); end
else
    anatDir = getAnatomyPath([]);
    if(exist(anatDir,'dir')) chdir(anatDir); end;
end

[f p] = uiputfile('*.*','Save mrGray overlay file...');
if isequal(f,0) | isequal(p,0)
    return;
else
    chdir(curDir);
    filename = fullfile(p,f);
end

if ~isempty(filename)
    recentDir = p;  % Save the most recently visited directory.

    fid = fopen(filename,'w');

    % Write out the number of colors
    fprintf(fid,'%d\n',numColors);

    % Write out the colormap [R,G,B]
    fprintf(fid,'%d %d %d\n', cmap');

    % Write out the location and index into the colormap for
    % each data point.
    fprintf(fid,'%d %d %d %d\n', [outLocs; cmapIdx(:)']);

    fclose(fid);
    fprintf('Finished saving functional data overlay.\n');
end

return;


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function fmhypt
%FMHYPT	Unit test for the function FMHYP.

%	O. Lemoine - February 1996.

N=256;
t=1:N;

% For a hyperbolic instantaneous frequency, mentionning F0 and C
F0=0.05; C=.3;
[x,iflaw]=fmhyp(N,[F0,C]);
if any(abs(iflaw-(F0+C./t)')>sqrt(eps))~=0,
  error('fmhyp test 1 failed');
end

% For a hyperbolic instantaneous frequency, mentionning P1 and P2
P1=[1,.1]; P2=[200,.4];
[x,iflaw]=fmhyp(N,P1,P2);
if abs(iflaw(P1(1))-P1(2))>sqrt(eps) | abs(iflaw(P2(1))-P2(2))>sqrt(eps),
  error('fmhyp test 2 failed');
end

% For a hyperbolic group delay, mentionning P1 and P2
P1=[10,.45]; P2=[236,.25];
[x,iflaw]=fmhyp(N,P1,P2);
if abs(iflaw(P1(1))-P1(2))>sqrt(eps) | abs(iflaw(P2(1))-P2(2))>sqrt(eps),
  error('fmhyp test 3 failed');
end

N=251;
t=1:N;

% For a hyperbolic instantaneous frequency, mentionning F0 and C
F0=0.05; C=.3;
[x,iflaw]=fmhyp(N,[F0,C]);
if any(abs(iflaw-(F0+C./t)')>sqrt(eps))~=0,
  error('fmhyp test 4 failed');
end

% For a hyperbolic instantaneous frequency, mentionning P1 and P2
P1=[1,.1]; P2=[200,.4];
[x,iflaw]=fmhyp(N,P1,P2);
if abs(iflaw(P1(1))-P1(2))>sqrt(eps) | abs(iflaw(P2(1))-P2(2))>sqrt(eps),
  error('fmhyp test 5 failed');
end

% For a hyperbolic group delay, mentionning P1 and P2
P1=[10,.45]; P2=[236,.25];
[x,iflaw]=fmhyp(N,P1,P2);
if abs(iflaw(P1(1))-P1(2))>sqrt(eps) | abs(iflaw(P2(1))-P2(2))>sqrt(eps),
  error('fmhyp test 6 failed');
end

author: HeLiangHIT


text:
function datapt = ma_GACentralBodyTasks(stateLogEntry, subTask, celBodyData)
%ma_GACentralBodyTasks Summary of this function goes here
%   Detailed explanation goes here

    switch subTask
        case 'CB_ID'
            datapt = stateLogEntry(8);
    end
end


text:
function syntheticExperiment(method)

if (nargin<1)
    method = 'SDRSAC';
end

method = 'SDRSAC';
%method = '4PCS';
% Synthetic experiment

config = readConfig();
%config = readConfig('redwood');

load(config.matPath);

%outRANSAC = pointCloudReg(D, M, config, 'RANSAC');
%outICP = pointCloudReg(M, D, config, 'ICP');

% if (strcmp(method,'SDRSAC'))
%
%     D = normalize(D,2, 'center');
%     M = normalize(M,2, 'center');
%     %
%     D = normalize(D,2, 'range');
%     M = normalize(M,2, 'range');
% end

%out = pointCloudReg(D, M, config, method);
%plotPointClouds(D, M, 'r.','b.');
log_file_path = './logs/logs.txt';
for nrun = 1:10
    for m = 3
        config.k = m;
        for nsample = 16
            config.pointPerSample = nsample;
            out = pointCloudReg(D, M, config, method);
            disp(out);
            if strcmp(method,'SDRSAC')
                log_results(log_file_path, out, config);
            end
        end
    end
end

disp(out);

end

%outSDR = pointCloudReg(M, D, config, '4PCS');
%outGlobGM = pointCloudReg(M, D, config, 'GlobGM');
%[out] = SDRReg(D, M, config);

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function [EVec, Eload, EVal] = pca_calc(vecs);
%
%  PCA_CALC   Principal Component Analysis
%
%   [EVEC, ELOAD, EVAL] = PCA_CALC(X) takes a column-wise de-meaned
%       data matrix X of size n-by-m and calculates the n
%       Eigenvectors (EVEC of size n-by-n) of the data covariance
%       matrix, the factor loadings (ELOAD of size n-by-m) and the
%       corresponding Eigenvalues (EVAL of size n-by-1).
%
%
% Author: Rami K. Niazy, FMRIB Centre, University of Oxford
%
% Copyright (c) 2004 University of Oxford.
%

[m,n]=size(vecs);
[Us,S,EVec] = svd(vecs,0);

if m == 1
    S = S(1);
else
    S = diag(S);
end
Eload = Us .* repmat(S',m,1);
S = S ./ sqrt(m-1);
if m <= n
    S(m:n,1) = 0;
    S(:,m:n) = 0;
end
EVal = S.^2;

author: PatternRecognition

repo: OpenBMI


text:
function [M,behavior,stim] = GetTimeIndexedData(hfig,isAllCells) %#ok<INUSD>
%{
% naming convention used:
M = GetTimeIndexedData(hfig);
M_0 = GetTimeIndexedData(hfig,'isAllCells');
%}

isZscore = getappdata(hfig,'isZscore');
isTrialRes = getappdata(hfig,'isTrialRes');
isClusRes = getappdata(hfig,'isClusRes');

% main data input
if ~isZscore,
    cellResp = getappdata(hfig,'CellResp');
    cellRespAvr = getappdata(hfig,'CellRespAvr');
else
    cellResp = getappdata(hfig,'CellRespZ');
    cellRespAvr = getappdata(hfig,'CellRespAvrZ');
end

isMotorseed = getappdata(hfig,'isMotorseed');
if ~isMotorseed,
    Behavior_full = getappdata(hfig,'Behavior_full');
    BehaviorAvr = getappdata(hfig,'BehaviorAvr');
else
    Behavior_full = getappdata(hfig,'Behavior_full_motorseed');
    BehaviorAvr = getappdata(hfig,'BehaviorAvr_motorseed');
end
stim_full = getappdata(hfig,'stim_full');
stimAvr = getappdata(hfig,'stimAvr');
% other params
isStimAvr = getappdata(hfig,'isStimAvr');
cIX = getappdata(hfig,'cIX');
tIX = getappdata(hfig,'tIX');
% absIX = getappdata(hfig,'absIX');

%% set data
isFullData = getappdata(hfig,'isFullData');
if isFullData

    if isStimAvr,
        if exist('isAllCells','var'),
            M = cellRespAvr(:,tIX);
        else
            M = cellRespAvr(cIX,tIX);
        end
        behavior = BehaviorAvr(:,tIX);
        stim = stimAvr(:,tIX);
    else
        if exist('isAllCells','var'),
            M = cellResp(:,tIX);
        else
            M = cellResp(cIX,tIX);
        end
        behavior = Behavior_full(:,tIX);
        stim = stim_full(:,tIX);
    end

    if isTrialRes
        [~,M] = GetTrialAvrLongTrace(hfig,M);
        [~,behavior] = GetTrialAvrLongTrace(hfig,behavior);
        %     cellRespAvr = getappdata(hfig,'CellRespAvrZ');
    end

    if isClusRes
        gIX = getappdata(hfig,'gIX');
        [~,M] = FindClustermeans(gIX,M);
    end
else
    M = [];
    if isStimAvr,
        behavior = BehaviorAvr(:,tIX);
        stim = stimAvr(:,tIX);
    else
        behavior = Behavior_full(:,tIX);
        stim = stim_full(:,tIX);
    end
end

setappdata(hfig,'behavior',behavior);
setappdata(hfig,'stim',stim);
end

author: xiuyechen

repo: FishExplorer


text:
function [m_database V_PCA V_Fisher ProjectedImages_Fisher] = FisherfaceCore(T)
% Use Principle Component Analysis (PCA) and Fisher Linear Discriminant (FLD) to determine the most
% discriminating features between images of faces.
%
% Description: This function gets a 2D matrix, containing all training image vectors
% and returns 4 outputs which are extracted from training database.
% Suppose Ti is a training image, which has been reshaped into a 1D vector.
% Also, P is the total number of MxN training images and C is the number of
% classes. At first, centered Ti is mapped onto a (P-C) linear subspace by V_PCA
% transfer matrix: Zi = V_PCA * (Ti - m_database).
% Then, Zi is converted to Yi by projecting onto a (C-1) linear subspace, so that
% images of the same class (or person) move closer together and images of difference
% classes move further apart: Yi = V_Fisher' * Zi = V_Fisher' * V_PCA' * (Ti - m_database)
%
% Argument:      T                      - (M*NxP) A 2D matrix, containing all 1D image vectors.
%                                         All of 1D column vectors have the same length of M*N
%                                         and 'T' will be a MNxP 2D matrix.
%
% Returns:       m_database             - (M*Nx1) Mean of the training database
%                V_PCA                  - (M*Nx(P-C)) Eigen vectors of the covariance matrix of the
%                                         training database
%                V_Fisher               - ((P-C)x(C-1)) Largest (C-1) eigen vectors of matrix J = inv(Sw) * Sb
%                ProjectedImages_Fisher - ((C-1)xP) Training images, which are projected onto Fisher linear space
%
% See also: EIG

% Original version by Amir Hossein Omidvarnia, October 2007
%                     Email: aomidvar@ece.ut.ac.ir

Class_number = ( size(T,2) )/2; % Number of classes (or persons)
Class_population = 2; % Number of images in each class
P = Class_population * Class_number; % Total number of training images

%%%%%%%%%%%%%%%%%%%%%%%% calculating the mean image
m_database = mean(T,2);

%%%%%%%%%%%%%%%%%%%%%%%% Calculating the deviation of each image from mean image
A = T - repmat(m_database,1,P);

%%%%%%%%%%%%%%%%%%%%%%%% Snapshot method of Eigenface algorithm
L = A'*A; % L is the surrogate of covariance matrix C=A*A'.
[V D] = eig(L); % Diagonal elements of D are the eigenvalues for both L=A'*A and C=A*A'.

%%%%%%%%%%%%%%%%%%%%%%%% Sorting and eliminating small eigenvalues
L_eig_vec = [];
for i = 1 : P-Class_number
    L_eig_vec = [L_eig_vec V(:,i)];
end

%%%%%%%%%%%%%%%%%%%%%%%% Calculating the eigenvectors of covariance matrix 'C'
V_PCA = A * L_eig_vec; % A: centered image vectors

%%%%%%%%%%%%%%%%%%%%%%%% Projecting centered image vectors onto eigenspace
% Zi = V_PCA' * (Ti-m_database)
ProjectedImages_PCA = [];
for i = 1 : P
    temp = V_PCA'*A(:,i);
    ProjectedImages_PCA = [ProjectedImages_PCA temp];
end

%%%%%%%%%%%%%%%%%%%%%%%% Calculating the mean of each class in eigenspace
m_PCA = mean(ProjectedImages_PCA,2); % Total mean in eigenspace
m = zeros(P-Class_number,Class_number);
Sw = zeros(P-Class_number,P-Class_number); % Initialization os Within Scatter Matrix
Sb = zeros(P-Class_number,P-Class_number); % Initialization of Between Scatter Matrix

for i = 1 : Class_number
    m(:,i) = mean( ( ProjectedImages_PCA(:,((i-1)*Class_population+1):i*Class_population) ), 2 )';

    S  = zeros(P-Class_number,P-Class_number);
    for j = ( (i-1)*Class_population+1 ) : ( i*Class_population )
        S = S + (ProjectedImages_PCA(:,j)-m(:,i))*(ProjectedImages_PCA(:,j)-m(:,i))';
    end

    Sw = Sw + S; % Within Scatter Matrix
    Sb = Sb + (m(:,i)-m_PCA) * (m(:,i)-m_PCA)'; % Between Scatter Matrix
end

%%%%%%%%%%%%%%%%%%%%%%%% Calculating Fisher discriminant basis's
% We want to maximise the Between Scatter Matrix, while minimising the
% Within Scatter Matrix. Thus, a cost function J is defined, so that this condition is satisfied.
[J_eig_vec, J_eig_val] = eig(Sb,Sw); % Cost function J = inv(Sw) * Sb
J_eig_vec = fliplr(J_eig_vec);

%%%%%%%%%%%%%%%%%%%%%%%% Eliminating zero eigens and sorting in descend order
for i = 1 : Class_number-1
    V_Fisher(:,i) = J_eig_vec(:,i); % Largest (C-1) eigen vectors of matrix J
end

%%%%%%%%%%%%%%%%%%%%%%%% Projecting images onto Fisher linear space
% Yi = V_Fisher' * V_PCA' * (Ti - m_database)
for i = 1 : Class_number*Class_population
    ProjectedImages_Fisher(:,i) = V_Fisher' * ProjectedImages_PCA(:,i);
end


text:
function b = isPointInTriangle(point, p1, p2, p3)
%ISPOINTINTRIANGLE Test if a point is located inside a triangle
%
%   B = isPointInTriangle(POINT, V1, V2, V3)
%   POINT is a 1-by-2 row vector containing coordinates of the test point,
%   V1, V2 and V3 are 1-by-2 row vectors containing coordinates of triangle
%   vertices. The function returns 1 is the point is inside or on the
%   boundary of the triangle, and 0 otherwise.
%
%   B = isPointInTriangle(POINT, VERTICES)
%   Specifiy the coordinates of vertices as a 3-by-2 array.
%
%   If POINT contains more than one row, the result B has as many rows as
%   the input POINT.
%
%
%   Example
%     % vertices of the triangle
%     p1 = [0 0];
%     p2 = [10 0];
%     p3 = [5 10];
%     tri = [p1;p2;p3];
%     % check if points are inside
%     isPointInTriangle([0 0], tri)
%     ans =
%         1
%     isPointInTriangle([5 5], tri)
%     ans =
%         1
%     isPointInTriangle([10 5], tri)
%     ans =
%         0
%     % check for an array of points
%     isPointInTriangle([0 0;1 0;0 1], tri)
%     ans =
%         1
%         1
%         0
%
%   See also
%   polygons2d, isPointInPolygon, isCounterClockwise
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-05-16,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% if triangle vertices are given as a single array, extract vertices
if nargin == 2
    p2 = p1(2, :);
    p3 = p1(3, :);
    p1 = p1(1, :);
end

% check triangle orientation
isDirect = isCounterClockwise(p1, p2, p3);

% check location of point with respect to each side
if isDirect
    b12 = isCounterClockwise(p1, p2, point) >= 0;
    b23 = isCounterClockwise(p2, p3, point) >= 0;
    b31 = isCounterClockwise(p3, p1, point) >= 0;
else
    b12 = isCounterClockwise(p1, p2, point) <= 0;
    b23 = isCounterClockwise(p2, p3, point) <= 0;
    b31 = isCounterClockwise(p3, p1, point) <= 0;
end

% combines the 3 results
b = b12 & b23 & b31;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [data_interleave] = interleaving(data_convolutional,Ncpc,Tx);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                       %
%%     Name: interleaving.m                                              %
%%                                                                       %
%%     Description: This process is for interleave the bits to transmit  %
%%      and to make sure two consecutive bits are not modulated to       %
%%      the same frequency. It is another form to fight against          %
%%      frequency selective fading and burst errors.                     %
%%                                                                       %
%%                                                                       %
%%     Parameters:                                                       %
%%      Sequence of bits to interleave.                                  %
%%                                                                       %
%%     Result: It gives back the chain of bits, properly interleaved.    %
%%                                                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The size of the block to use to interleave the sequence is needed.
% This value is taken from table 223 of the standard.
% The value of the used block depends on the mapping that we want.

switch Ncpc
    case 1                          % for BPSK
        Ncbps=192;
    case 2                          % for QPSK
        Ncbps=384;
    case 4                          % for 16-QAM
        Ncbps=768;
    case 6                          % for 64-QAM
        Ncbps=1152;
end

s=ceil(Ncpc/2);

if Tx==1
  % In the following operations, it is necessary to consider the following values:
  % k-->index of the bit encoded BEFORE the first permutation
  % mk-->index of this bit BEFORE the second permutation and AFTER the first one
  % jk-->index after the SECOND permutation, just before the mapping of the sign.

    k = 0:Ncbps-1;
    mk = ((Ncbps/12)*mod(k,12))+floor(k/12);                     % First permutation
    jk = s*floor(mk/s)+mod(mk+Ncbps-floor(12*mk/Ncbps),s);       % Second permutation

    % Now I must arrange the indices to know in what order I must take the bits of the entering sequence.
    [a c] = sort(jk);

elseif Tx==0
  % Again, it is necessary to understand the values:
  % j-->index of the bit encoded BEFORE the first permutation
  % mj-->index of this bit BEFORE the second permutation and AFTER the first one
  % kj-->index after the SECOND permutation, just before the mapping of the sign.

    j = 0:Ncbps-1;
    mj = s*floor(j/s) + mod((j + floor(12*j/Ncbps)),s);          % First permutation
    kj = 12*mj-(Ncbps-1)*floor(12*mj/Ncbps);                     % Second permutation

    % The indices are ordered to know what must be taken.
    [a c]= sort(kj);

end

% finally the bits are rearranged.

 i = 1:Ncbps-1;
 data_interleave = zeros(1,Ncbps);
 data_interleave(i) = data_convolutional(c(i));