text:
% eegplotgold() - display EEG data in a clinical format
%
% Usage:
% >> eegplotgold('dataname', samplerate, 'chanfile', 'title', yscaling, range)
%
% Inputs:
%   'dataname' - quoted name of a desktop global variable (see Ex. below)
%   samplerate - EEG sampling rate in Hz (0 -> default 256 Hz)
%   'chanfile' - file of channel info in topoplot() style
%                                        (0 -> channel numbers)
%   'title'    - plot title string       (0 -> 'eegplotgold()')
%   yscaling   - initial y scaling factor (0 - default is 300)
%   range      - how many seconds to display in window (0 -> 10)
%
% Note: this version of eegplotgold() reguires that your data matrix
%       be defined as a global variable before running this routine.
%
% Example:  >> global dataname
%           >> eegplotgold('dataname')
%
% Author: Colin Humphries, CNL, Salk Institute, La Jolla, 3/97
%
% See also: eegplot(), eegplotold(), eegplotsold()

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

% 4-4-97 shortened name to eegplotgold() -sm
% 5-20-97 added read of icadefs.m for MAXEEGPLOTCHANS -sm
% 8-10-97 Clarified chanfile type -sm
% 01-25-02 reformated help & license, added links -ad

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

function y = eegplotold(dataname, samplerate, channamefile, titleval, yscaling, range)

eval (['global ',dataname])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define defaults
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set initial spacing

eval(['DEFAULT_SPACING = max(max(',dataname,''')-min(',dataname,'''));'])

%  spacing_var/20 = microvolts/millimeter with 21 channels
%  for n channels: 21/n * spacing_var/20 = microvolts/mm
%  for clinical data.

DEFAULT_SAMPLERATE = 256;	% default rate in Hz.
DEFAULT_PLOTTIME = 10;		% default 10 second window
DEFAULT_TITLE = 'eegplotgold()';
errorcode=0;			    % initialize error indicator

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Allow for different numbers of arguments
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 6
   PLOT_TIME = DEFAULT_PLOTTIME;
else
   PLOT_TIME = range;
end
if nargin < 5
   spacing_var = DEFAULT_SPACING;
else
   spacing_var = yscaling;
end
if spacing_var == 0
   spacing_var = DEFAULT_SPACING;
end

if nargin < 4
   titleval = DEFAULT_TITLE;
end
if titleval == 0
   titleval = DEFAULT_TITLE;
end
if nargin < 3
	channamefile = 0;
end
if nargin < 2
   samplerate = DEFAULT_SAMPLERATE;
end
if samplerate == 0,
	samplerate = DEFAULT_SAMPLERATE;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define internal variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

SAMPLE_RATE = samplerate;
time = 0;
eval(['[chans,frames] = size(',dataname,');'])		%size of data matrix

maxtime = frames / samplerate;       %size of matrix in seconds

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read the channel names
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if channamefile ~=0,		% read file of channel names
	chid = fopen(channamefile,'r');
	if chid <3,
		fprintf('plotdata: cannot open file %s.\n',channamefile);
		errorcode=2;
		channamefile = 0;
	else
		fprintf('Chan info file %s opened\n',channamefile);
	end;

    icadefs;   % read MAXEEGPLOTCHANS from icadefs.m
	if errorcode==0,
		channames = fscanf(chid,'%s',[6 MAXEEGPLOTCHANS]);
		channames = channames';
    	[r c] = size(channames);
		for i=1:r
			for j=1:c
				if channames(i,j)=='.',
					channames(i,j)=' ';
				end;
			end;
		end;
		% fprintf('%d channel names read from file.\n',r);
		if (r>chans)
			fprintf('Using first %d names.\n',chans);
			channames = channames(1:chans,:);
		end;
		if (r<chans)
			fprintf('Only %d channel names read.\n',r);
		end;
	end;
  end
  if channamefile ==0, % plot channel numbers
	channames = [];
	for c=1:chans
		if c<10,
			numeric = ['   ' int2str(c)];	% four-character fields
		else
			numeric = ['  '  int2str(c)];
		end
		channames = [channames;numeric];
	end;
  end; % setting channames

channames = str2mat(channames, ' ');	% add padding element to Y labels

Xlab = num2str(time);
for j = 1:1:PLOT_TIME
   Q = num2str(time+j);
   Xlab = str2mat(Xlab, Q);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set Graph Characteristics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   figure;				% plot a new figure
   fighandle = gcf;
   orient landscape		% choose landscape printer mode
   hold on;
   set(gcf,'NumberTitle','off')
   set(gcf,'Name',['EEGPLOTOLD #',num2str(gcf)])
   set (gca, 'xgrid', 'on')				%Xaxis gridlines only
   set (gca, 'GridLineStyle','-')			%Solid grid lines
   set (gca, 'XTickLabels', Xlab)			%Use Xlab for tick labels
   set (gca, 'Box', 'on')
   set (gca, 'XTick', time*samplerate:1.0*samplerate:PLOT_TIME*samplerate)
   set (gca, 'Ytick', 0:spacing_var:chans*spacing_var)  % ytickspacing on channels
   set (gca, 'TickLength', [0.001 0.001])
   title(titleval)					% title is titleval
   axis([0 PLOT_TIME*samplerate 0 (chans+1)*spacing_var]);       % set axis values
   set (gca, 'YTickLabels', flipud(channames))   	% write channel names

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the selected EEG data epoch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:chans
   if (maxtime-time>PLOT_TIME)
      eval(['F = ',dataname,'(chans-i+1,(time*samplerate)+1:(time+PLOT_TIME*samplerate));'])
   else
      eval(['F = ',dataname,'(chans-i+1,(time*samplerate)+1:(maxtime*samplerate));'])
   end
   F = F - mean(F) + i*spacing_var;
   plot (F,'clipping','off')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot Scaling I
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

line([(PLOT_TIME+.3)*samplerate,(PLOT_TIME+.3)*samplerate],[1.5*spacing_var 2.5*spacing_var],'clipping','off','color','w')
line([(PLOT_TIME+.3)*samplerate-10,(PLOT_TIME+.3)*samplerate+10],[2.5*spacing_var,2.5*spacing_var],'clipping','off','color','w')
line([(PLOT_TIME+.3)*samplerate-10,(PLOT_TIME+.3)*samplerate+10],[1.5*spacing_var,1.5*spacing_var],'clipping','off','color','w')
text((PLOT_TIME+.5)*samplerate,2*spacing_var,num2str(round(spacing_var)),'clipping','off')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% User Control Routines
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   slider_position = [.125 .030 .3 .024];   % position of user-controlled slider
   edit_position = [.65 .025 .1 .05];       % position of edit box
   slider_position2 = [.8 .03 .1 .024];
   b1_position = [.175 .022 .09 .047];
   b2_position = [.29 .022 .09 .047];
   b3_position = [.125 .022 .045 .047];
   b4_position = [.385 .022 .045 .047];

   Max_Space = 1;
   Min_Space = DEFAULT_SPACING*2;
%   User_Data_Mat = [data;zeros(1,length(data))];

axhandle = gca;
   User_Data_Mat(1) = samplerate;
   User_Data_Mat(2) = PLOT_TIME;
   User_Data_Mat(3) = spacing_var;
   User_Data_Mat(4) = time;
   User_Data_Mat(5) = maxtime;
   User_Data_Mat(6) = axhandle;
   User_Data_Mat(7) = 1;  % color
   User_Data_Mat(8) = frames;
   User_Data_Mat(9) = chans;
   User_Data_Mat(12) = 1;

   tstring1 = 'data1973 = get(gcf,''UserData'');';
   tstring2 = 'set(gcf,''UserData'',data1973);';
   tstring3 = 'eegdrawgv(gcf);';

   TIMESTRING = [tstring1,'if (data1973(4)-data1973(2))<0;','data1973(4) = 0;','else;','data1973(4) = data1973(4) - data1973(2);','end;',tstring2,tstring3,'clear data1973'];

   hb = uicontrol('Style','PushButton','Units','Normalized','position',b1_position,'String','PREV','Callback',TIMESTRING);

   TIMESTRING = [tstring1,'if (data1973(4)+data1973(2))>=data1973(5);','data1973(4)=data1973(4);','else;','data1973(4) = data1973(4) + data1973(2);','end;',tstring2,tstring3,'clear data1973'];

   hf = uicontrol('Style','PushButton','Units','Normalized','position',b2_position,'String','NEXT','Callback',TIMESTRING);

   TIMESTRING = [tstring1,'if (data1973(4)-1)<0;','data1973(4) = 0;','else;','data1973(4) = data1973(4) - 1;','end;',tstring2,tstring3,'clear data1973'];

   hbos = uicontrol('Style','PushButton','Units','Normalized','position',b3_position,'String','<<','Callback',TIMESTRING);

   TIMESTRING = [tstring1,'if (data1973(4)+1)>=data1973(5);','data1973(4)=data1973(4);','else;','data1973(4) = data1973(4) + 1;','end;',tstring2,tstring3,'clear data1973'];

   hfos = uicontrol('Style','PushButton','Units','Normalized','position',b4_position,'String','>>','Callback',TIMESTRING);

   TIMESTRING = [tstring1,'time1973 = get(gco,''string'');','time1973 = str2num(time1973);','data1973(4) = time1973;',tstring2,tstring3,'clear time1973 data1973'];

   w=uicontrol('style','edit','units','normalized','HorizontalAlignment','left','position',edit_position,'UserData',axhandle,'callback',TIMESTRING);

TIMESTRING = [tstring1,'data1973(3) = get(gco,''value'');',tstring2,tstring3,'clear time1973 data1973'];

u=uicontrol('style','slider','units','normalized','position',slider_position2,'Max',Max_Space,'Min',Min_Space,'value',spacing_var,'UserData',axhandle,'callback',TIMESTRING);

%%%%%%%%%%%%%%%%%%%%%%%%%
%Set up ui menus
%%%%%%%%%%%%%%%%%%%%%%%%%

%Window menu:

TIMESTRING = ['fighand1973 = gcf;','delete(fighand1973);','clear fighand1973;'];

fm1 = uimenu('Label','Window');
fm2 = uimenu(fm1,'Label','Close ','UserData',fighandle,'Callback',TIMESTRING);

%Display menu:

TIMESTRING = [tstring1,'out1973 = gettext(''Input new windowlength(sec).'');','if isempty(out1973);','out1973 = 0;','else;','data1973(2) = str2num(out1973);',tstring2,tstring3,'end;','clear data1973 out1973'];

dm1 = uimenu('Label','Display');
dm2 = uimenu(dm1,'Label',' Window Length','Interruptible','yes','Callback',TIMESTRING);
dm3 = uimenu(dm1,'Label',' Color');

TIMESTRING = [tstring1,'data1973(7) = 1;',tstring2,tstring3,'clear data1973;'];

dm4 = uimenu(dm3,'Label','Yellow ','UserData',axhandle,'Interruptible','yes','Callback',TIMESTRING);

TIMESTRING = [tstring1,'data1973(7) = 2;',tstring2,tstring3,'clear data1973;'];

dm5 = uimenu(dm3,'Label','White ','UserData',axhandle,'Interruptible','yes','Callback',TIMESTRING);

TIMESTRING = ['label1973 = gettext(''Enter new title.'');','if isempty(label1973);','label1973 = 0;','else;','title(label1973);','end;','clear label1973;'];

dm6 = uimenu(dm1,'Label','Title ','Interruptible','yes','Callback',TIMESTRING);

TIMESTRING = [tstring1,'Check1973 = get(data1973(10),''checked'');','if (Check1973(1:2) == ''on'');','set(data1973(10),''Checked'',''off'');','set(data1973(6),''XGrid'',''off'');','else;','set(data1973(10),''Checked'',''on'');','set(data1973(6),''XGrid'',''on'');','end;','clear data1973 Check1973;'];

dm7 = uimenu(dm1,'Label','Grid','Checked','on','Callback',TIMESTRING);

User_Data_Mat(10) = dm7;

TIMESTRING = [tstring1,'Check1973 = get(data1973(11),''checked'');','if (Check1973(1:2) == ''on'');','set(data1973(11),''Checked'',''off'');','data1973(12)= 0;','else;','set(data1973(11),''Checked'',''on'');','data1973(12) = 1;','end;',tstring2,tstring3,'clear data1973 Check1973;'];

dm8 = uimenu(dm1,'Label','Scaling I','Checked','on','Callback',TIMESTRING);

User_Data_Mat(11) = dm8;

%Settings menu:

sm1 = uimenu('Label','Settings');

TIMESTRING = [tstring1,'Srate1973 = gettext(''Enter new samplerate'');','if isempty(Srate1973);','Srate1973 = 0;','else;','data1973(1) = str2num(Srate1973);','data1973(5) = data1973(8)/data1973(1);','data1973(4) = 0;',tstring2,tstring3,'end;','clear Srate1973 data1973'];

sm2 = uimenu(sm1,'Label','Samplerate','Interruptible','yes','Callback',TIMESTRING);

%Electrodes menu:

em1 = uimenu('Label','Electrodes');

TIMESTRING = [tstring1,'ChanNamefile1973 = gettext(''Enter Electrode file to load.'');','if isempty(ChanNamefile1973);','ChanNamefile1973=0;','else;','ChanNames1973 = loadelec(ChanNamefile1973);','set(data1973(6),''YTickLabels'',flipud(ChanNames1973));','end;','clear data1973 ChanNamefile1973 ChanNames1973'];

em2 = uimenu(em1,'Label','Load Electrode File ','Interruptible','yes','Callback',TIMESTRING);

TIMESTRING = [tstring1,'ChanNames1973 = makeelec(data1973(9));','if isempty(ChanNames1973);','ChanNames1973=0;','else;','set(data1973(6),''YTickLabels'',flipud(ChanNames1973));','end;','clear data1973 ChanNames1973'];

em3 = uimenu(em1,'Label','Make Electrode File ','Interruptible','yes','Callback',TIMESTRING);

set(axhandle,'UserData',dataname)
set(fighandle,'UserData',User_Data_Mat)

author: PatternRecognition

repo: OpenBMI


text:
function [po,pin] = spm_prep2sn(p)
% Convert the output from spm_preproc into an sn.mat file
% FORMAT [po,pin] = spm_prep2sn(p)
% p   - the results of spm_preproc
%
% po  - the output in a form that can be used by spm_write_sn
% pin - the inverse transform in a form that can be used by spm_write_sn
%
% The outputs are saved in sn.mat files only if they are not requested LHS.
%__________________________________________________________________________
% Copyright (C) 2005-2011 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_prep2sn.m 4873 2012-08-30 19:06:26Z john $

if ischar(p), p = load(p); end

VG          = p.tpm;
VF          = p.image;
[Y1,Y2,Y3]  = create_def(p.Twarp,VF,VG(1),p.Affine);
[Y1,Y2,Y3]  = spm_invdef(Y1,Y2,Y3,VG(1).dim(1:3),eye(4),eye(4));
MT          = procrustes(Y1,Y2,Y3,VG(1),VF.mat);
d           = size(p.Twarp);
[Affine,Tr] = reparameterise(Y1,Y2,Y3,VG(1),VF.mat,MT,max(d(1:3)+2,[8 8 8]));
flags       = struct(...
    'ngaus',    p.ngaus,...
    'mg',       p.mg,...
    'mn',       p.mn,...
    'vr',       p.vr,...
    'warpreg',  p.warpreg,...
    'warpco',   p.warpco,...
    'biasreg',  p.biasreg,...
    'biasfwhm', p.biasfwhm,...
    'regtype',  p.regtype,...
    'fudge',    p.fudge,...
    'samp',     p.samp,...
    'msk',      p.msk,...
    'Affine',   p.Affine,...
    'Twarp',    p.Twarp,...
    'Tbias',    p.Tbias,...
    'thresh',   p.thresh);

% Parameterisation for the forward transform
po          = struct(...
    'VG',       VG,...
    'VF',       VF,...
    'Tr',       Tr,...
    'Affine',   Affine,...
    'flags',    flags);

% Parameterisation for the inverse transform
pin         = struct(...
    'VG',       p.image,...
    'VF',       p.tpm(1),...
    'Tr',       p.Twarp,...
    'Affine',   p.tpm(1).mat\p.Affine*p.image.mat,...
    'flags',    flags);

if ~nargout
    [pth,nam] = spm_fileparts(VF.fname);

    fnam_out  = fullfile(pth,[nam '_seg_sn.mat']);
    fnam_inv  = fullfile(pth,[nam '_seg_inv_sn.mat']);
    save(fnam_out,'-struct','po', spm_get_defaults('mat.format'));
    save(fnam_inv,'-struct','pin', spm_get_defaults('mat.format'));
end
return;
%==========================================================================

%==========================================================================
function [Affine,Tr] = reparameterise(Y1,Y2,Y3,B,M2,MT,d2)
% Take a deformation field and reparameterise in the same form
% as used by the spatial normalisation routines of SPM
d          = [size(Y1) 1];
[x1,x2,o]  = ndgrid(1:d(1),1:d(2),1);
x3         = 1:d(3);
Affine     = M2\MT*B(1).mat;
A          = inv(Affine);

B1  = spm_dctmtx(d(1),d2(1));
B2  = spm_dctmtx(d(2),d2(2));
B3  = spm_dctmtx(d(3),d2(3));
pd  = prod(d2(1:3));
AA  = eye(pd)*0.01;
Ab  = zeros(pd,3);
spm_progress_bar('init',length(x3),'Reparameterising','Planes completed');
for z=1:length(x3)
    y1       = double(Y1(:,:,z));
    y2       = double(Y2(:,:,z));
    y3       = double(Y3(:,:,z));
    msk      = isfinite(y1);
    w        = double(msk);
    y1(~msk) = 0;
    y2(~msk) = 0;
    y3(~msk) = 0;
    z1       = A(1,1)*y1+A(1,2)*y2+A(1,3)*y3 + w.*(A(1,4) - x1);
    z2       = A(2,1)*y1+A(2,2)*y2+A(2,3)*y3 + w.*(A(2,4) - x2);
    z3       = A(3,1)*y1+A(3,2)*y2+A(3,3)*y3 + w *(A(3,4) - z );
    b3       = B3(z,:)';
    Ab(:,1)  = Ab(:,1) + kron(b3,spm_krutil(z1,B1,B2,0));
    Ab(:,2)  = Ab(:,2) + kron(b3,spm_krutil(z2,B1,B2,0));
    Ab(:,3)  = Ab(:,3) + kron(b3,spm_krutil(z3,B1,B2,0));
    AA       = AA  + kron(b3*b3',spm_krutil(w, B1,B2,1));
    spm_progress_bar('set',z);
end
spm_progress_bar('clear');
Tr  = reshape(AA\Ab,[d2(1:3) 3]);
return;
%==========================================================================

%==========================================================================
function MT = procrustes(Y1,Y2,Y3,B,M2)
% Take a deformation field and determine the closest rigid-body
% transform to match it, with weighing.
%
% Example Reference:
% F. L. Bookstein (1997).  "Landmark Methods for Forms Without
% Landmarks: Morphometrics of Group Differences in Outline Shape"
% Medical Image Analysis 1(3):225-243

M1        = B.mat;
d         = B.dim(1:3);
[x1,x2,o] = ndgrid(1:d(1),1:d(2),1);
x3        = 1:d(3);
c1        = [0 0 0];
c2        = [0 0 0];
sw        =  0;
spm_progress_bar('init',length(x3),'Procrustes (1)','Planes completed');
for z=1:length(x3)
    y1    = double(Y1(:,:,z));
    y2    = double(Y2(:,:,z));
    y3    = double(Y3(:,:,z));
    msk   = find(isfinite(y1));
    w     = spm_sample_vol(B(1),x1(msk),x2(msk),o(msk)*z,0);
    swz   = sum(w(:));
    sw    = sw+swz;
    c1    = c1 + [w'*[x1(msk) x2(msk)] swz*z ];
    c2    = c2 +  w'*[y1(msk) y2(msk) y3(msk)];
    spm_progress_bar('set',z);
end
spm_progress_bar('clear');
c1        = c1/sw;
c2        = c2/sw;
T1        = [eye(4,3) M1*[c1 1]'];
T2        = [eye(4,3) M2*[c2 1]'];
C         = zeros(3);
spm_progress_bar('init',length(x3),'Procrustes (2)','Planes completed');
for z=1:length(x3)
    y1    = double(Y1(:,:,z));
    y2    = double(Y2(:,:,z));
    y3    = double(Y3(:,:,z));
    msk   = find(isfinite(y1));
    w     = spm_sample_vol(B(1),x1(msk),x2(msk),o(msk)*z,0);
    C = C + [(x1(msk)-c1(1)).*w (x2(msk)-c1(2)).*w (    z-c1(3))*w ]' * ...
            [(y1(msk)-c2(1))    (y2(msk)-c2(2))    (y3(msk)-c2(3)) ];
    spm_progress_bar('set',z);
end
spm_progress_bar('clear');
[u,s,v]    = svd(M1(1:3,1:3)*C*M2(1:3,1:3)');
R          = eye(4);
R(1:3,1:3) = v*u';
MT         = T2*R*inv(T1);
return;
%==========================================================================

%==========================================================================
function [Y1,Y2,Y3] = create_def(T,VG,VF,Affine)
% Generate a deformation field from its parameterisation.
d2   = size(T);
d    = VG.dim(1:3);
M    = VF.mat\Affine*VG.mat;
[x1,x2,o] = ndgrid(1:d(1),1:d(2),1);
x3   = 1:d(3);
B1   = spm_dctmtx(d(1),d2(1));
B2   = spm_dctmtx(d(2),d2(2));
B3   = spm_dctmtx(d(3),d2(3));
[pth,nam] = spm_fileparts(VG.fname);
spm_progress_bar('init',length(x3),['Creating Def: ' nam],'Planes completed');
for z=1:length(x3)
    [y1,y2,y3] = defs(T,z,B1,B2,B3,x1,x2,x3,M);
    Y1(:,:,z)  = single(y1);
    Y2(:,:,z)  = single(y2);
    Y3(:,:,z)  = single(y3);
    spm_progress_bar('set',z);
end
spm_progress_bar('clear');
return;
%==========================================================================

%==========================================================================
function [x1,y1,z1] = defs(sol,z,B1,B2,B3,x0,y0,z0,M)
if ~isempty(sol)
    x1a = x0    + transf(B1,B2,B3(z,:),sol(:,:,:,1));
    y1a = y0    + transf(B1,B2,B3(z,:),sol(:,:,:,2));
    z1a = z0(z) + transf(B1,B2,B3(z,:),sol(:,:,:,3));
else
    x1a = x0;
    y1a = y0;
    z1a = z0;
end
x1  = M(1,1)*x1a + M(1,2)*y1a + M(1,3)*z1a + M(1,4);
y1  = M(2,1)*x1a + M(2,2)*y1a + M(2,3)*z1a + M(2,4);
z1  = M(3,1)*x1a + M(3,2)*y1a + M(3,3)*z1a + M(3,4);
return;
%==========================================================================

%==========================================================================
function t = transf(B1,B2,B3,T)
d2 = [size(T) 1];
t1 = reshape(T, d2(1)*d2(2),d2(3));
t1 = reshape(t1*B3', d2(1), d2(2));
t  = B1*t1*B2';
return;
%==========================================================================

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = mpower_internal(varargin)
switch class(varargin{1})
    case 'double'
        varargout{1} = varargin{1};
    case 'char'
        % Note, the order here is assumed in GP formulation
        varargout{1} = (varargin{3} == varargin{2});
        varargout{2} = struct('convexity','none','monotonicity','none','definiteness','none','model','exact');
        varargout{3} = varargin{3};
    otherwise
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [U_final, V_final, nIter_final, objhistory_final] = LCCF(X, k, W, options, U, V)
% Locally Consistant Concept Factorization (LCCF)
%
% where
%   X
% Notation:
% X ... (mFea x nSmp) data matrix
%       mFea  ... number of words (vocabulary size)
%       nSmp  ... number of documents
% k ... number of hidden factors
% W ... weight matrix of the affinity graph
%
% options ... Structure holding all settings
%
% You only need to provide the above four inputs.
%
% X = X*U*V'
%
% References:
% [1] Deng Cai, Xiaofei He, Jiawei Han, "Locally Consistent Concept
%     Factorization for Document Clustering", IEEE Transactions on Knowledge
%     and Data Engineering, Vol. 23, No. 6, pp. 902-913, 2011.
%
%
%   version 2.0 --April/2010
%   version 1.0 --Dec./2008
%
%   Written by Deng Cai (dengcai AT gmail.com)
%

nSmp=size(X,2);

if ~isfield(options,'error')
    options.error = 1e-5;
end
if ~isfield(options, 'maxIter')
    options.maxIter = [];
end

if ~isfield(options,'nRepeat')
    options.nRepeat = 10;
end

if ~isfield(options,'minIter')
    options.minIter = 30;
end

if ~isfield(options,'meanFitRatio')
    options.meanFitRatio = 0.1;
end

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

if isfield(options,'alpha_nSmp') && options.alpha_nSmp
    options.alpha = options.alpha*nSmp;
end

if ~exist('U','var')
    U = [];
    V = [];
end

K = constructKernel(X',[],options);

if isfield(options,'weight') && strcmpi(options.weight,'NCW')
    D_mhalf = sum(K,2).^-.5;
    D_mhalf = spdiags(D_mhalf,0,nSmp,nSmp);
    K = D_mhalf*K*D_mhalf;
end

if ~isfield(options,'Optimization')
    options.Optimization = 'Multiplicative';
end

switch lower(options.Optimization)
    case {lower('Multiplicative')}
        [U_final, V_final, nIter_final,objhistory_final] = LCCF_Multi(K, k, W, options, U, V);
    otherwise
        error('optimization method does not exist!');
end

author: ZJULearning

repo: MatlabFunc


text:
function [ b, m1, m2, r11, r12, r21, r22, seed ] = p36_p_get (  )

%*****************************************************************************80
%
%% P36_P_GET gets the values of the parameters for problem 36.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real B, M1, M2, R11, R12, R21, R22,
%    problem parameters.
%
%    Output, integer SEED, a seed for the random number generator.
%
  b = 0.0;
  m1 = 0;
  m2 = 0;
  r11 = 0.0;
  r12 = 0.0;
  r21 = 0.0;
  r22 = 0.0;
  seed = 0;

  [ b, m1, m2, r11, r12, r21, r22, seed ] = ...
    p36_p_val ( 'GET', b, m1, m2, r11, r12, r21, r22, seed );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [results_table, results_table_pos, results_table_neg, table_legend_text] = table_simple(r)
% Make a simple table of fMRI results with one row per region
%
% [results_table, results_table_pos, results_table_neg, table_legend_text] = table_simple(r)
%
% - This is most useful if you've divided an activation map into blobs (regions)
% separated by defined anatomical parcels.
% - e.g., see image_vector.subdivide_by_atlas
%
% see also region.table, for autolabeling of regions with an atlas and more
% detailed table output.
%
% Tor Wager, April 2021

n_cols = 140;                       % 140 good for HTML reports
sep_str = repmat('_', 1, n_cols);   % see textwrap

dosep = true;
dolegend = true;
table_legend_text = '';

[results_table_pos, results_table_neg] = deal(table());

if isempty(r) || length(r) == 0
    disp('No regions to display');
    return
end

% Separate subregions with positive and negative values if requested
% -------------------------------------------------------------------------
if dosep
    % separate pos and neg
    [poscl, negcl] = posneg_separate(r);

    r = [poscl negcl];
    ispos = [true(1, length(poscl)) false(1, length(negcl))]; % logical for splitting combined cl later

    clear poscl negcl

    fprintf('\n%s\nPositive Effects\n', sep_str)

else

    fprintf('\n%s\nTable of all regions\n', sep_str)
end

% Volume in mm^3
% count per cubic mm - voxel count * voxel volume
get_vox_volume_in_mm = @(r) size(r.XYZ, 2) .* prod(abs(diag(r.M(1:3, 1:3))));

Volume = zeros(length(r), 1);
for i = 1:length(r)

    Volume(i, 1) = get_vox_volume_in_mm(r(i));

end

results_table = table({r.shorttitle}', 'VariableNames', {'Region'});

results_table.Volume = Volume;

results_table.XYZ = round(cat(1, r.mm_center));

results_table = [results_table get_signed_max(r, 'Z', 'maxZ')];  % use function because may be empty, handle if so

results_table.Properties.VariableNames{4} = 'maxVal';

if dosep

    results_table_pos = results_table(ispos, :);
    results_table_neg = results_table(~ispos, :);

end

% Legend info

voxvol = prod(abs(diag(r(1).M(1:3, 1:3))));

table_legend_text = {'Note: XYZ values are Montreal Neurologic Institute coordinates in mm.'};
table_legend_text{2} = sprintf('Voxel volume is %3.2fmm^3', voxvol);
table_legend_text{3} = 'maxVal is maximum positive or negatively-valued value in the .Z field. Its meaning varies across analyses. For a t-test, it usually reflects the maximum t-value.';

%
%     % Sort, if asked for (default = yes)
%     if dosortrows
%
%         % Replace empty strings so sort will work
%         whempty = cellfun(@isempty, results_table_pos.modal_label_descriptions);
%         results_table_pos.modal_label_descriptions(whempty) = {'X_No_label'};
%
%         whempty = cellfun(@isempty, results_table_neg.modal_label_descriptions);
%         results_table_neg.modal_label_descriptions(whempty) = {'X_No_label'};
%
%         results_table_pos = sortrows(results_table_pos, 'modal_label_descriptions');
%         results_table_neg = sortrows(results_table_neg, 'modal_label_descriptions');
%
%     end
%

% Now split into positive and neg sub-tables and display
if dosep
    % Tables for pos and neg results separately

    if any(ispos)
        disp(results_table_pos)
    else
        disp('No regions to display');
    end

    fprintf('\nNegative Effects\n')
    if any(~ispos)
        disp(results_table_neg)
    else
        disp('No regions to display');
    end

else
    % No separation

    if size(results_table, 1) > 0
        disp(results_table)
    else
        disp('No regions to display');
    end

end % dosep

if dolegend == false || isempty(table_legend_text)

    return

else
    % Print legend text

    canlab_print_legend_text(table_legend_text{:});
end

% clean up text for legend
% if length(table_legend_text) > 2
%     table_legend_text(2:3) = [];
% end

end % main function

function val_table = get_signed_max(cl, myfield, tablevarname)
% Returns a table with var "MaxZ" nregions x 1, or empty if cl.Z is empty

%maxZ = @(i) max(double(cl(i).Z));

smax = @(i) signedmax(cl(i).(myfield));

for i = 1:length(cl)

    if ~isempty(cl(i).(myfield))
        myZ(i, 1) = smax(i);
    else
        myZ(i, 1) = NaN;
    end

    %if isinf(smax(i)), keyboard, end

end

% Fix infinite vals - only for .Z . so this is not generalizable beyond
% this function without modifications:
maxZ = norminv(1 - 1E-12);
myZ(myZ > maxZ) = maxZ;

if all(isnan(myZ))
    val_table = [];
else
    val_table = table(myZ,  'VariableNames', {tablevarname});
end

end % function

function val = signedmax(vals)

vals = double(vals);
[maxabs, wh] = max(abs(vals));

val = sign(vals(wh)) .* maxabs;

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function A = meanZero(hyp, x, i)

% Zero mean function. The mean function does not have any parameters.
%
% m(x) = 0
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2010-01-10.
%
% See also MEANFUNCTIONS.M.

if nargin<2, A = '0'; return; end             % report number of hyperparameters
A = zeros(size(x,1),1);                                    % derivative and mean

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [ coor, tri ] = getUniformVector( level )
%GETUNIFORMVECTOR Summary of this function goes here
%   Detailed explanation goes here
if exist(sprintf('./BasicFuncPano/icosahedron2sphere/uniformvector_lvl%d.mat', level), 'file')
    load(sprintf('./BasicFuncPano/icosahedron2sphere/uniformvector_lvl%d.mat', level));
else
    [coor,tri] = icosahedron2sphere(level);
end

end

author: zouchuhang

repo: LayoutNet


text:
% computes compactness of ROIs and fits ellipses
function stat = anatomize(ops, mPix, mLam, stat)

di = ops.diameter;
d0   = ceil(ops.diameter); % expected cell diameter

% data Ucell is nMaps by Ly by Lx

dx = repmat([-d0:d0], 2*d0+1, 1);
dy = dx';

rs = dx.^2 + dy.^2;
dx = dx(rs<=d0^2);
dy = dy(rs<=d0^2);

d2p = (bsxfun(@minus, dx, dx').^2 + bsxfun(@minus, dy, dy').^2).^.5;

xlx         = repmat(-ceil(2*d0):1:ceil(2*d0), 2*ceil(2*d0)+1, 1);
rgrid       = sqrt(xlx.^2 + xlx'.^2);
[rgridsort, isort]  = sort(rgrid(:), 'ascend');
xlxt        = xlx';

d2p0 = (bsxfun(@minus, xlx(:), xlx(:)').^2 + bsxfun(@minus, xlxt(:), xlxt(:)').^2).^.5;
d2p0 = d2p0(isort, isort);

%%
rd = zeros(size(mPix,2), 1);
rd0 = zeros(size(mPix,2), 1);

%% compute compactness and aspect ratio
for j = 1:size(mPix,2)

    lam  = mLam(:,j);

    gpix = lam>1e-3;

    dd = d2p(gpix, gpix);

    stat(j).mrs(1) = mean(dd(:))/d0;
    dd = d2p0(1:sum(gpix), 1:sum(gpix));
    stat(j).mrs0(1) = mean(dd(:))/di;

    stat(j).cmpct = stat(j).mrs(1)/stat(j).mrs0(1);

    params = FitMVGaus(stat(j).ypix,stat(j).xpix, stat(j).lam);

    % save ellipse information
    stat(j).aspect_ratio = sqrt(max(params.eval) / min(params.eval));
    stat(j).ellipse      = params.xy;

end

%%

author: cortex-lab


text:
function Obj = SOFAremoveVariable(Obj,Name)
%SOFAremoveVariable - Remove a variable, its dimension, and its attributes from a SOFA object
%  Usage: Obj = SOFAremoveVariable(Obj,Name)
%
%  SOFAremoveVariable(Obj,Name) removes the variable Name
%  from the SOFA structure Obj.
%

% #Author: Piotr Majdak: adapted from SOFAaddVariable (19.06.2019)
% #Author: Michael Mihocic: header documentation updated (28.10.2021)
% #Author: Michael Mihocic: removing variable attributes; minor bug fixed (20.09.2022)
%
% SOFA Toolbox - function SOFAremoveVariable
% Copyright (C) Acoustics Research Institute - Austrian Academy of Sciences
% Licensed under the EUPL, Version 1.2 or - as soon they will be approved by the European Commission - subsequent versions of the EUPL (the "License")
% You may not use this work except in compliance with the License.
% You may obtain a copy of the License at: https://joinup.ec.europa.eu/software/page/eupl
% Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing  permissions and limitations under the License.

switch Name
    case {'API','PRIVATE','GLOBAL','Data'}
        error('This variable name is reserved.');
    otherwise
        if isfield(Obj,Name)
            Obj=rmfield(Obj,Name);
            if isfield(Obj.API.Dimensions,Name)
                Obj.API.Dimensions=rmfield(Obj.API.Dimensions,Name);
            end
            % check if fields "Name_*" are existing and remove them (variable attributes)
            ObjFields=fieldnames(Obj);

            if exist('OCTAVE_VERSION','builtin')
              % We're in Octave
              IndexC = strfind(ObjFields,[Name '_']);
              IndexContainsName = find(not(cellfun('isempty',IndexC))); % do not replace by 'contains' (not supported in Octave)
            else
              % We're in Matlab
              IndexContainsName = find(contains(ObjFields,[Name '_']));
            end

            for ii=length(IndexContainsName):-1:1 % loop is not entered if empty IndexContainsName; run backwards because of removed indices
                % check if field name starts with "Name_*"
                if startsWith(char(ObjFields(IndexContainsName(ii))),Name)
%                     disp(['removing field: ' char(ObjFields(IndexContainsName(ii)))]) % debug
                    Obj=rmfield(Obj,ObjFields(IndexContainsName(ii)));
                end
            end

        end
end

author: sofacoustics

repo: SOFAtoolbox


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

%Fourier Transfrom properties

%	Time shifting

%x(t)=cos(t)
syms t w
x=cos(t);
t0=2;
xt0=cos(t-t0);
Left=fourier(xt0,w)
X=fourier(x,w);

Right=exp(-j*w*t0)*X

author: Sable

repo: mcbench-benchmarks


text:
function G = createGroupIncidenceMatrix_old(model, gcmOutputFile, gcmMetList, jankowskiGroupData)
% Creates `groupData` struct to calculate reaction Gibbs energies with reduced
% error in `vonB`.
%
% USAGE:
%
%    G = createGroupIncidenceMatrix_old(model, gcmOutputFile, gcmMetList, jankowskiGroupData)
%
% INPUTS:
%    model:
%    gcmOutputFile:
%    gcmMetList:
%    jankowskiGroupData:
%
% OUTPUT:
%    G:

groups = model.jankowskiGroupData.groups;
metList = model.gcmMetList;
fidin = fopen(model.gcmOutputFile, 'r');

tline1 = fgetl(fidin); % Header row

G = zeros(length(model.mets),length(groups));
counter = 0;

while 1
    tline1 = fgetl(fidin);
    counter = counter + 1;

    if ~ischar(tline1)
        break;
    end

    if ~strcmp('NONE;NONE;', tline1(1:10))
        semiColonIdx = strfind(tline1, ';');
        barIdx = strfind(tline1, '|');
        colonIdx = strfind(tline1, ':');

        thisMetGroup = tline1((semiColonIdx(2) + 1):(colonIdx(1) - 1));
        thisMetGroupCount = str2double(tline1((colonIdx(1) + 1):(barIdx(1) - 1)));

        if any(ismember(groups,thisMetGroup))
            G(ismember(model.mets,metList(counter)),ismember(groups,thisMetGroup)) = thisMetGroupCount;
        else
            error([thisMetGroup, ' not in group list from GCM paper.'])
        end

        for n = 1:(length(barIdx)-1)
            thisMetGroup = tline1((barIdx(n) + 1):(colonIdx(n+1) - 1));
            thisMetGroupCount = str2double(tline1((colonIdx(n+1) + 1):(barIdx(n+1) - 1)));
            if any(ismember(groups,thisMetGroup))
                G(ismember(model.mets,metList(counter)),ismember(groups,thisMetGroup)) = thisMetGroupCount;
            else
                error([thisMetGroup, ' not in group list from GCM paper.'])
            end
        end

    end

end

fclose(fidin);

author: opencobra

repo: cobratoolbox


text:
function psi_values_test ( )

%*****************************************************************************80
%
%% PSI_VALUES_TEST demonstrates the use of PSI_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PSI_VALUES_TEST:\n' );
  fprintf ( 1, '  PSI_VALUES stores values of\n' );
  fprintf ( 1, '  the PSI function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X            PSI(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = psi_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end


text:
%  INTERNAL FUNCTION: this function is a decoy
%
%
%    Other functions named pages2struct
%
%       ts/pages2struct
%

author: jmaih


text:
function pde = DataDivCurlConstant
%%
%           curl u = g    in \Omega,
%           div (eps* u) = f    in \Omega,
%           n \cdot (eps * u) = g_N  on \Gamma_0
% Reference: a constant vector field

alpha = 1;

pde.Eps = @(p) [1+p(:,1)*0 p(:,3)*0 p(:,3)*0; ...
    p(:,3)*0 1+p(:,2)*0 p(:,3)*0; ...
    0*p(:,1) 0*p(:,2) 1+p(:,3)*0];

pde.exactu = @(p) alpha*[ones(size(p,1),1),...
                         ones(size(p,1),1),...
                         ones(size(p,1),1)];
pde.f = @(p) 0*p(:,1);
pde.g = @(p) [0*p(:,1), 0*p(:,2), 0*p(:,3)];

end

repo: ifem


text:
function y=rmDecimate(x,r)
% rmDecimate - decimate the signal by factor (r) along first dimension.
%
%  y=rmDecimate(x,r);
%
% 2007/01 SOD: wrapper for matlab's decimate

% input check
if ~exist('x','var') || isempty(x),
    disp('Need x');
    return;
end

% no decimation if: r is not given, empty or smaller than 2.
if ~exist('r','var') || isempty(r) || r<2,
    y=x;
    return;
end

% get size of x
[s1, s2] = size(x);

% sometimes we input single precision
convert_back_to_single = false;
if isa(x,'single')
    convert_back_to_single = true;
    x=double(x);
end

% filter along first dimension, i should really vectorize this....
y = zeros(ceil(s1./r),s2);
for ii=1:s2
    % The function decimate appears to be gone from later versions of
    % Matlab.  Check what's going on here ...
    y(:,ii) = decimate(x(:,ii),r);
end

if convert_back_to_single,
    y = single(y);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


function Lbw = Lbw(angle_of_attack, angle_of_sideslip)
% Rotation matrix from Wind-axis to Aircraft's body axis
sa = sind(angle_of_attack);
ca = cosd(angle_of_attack);
sb = sind(angle_of_sideslip);
cb = cosd(angle_of_sideslip);
Lbw = [...
    ca*cb -ca*sb -sa
    sb cb 0
    sa*cb -sa*sb ca];
end


text:
% Copy this template configuration file to your VOT workspace.
% Enter the full path to the ECO repository root folder.

ECO_repo_path = 'E:/ForwardBackward/Tracker/FBT_ECO';

tracker_label = 'ECO';
tracker_command = generate_matlab_command('benchmark_tracker_wrapper(''ECO'', ''eco_vot_deep'', true)', {[ECO_repo_path '/VOT_integration/benchmark_wrapper']});
tracker_interpreter = 'matlab';

repo: UDT


text:
% submitWeb Creates files from your code and output for web submission.
%
%   If the submit function does not work for you, use the web-submission mechanism.
%   Call this function to produce a file for the part you wish to submit. Then,
%   submit the file to the class servers using the "Web Submission" button on the
%   Programming Assignments page on the course website.
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function submitWeb(partId)
  if ~exist('partId', 'var') || isempty(partId)
    partId = [];
  end

  submit(partId, 1);
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function point = intersectThreePlanes(plane1, plane2, plane3)
%INTERSECTTHREEPLANES Return intersection point between 3 planes in space.
%
%   LINE = intersectThreePlanes(PLANE1, PLANE2, PLANE3)
%   Returns the point or straight line belonging to three planes.
%   PLANE:  [x0 y0 z0  dx1 dy1 dz1  dx2 dy2 dz2]
%   POINT:  [x0 y0 z0]
%   IF rank of the coefficient matrix r1 = 3 and
%   Rank of the augmented matrix r2 = 3 return point
%   Otherwise returns point with NaN values.
%
%   See also
%   planes3d, intersectPlanes, intersectLinePlane

% ------
% Author: Roozbeh Geraili Mikola
% E-mail: roozbehg@berkeley.edu or roozbehg@live.com
% Created: 2017-09-20
% Copyright 2017-2022

% plane normal
n1 = normalizeVector3d(cross(plane1(:,4:6), plane1(:, 7:9), 2));
n2 = normalizeVector3d(cross(plane2(:,4:6), plane2(:, 7:9), 2));
n3 = normalizeVector3d(cross(plane3(:,4:6), plane3(:, 7:9), 2));

% Uses Hessian form, ie : N.p = d
% I this case, d can be found as : -N.p0, when N is normalized
d1 = dot(n1, plane1(:,1:3), 2);
d2 = dot(n2, plane2(:,1:3), 2);
d3 = dot(n3, plane3(:,1:3), 2);

% create coefficient and augmented matrices
A = [n1;n2;n3];
D = [d1;d2;d3];
AD = [n1,d1;n2,d2;n3,d3];

% calculate rank of the coefficient and augmented matrices
r1 = rank(A);
r2 = rank(AD);

% if rank of the coefficient matrix r1 = 3 and
% rank of the augmented matrix r2 = 3 return point
% and if r1 = 2 and r2 = 2 return line,
% otherwise returns point with NaN values.
if r1 == 3 && r2 == 3
    % Intersecting at a point
    point = (A\D)';
else
    point = [NaN NaN NaN];
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function x = esvm_hog(I, sbin)
%Return the current feature function, same as esvm_features

if nargin == 0
  x = 31;
  return
end

x = features_pedro(I,sbin);
%x = features_raw(I,sbin);

author: quantombone

repo: exemplarsvm


text:
function [X,Y,Z] = cov3elli(x,P,ns,NP)

% COV3ELLI  3D ellipsoid from Gaussian mean and covariance.
%   [X,Y,Z] = COV3ELLI(x0,P,ns,NP) gives X, Y and Z coordinates of the
%   points corresponding to the 2 biggest semi-diametres of the ellipsoid
%   defined by the covariances matrix P and centered at x0:
%
%        (x-x0)'*(P^-1)*(x-x0) = ns^2.
%
%   The ellipsoid can be plotted in a 3D graphic by just creating a line
%   with line(X,Y,Z).
%
%   See also COV2ELLI, IDP3ELLI, LINE.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

persistent spheroid

% Basic shape: 2 circles at 90 degrees, with the pole at the major axis X.
if isempty(spheroid)
    alpha = 2*pi/NP*(0:NP);
    spheroid = [...
        cos(alpha)    cos(alpha)
        sin(alpha)    zeros(1,NP+1)
        zeros(1,NP+1) sin(alpha)];
end

% Rotation R and semi-diameters d, obtained from P by SVD
[R,D] = svd(P);  % rotation matrix R
d     = sqrt(D); % semi-axes

% spheroid -> aligned ellipsoid -> rotated ellipsoid -> ns-ellipsoid
ellip = ns*R*d*spheroid;

% Choleski is not good here: ellipsoid cross not aligned w/ axis
% C = chol(P)';
% ellip = ns*C*spheroid;

% output ready for plotting (X, Y and Z line vectors), offset by x
X = x(1)+ellip(1,:);
Y = x(2)+ellip(2,:);
Z = x(3)+ellip(3,:);

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [mi3] = oz2mi3(oz)
% Convert volume from US liquid ounces to cubic miles.
% Chad Greene 2012
mi3 = oz*7.0950670584e-15;


text:
function jed = epoch_to_jed_zoroastrian ( )

%*****************************************************************************80
%
%% EPOCH_TO_JED_ZOROASTRIAN: epoch of the Zoroastrian calendar as a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real JED, the Julian Ephemeris Date of the epoch.
%
  jed = 1862836.5;

  return
end


text:
%% Tracking of rotating point using Kalman filter
%
% Tracking of rotating point.
% Rotation speed is constant.
% Both state and measurements vectors are 1D (a point angle),
% Measurement is the real point angle + gaussian noise.
% The real and the estimated points are connected with yellow line segment,
% the real and the measured points are connected with red line segment.
% (if Kalman filter works correctly,
% the yellow segment should be shorter than the red one).
%
% Pressing any key will reset the tracking with a different speed.
% Close the window to stop the program.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/kalman.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/kalman.py>
%

%%
% set up display window
img = zeros(500, 500, 3, 'uint8');
hFig = figure('KeyPressFcn',@(o,e)setappdata(o, 'flag',true), ...
    'Menubar','none', 'Name','Kalman Filter demo');
setappdata(hFig, 'flag',false);
hImg = imshow(img);

%%
% helper anonymous functions
calcPoint = @(center,R,angle) center + [cos(angle), -sin(angle)]*R;
drawCross = @(img,center,clr,d) cv.line(...
    cv.line(img, center-d, center+d, 'Color',clr, ...
        'Thickness',1, 'LineType','AA'), ...
    center+[d -d], center+[-d d], 'Color',clr, ...
        'Thickness',1, 'LineType','AA');

%%
% create and initialize Kalman filter
KF = cv.KalmanFilter(2, 1);
state = zeros(2,1);  % [phi; delta_phi]
processNoise = zeros(2,1);
measurement = zeros(1,1);

%%
% keep repeating until figure is closed
while ishghandle(hFig)
    setappdata(hFig, 'flag',false);

    % initialize KF
    state = randn(size(state))*0.1;
    KF.transitionMatrix = [1 1; 0 1];
    KF.measurementMatrix = eye(size(KF.measurementMatrix));
    KF.processNoiseCov = eye(size(KF.processNoiseCov))*1e-5;
    KF.measurementNoiseCov = eye(size(KF.measurementNoiseCov))*1e-1;
    KF.errorCovPost = eye(size(KF.errorCovPost));
    KF.statePost = randn(size(KF.statePost))*0.1;

    % main loop
    while ishghandle(hFig)
        center = [size(img,2) size(img,1)]/2;
        R = size(img,2)/3;
        stateAngle = state(1);
        statePt = calcPoint(center, R, stateAngle);

        prediction = KF.predict();
        predictAngle = prediction(1);
        predictPt = calcPoint(center, R, predictAngle);

        measurement = randn(size(measurement))*KF.measurementNoiseCov(1);

        % generate measurement
        measurement = measurement + KF.measurementMatrix*state;

        measAngle = measurement(1);
        measPt = calcPoint(center, R, measAngle);

        % plot points
        img(:) = 0;
        img = drawCross(img, statePt, [255 255 255], 3);
        img = drawCross(img, measPt, [255 0 0 0], 3);
        img = drawCross(img, predictPt, [0 255 0], 3);
        img = cv.line(img, statePt, measPt, 'Color',[255 0 0], ...
            'Thickness',3, 'LineType','AA');
        img = cv.line(img, predictPt, measPt, 'Color',[255 255 0], ...
            'Thickness',3, 'LineType','AA');

        if rand > 0.75
            KF.correct(measurement);
        end

        processNoise = randn(size(processNoise))*sqrt(KF.processNoiseCov(1,1));
        state = KF.transitionMatrix*state + processNoise;

        % update display
        set(hImg, 'CData',img);

        % break of inner loop on any key press
        flag = getappdata(hFig, 'flag');
        if isempty(flag)||flag, break; end
        pause(0.1)
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function r = uminus(c)

%tstoolbox/@core/uminus
%   Syntax:
%     * r = uminus(c)
%
%   Input Arguments:
%     * c - core object
%
%   negate time series
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

r = core(-data(c));


text:
% Script file to find parameters for simplified valve actuator model
% Copyright 2010 MathWorks, Inc.

% This script file invokes optimization process to find out a suitable
% approximation to the characteristic of the 2-Position Valve Actuator with
% a simple, Simulink-based model of the actuator

% Required characteristics
t_on = 0.2;         % [s]
stroke = 0.005;     % [m]
gain = stroke/12;   % 12 is a command signal

% Two actuator parameters, time constant and delay time are determined
% during optimization

% Set initial values for variable parameters
x0 = [0.1 0.2/3];

% Perform optimization
[x,fval,exitflag,output] = ...
    fminsearch(@obj_find_simple_valve_driver_param,x0, ...
    optimset('Tolx',1e-6,'Display','iter'));


text:
function [mm] = m2mm(m)
% Convert length from meters to millimeters.  Trivial, yes, but using this
% function may help make your code more clear.  It'll save you a year from
% now when you revisit the code you're writing and you say, "why did I
% multiply this number by a thousand?"
mm = m*1000;


text:
function res = getEigenVals(a)
res = a.eigenVals;

author: thomaskuestner


text: addpath(genpath(pwd));

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function [Eft1, Eft2, Covf, lpyt] = pred_coxph(gp, x, y, xt, varargin)
% PRED_COXPH Wrapper for returning useful values for coxph likelihood
%
%  Description
%     [EF1, EF2, COVF, LPYT] = PRED_COXPH(GP,X,Y,XT, OPTIONS)
%     Returns predictions for the two latent processes EF1 and EF2,
%     covariance matrix COVF and logarithms of predictive densities at XT.

% Copyright (c) 2012-2013 Ville Tolvanen

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

ip=inputParser;
ip.FunctionName = 'PRED_COXPH';
ip.addRequired('gp',@isstruct);
ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('xt',  @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('yt', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('zt', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('predcf', [], @(x) isempty(x) || ...
                 isvector(x) && isreal(x) && all(isfinite(x)&x>0))
ip.addParamValue('tstind', [], @(x) isempty(x) || iscell(x) ||...
                 (isvector(x) && isreal(x) && all(isfinite(x)&x>0)))
ip.parse(gp, x, y, xt, varargin{:});

if ~strcmp(gp.lik.type, 'Coxph')
  error('Likelihood not Coxph')
end
if nargout > 3
  if numel(gp.jitterSigma2)==1
    [Ef, Covf, lpyt] = gp_pred(gp, x, y, xt, varargin{:});
  else
    [Ef, Covf, lpyt] = gpmc_preds(gp, x, y, xt, varargin{:});
  end
else
  if numel(gp.jitterSigma2)==1
    [Ef, Covf] = gp_pred(gp, x, y, xt, varargin{:});
  else
    [Ef, Covf] = gpmc_preds(gp, x, y, xt, varargin{:});
  end
end
ntime=size(gp.lik.xtime,1);
if isfield(gp.lik, 'stratificationVariables')
  ind_str=gp.lik.stratificationVariables;
  nf1=ntime.*unique([x(:,ind_str); xt(:,ind_str)], 'rows');
else
  nf1=ntime;
end

Eft1 = Ef(1:nf1,:); Ef(1:nf1,:) = []; Eft2 = Ef;

end

author: fieldtrip

repo: fieldtrip


text:
function types = guitypes()
% GUITYPES returns a cell array of MATLAB/Octave GUI object types
%
% Syntax
%   types = guitypes()
%
% These types are ignored by matlab2tikz and figure2dot.
%
% See also: matlab2tikz, figure2dot

types = {'uitoolbar', 'uimenu', 'uicontextmenu', 'uitoggletool',...
         'uitogglesplittool', 'uipushtool', 'hgjavacomponent'};
end

author: thomas-koehler

repo: SupER


text:
function p = sup(p);
%SUP          Polynomial of upper bounds for (interval) polynomial (same as p.sup)
%
%   r = sup(p)
%

% written  08/28/00     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  p.c = sup(p.c);
  p = normalize(p);

repo: OpenMAS


text:
%{

This code is to render a Mesh given a 3x4 camera matrix with an image resolution widthxheight. The rendering result is an ID map for facets, edges and vertices. This can usually used for occlusion testing in texture mapping a model from an image, such as the texture mapping in the following two papers.

--Jianxiong Xiao http://mit.edu/jxiao/

Citation:

[1] J. Xiao, T. Fang, P. Zhao, M. Lhuillier, and L. Quan
Image-based Street-side City Modeling
ACM Transaction on Graphics (TOG), Volume 28, Number 5
Proceedings of ACM SIGGRAPH Asia 2009

[2] J. Xiao, T. Fang, P. Tan, P. Zhao, E. Ofek, and L. Quan
Image-based Facade Modeling
ACM Transaction on Graphics (TOG), Volume 27, Number 5
Proceedings of ACM SIGGRAPH Asia 2008

%}

clear
clc

syms P00 P01 P02 P03 P10 P11 P12 P13 P20 P21 P22 P23
syms m_far m_near m_width m_height scale

projection = [P00 P01 P02 P03; P10 P11 P12 P13; P20 P21 P22 P23; 0 0 0 1];

% handle the near and far clip plane in OpenGL
protr = [1.0, 0.0, 0.0, 0.0;
    0.0, 1.0, 0.0, 0.0;
    0.0, 0.0, m_far / (m_far - m_near), - m_near * m_far /(m_far - m_near);
    0.0, 0.0, 1.0, 0.0];

% handle half pixel inconsistency
offset = [1.0, 0.0, 0.5, 0.0;
    0.0, 1.0, 0.5, 0.0;
    0.0, 0.0, 1.0, 0.0;
    0.0, 0.0, 0.0, 1.0];

% undo image aspect ratio and size
m0 = [m_width / 2, 0.0,         0.0, 0 + m_width / 2.0;
    0.0,        m_height / 2, 0.0, 0 + m_height / 2.0;
    0.0,        0.0,         0.5, 0.5;
    0.0,        0.0,         0.0, 1.0];

% handle scaling
m1=[1 / scale,   0.0,         0.0,                0.0;
    0.0,        1 / scale,    0.0,                0.0;
    0.0,        0.0,         1.0,                0.0;
    0.0,        0.0,         0.0,                1.0];

% handle upside down in vertical direction in image
m2=[1,           0.0,         0.0,                0.0;
    0.0,        -1.0,        m_height,            0.0;
    0.0,        0.0,         1.0,                0.0;
    0.0,        0.0,         0.0,                1.0];

M = (m0 \ protr * m2 * m1 * offset * projection);

M =

[     P20*(1/(m_width*scale) - 1) + (2*P00)/(m_width*scale),     P21*(1/(m_width*scale) - 1) + (2*P01)/(m_width*scale),     P22*(1/(m_width*scale) - 1) + (2*P02)/(m_width*scale),                       P23*(1/(m_width*scale) - 1) + (2*P03)/(m_width*scale)]
[ - P20*(1/(m_height*scale) - 1) - (2*P10)/(m_height*scale), - P21*(1/(m_height*scale) - 1) - (2*P11)/(m_height*scale), - P22*(1/(m_height*scale) - 1) - (2*P12)/(m_height*scale),                   - P23*(1/(m_height*scale) - 1) - (2*P13)/(m_height*scale)]
[                   (P20*(m_far + m_near))/(m_far - m_near),                   (P21*(m_far + m_near))/(m_far - m_near),                   (P22*(m_far + m_near))/(m_far - m_near), (P23*(m_far + m_near))/(m_far - m_near) - (2*m_far*m_near)/(m_far - m_near)]
[                                                       P20,                                                       P21,                                                       P22,                                                                         P23]

inv_width_scale  = 1/(m_width*scale);
inv_height_scale = 1/(m_height*scale);
inv_width_scale_1 =inv_width_scale - 1;
inv_height_scale_1_s = -(inv_height_scale - 1);
inv_width_scale_2 = inv_width_scale*2;
inv_height_scale_2_s = -inv_height_scale*2;
m_far_a_m_near = m_far + m_near;
m_far_s_m_near = m_far - m_near;
m_far_d_m_near = m_far_a_m_near/m_far_s_m_near;

M =
[     P20*(inv_width_scale - 1) + (2*P00)*inv_width_scale,   P21*(inv_width_scale - 1) + (2*P01)*inv_width_scale  , P22*(inv_width_scale - 1) + (2*P02)*inv_width_scale    ,                       P23*(inv_width_scale - 1) + (2*P03)*inv_width_scale]
[ - P20*(inv_height_scale - 1) - (2*P10)*inv_height_scale, - P21*(inv_height_scale - 1) - (2*P11)*inv_height_scale, - P22*(inv_height_scale - 1) - (2*P12)*inv_height_scale,                   - P23*(inv_height_scale - 1) - (2*P13)*inv_height_scale]
[                   (P20*m_far_a_m_near)/m_far_s_m_near,                   (P21*m_far_a_m_near)/m_far_s_m_near,                   (P22*m_far_a_m_near)/m_far_s_m_near, (P23*m_far_a_m_near)/m_far_s_m_near - (2*m_far*m_near)/m_far_s_m_near]
[                                                       P20,                                                       P21,                                                       P22,                                                                         P23]

% final OpenGL matrix = M'
% column 1
P20*inv_width_scale_1 + P00*inv_width_scale_2,
P21*inv_width_scale_1 + P01*inv_width_scale_2,
P22*inv_width_scale_1 + P02*inv_width_scale_2,
P23*inv_width_scale_1 + P03*inv_width_scale_2,

% column 2
P20*inv_height_scale_1_s + P10*inv_height_scale_2_s,
P21*inv_height_scale_1_s + P11*inv_height_scale_2_s,
P22*inv_height_scale_1_s + P12*inv_height_scale_2_s,
P23*inv_height_scale_1_s + P13*inv_height_scale_2_s,

% column 3
P20*m_far_d_m_near,
P21*m_far_d_m_near,
P22*m_far_d_m_near,
P23*m_far_d_m_near - (2*m_far*m_near)/m_far_s_m_near,
% column 4
 P20
 P21
 P22
 P23

% final matrix to opengl

P20*inv_width_scale_1 + P00*inv_width_scale_2,  P20*inv_height_scale_1_s + P10*inv_height_scale_2_s,    P20*m_far_d_m_near,     P20,
P21*inv_width_scale_1 + P01*inv_width_scale_2,  P21*inv_height_scale_1_s + P11*inv_height_scale_2_s,    P21*m_far_d_m_near,     P21,
P22*inv_width_scale_1 + P02*inv_width_scale_2,  P22*inv_height_scale_1_s + P12*inv_height_scale_2_s,    P22*m_far_d_m_near,     P22,
P23*inv_width_scale_1 + P03*inv_width_scale_2,  P23*inv_height_scale_1_s + P13*inv_height_scale_2_s,    P23*m_far_d_m_near - (2*m_far*m_near)/m_far_s_m_near,   P23

inv_width_scale  = 1/(m_width*scale);
inv_height_scale = 1/(m_height*scale);
inv_width_scale_1 =inv_width_scale - 1;
inv_height_scale_1_s = -(inv_height_scale - 1);
inv_width_scale_2 = inv_width_scale*2;
inv_height_scale_2_s = -inv_height_scale*2;
m_far_a_m_near = m_far + m_near;
m_far_s_m_near = m_far - m_near;
m_far_d_m_near = m_far_a_m_near/m_far_s_m_near;

author: jianxiongxiao

repo: ProfXkit


text:
function [ labels, valid_ids, vid_ids, vid_names  ] = extract_SEMAINE_labels( SEMAINE_dir, recs, aus )
%EXTRACT_SEMAINE_LABELS Summary of this function goes here
%   Detailed explanation goes here

    % Get the right eaf file

    aus_SEMAINE = [2 12 17 25 28 45];

    inds_to_use = [];

    for i=1:numel(aus)

        inds_to_use = cat(1, inds_to_use, find(aus_SEMAINE == aus(i)));

    end

    labels = cell(numel(recs), 1);
    valid_ids = cell(numel(recs), 1);
    vid_names = cell(numel(recs), 1);
    vid_ids = zeros(numel(recs), 2);

    for i=1:numel(recs)

        file = dir([SEMAINE_dir, '/', recs{i}, '/*.eaf']);

        vid_ids(i,:) = dlmread([SEMAINE_dir, '/', recs{i}, '.txt'], ' ');

        vid_names_c = dir([SEMAINE_dir, '/', recs{i}, '/*.avi']);
        [~, vid_names{i},~] = fileparts(vid_names_c.name);

        xml_file = [SEMAINE_dir, recs{i}, '\' file.name];
        [root_xml, name_xml, ~] = fileparts(xml_file);

        activations = ParseSEMAINEAnnotations([SEMAINE_dir, recs{i}, '/' file.name]);

        if(size(activations,1) < vid_ids(i,2))
            vid_ids(i,2) = size(activations,1);
            if(vid_ids(i,2) > 2999)
                vid_ids(i,1) = vid_ids(i,2) - 2999;
            end
        end

        labels{i} = activations(vid_ids(i,1)+1:vid_ids(i,2), 1 + inds_to_use);

        % all indices in SEMAINE are valid
        valid_ids{i} = ones(size(labels{i},1),1);

    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function MakefMasksForAllFish(hfig)

% batch input data setup:
isFullData = 1;
data_masterdir = GetCurrentDataDir();

% M_fishset = GetFishStimset();
% M_stimrange = GetStimRange();

range_fish =  11:18; % range_fish = GetFishRange();

%% custom params here:
thres_minsize = 10000;
thres_maxsize = 6*10^5;

%%
for i = 1:length(range_fish),
    i_fish = range_fish(i);
    disp(i_fish);

        option
    [cIX,gIX,M_xyz_norm] = GetDefaultClustersFromLoad(hfig,i_fish,option);

    % Load fish
    LoadFullFish(hfig,i_fish,-1);
    absIX = getappdata(hfig,'absIX');
    ClusterIDs = GetClusterIDs();
    [cIX,gIX] = LoadCluster_Direct(i_fish,ClusterIDs(1),ClusterIDs(2),absIX);
    anat_stack_norm = getappdata(hfig,'anat_stack_norm');
    CellXYZ_norm = getappdata(hfig,'CellXYZ_norm');

    % make fMask for each auto-cluster
    fMASKs = [];
    MaskDatabase = sparse(0,0);
    mask_clusID = [];
    gIX = SqueezeGroupIX(gIX);
    for i_clus = 1:length(unique(gIX)),
        IX = find(gIX==i_clus);
        newMask = MakeFunctionalMask(anat_stack_norm,CellXYZ_norm,absIX,cIX(IX),'isBatch');
        sz = length(find(newMask));
        if sz>thres_minsize && sz<thres_maxsize,
            MaskDatabase(:, end+1) = newMask;
            mask_clusID = [mask_clusID;i_clus];
        end
    end
    fMASKs.MaskDatabase = MaskDatabase;
    fMASKs.mask_clusID = mask_clusID;
    fMASKs.i_fish = i_fish;
    fMASKs.note = 'Auto_defS_woA_Master0.5';
    % save:
    name = ['fMASKs_',num2str(i_fish),'.mat'];
    save(fullfile(data_masterdir,'fMASKs',name),'fMASKs');
end

author: xiuyechen

repo: FishExplorer

path: github-repos/MATLAB/xiuyechen-FishExplorer/FishExplorer-c61392cf0835480d64fc03c15f1992935fdc7106/unused analyses/conserved cluster screens/obsl mask screen/MakefMasksForAllFish.m


text:
function chrpak_test046 ( )

%*****************************************************************************80
%
%% TEST046 tests HEX_TO_I4 and I4_TO_HEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 June 2007
%
%  Author:
%
%    John Burkardt
%
  test_num = 3;

  i4_test = [ 21, -32, 1776 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST046\n' );
  fprintf ( 1, '  HEX_TO_I4, hexadecimal->integer.\n' );
  fprintf ( 1, '  I4_TO_HEX, integer->hexadecimal\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I  I4_TO_HEX(I)  HEX_TO_I4(I4_TO_HEX(I))\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    i4 = i4_test(test);
    hex = i4_to_hex ( i4 );
    j4 = hex_to_i4 ( hex );

    fprintf ( 1, '  %8d  %8s  %8d\n', i4, hex, j4 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function linplus_test56 ( )

%*****************************************************************************80
%
%% TEST56 tests R8SM_ML.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  m = 7;
  n = m;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST56' );
  fprintf ( 1, '  R8SM_ML computes A*x or A''*X\n' );
  fprintf ( 1, '    where A is a Sherman Morrison matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M =    %d\n', m );
  fprintf ( 1, '  Matrix columns N = %d\n', n );

  for job = 0 : 1
%
%  Set the matrix.
%
    [ a, u, v, seed ] = r8sm_random ( m, n, seed );

    r8sm_print ( m, n, a, u, v, '  The Sherman Morrison matrix:' );
%
%  Set the desired solution.
%
    x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
    if ( job == 0 )
      b = r8sm_mxv ( m, n, a, u, v, x );
    else
      b = r8sm_vxm ( m, n, a, u, v, x );
    end
%
%  Factor the matrix.
%
    [ a_lu, pivot, info ] = r8ge_fa ( n, a );

    if ( info ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Fatal error!\n' );
      fprintf ( 1, '  R8GE_FA declares the matrix is singular!\n' );
      fprintf ( 1, '  The value of INFO is %d\n', info );
      return
    end
%
%  Now multiply factored matrix times solution to get right hand side again.
%
    b2 = r8sm_ml ( n, a_lu, u, v, pivot, x, job );

    if ( job == 0 )
      r8vec2_print_some ( n, b, b2, 10, '  A*x and PLU*x' );
    else
      r8vec2_print_some ( n, b, b2, 10, '  A''*x and (PLU)''*x' );
    end

  end

  return
end


text: StartSIFT;

repo: REST


text:
function g = modelOutputGradX(model, X)

% MODELOUTPUTGRADX Compute derivatives with respect to model inputs of model outputs.
% FORMAT
% DESC gives the gradients of the outputs from the model with
% respect to the inputs.
% ARG model : the model structure for which gradients are computed.
% ARG X : input locations where gradients are to be computed.
% RETURN g : gradients of the model output with respect to the
% model parameters for the given input locations.
%
% SEEALSO : modelCreate, modelOutputGrad, modelLogLikelihood, modelLogLikeGradients
%
% COPYRIGHT : Neil D. Lawrence, 2006

% MLTOOLS

fhandle = str2func([model.type 'OutputGradX']);
g = fhandle(model, X);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% ACOT   Inverse cotangent, result in radian.
%    ACOT(X) is the inverse cotangent of the elements of X.
%
%    Class support for input X:
%       float: double, single
%
%    See also COT, ACOTD.
%
%    Reference page in Doc Center
%       doc acot
%
%    Other functions named acot
%
%       codistributed/acot    gpuArray/acot    sym/acot    ts/acot
%


text:
% make_avi_movie_example1.m
%
%  Discussion:
%
%    This is a simple example program to create an
%    Audio Video Interleaved (AVI) movie that can be played
%    independently of Matlab, using, for example, the XINE player
%    for Linux. We capture NUMFRAMES frames for the 1-D sine
%    wave with phase changing from 0 to 2 pi.
%
%    Before running this program close all figures and don't interfere
%    with the generated figures until the recording process is done.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2004
%
%  Author:
%
%    Marcus Garvie
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MAKE_AVI_MOVIE_EXAMPLE1\n' );
  fprintf ( 1, '  Create an AVI animation of simple graphics,\n' );
  fprintf ( 1, '  generating one frame at a time.\n' )
%
%  Set the total number of frames that we will generate.
%
  numframes = 100;
  hp = 2 * pi / numframes;
%
%  Control how fast final movie will play; total time
%  of movie in seconds will be about (numframes)/(num_frames_per_second)
%
  num_frames_per_second = 10;
  dur = numframes / num_frames_per_second;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  This movie will contain %d frames.\n', numframes );
  fprintf ( 1, '  The number of frames per second will be %f\n', num_frames_per_second );
  fprintf ( 1, '  so the movie should take %f seconds to play.\n', dur );
%
%  The AVIFILE function creates a new Audio Video Interleaved (AVI) file.
%  We specify its name, and we also set the value of "FPS", which is the
%  rate at which we will display frames (individual snapshots) per second.
%
  aviobj = avifile ( 'sinwave.avi', 'fps', num_frames_per_second );
%
%  Control resolution of wave
%
  num_nodes = 100;
  h = 2 * pi / num_nodes;
  i = 1 : num_nodes + 1;
%
%  Set the X values at which the function will be evaluated.
%
  x = ( i - 1 ) * h;
%
%  Now we make plots of sine waves.  Each new plot is a frame which
%  can be added to the AVI movie using the ADDFRAME function.
%
  for j =  1 : numframes+1
    y = sin ( x - ( j - 1 ) * hp );
    plot ( x, y )
    axis tight
    frame = getframe ( gca );
    aviobj = addframe ( aviobj, frame );
  end
%
%  Tell MATLAB we have completed the movie.
%
  aviobj = close ( aviobj );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'MAKE_AVI_MOVIE_EXAMPLE1\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '  The movie file "sinwave.avi" has been created.\n' );


text:
function [] = visGaussianClassLearning(RBM);
%-----------------------------------------
%  [] =  visGaussianClassLearning(RBM);
%-----------------------------------------
% DES

nVis = floor(sqrt(size(RBM.W,1))).^2;

subplot(331);
scatter(RBM.auxVars.batchX(:,1),RBM.auxVars.batchX(:,2),'r.');

title('Batch Data');
xlim([-10, 10]); ylim([-20, 0]);
%  xlim([-5, 5]); ylim([-5, 5]);

subplot(332);
scatter(RBM.pVis(:,1),RBM.pVis(:,2),'b.');
title('Fantasies');
%  xlim([-5, 5]); ylim([-5, 5]);
xlim([-10, 10]); ylim([-20, 0]);

subplot(333);
visWeights(flipud(RBM.dW(1:nVis,:)));
title ('Weight Gradients');

subplot(334);
bar(RBM.b);
title('Visible Bias');

%
%  subplot(335)
%  bar((RBM.sigma2));
%  title('Visible Variance');

subplot(336);
visWeights(flipud(RBM.W(1:nVis,:)));
title('Connection Weights');

subplot(337);
plot(RBM.auxVars.error);
title('Reconstruction errors');

subplot(339);
semilogy(RBM.auxVars.lRate);
title('Learning Rate');

drawnow

author: dustinstansbury

repo: medal


text:
%This directory contains functions for statistical analysis.
%
%STAT_NORMAL_CDF - evaluates the normal cumulative distribution function.
%STAT_RMANOVA - repeated-measures ANOVA
%STAT_CALCTCRIT - calculates the critical t-value given significance threshold and degrees of freedom
%STAT_PERCENTILES - returns the percentiles for data samples.

author: bbci


text:
function Y=hooklaw(X, td, E, ni);

%Y=hooklaw(X, td, E, ni)
%
%This function calculates the strain or stress according to Hook's Law.
%Transform direction determine the input variable td.
%The proper values of td variable:
%td='stress_strain' - if input variable X is a stress tensor
%td='strain_stress' - if input variable X is a strain tensor
%
%But, for the purpose of computations we use the vector representation
%of a tensor X(ij), so that:
%X - is a vector representation of a tensor X(ij)
%    with the components in the following order:
%
%X=[Xxx Xyy Xzz Xxy Xxz Xzx Xyx Xyz Xxz]
%or
%X=[Xxx Xyy Xzz Xxy Xxz Xyz]
%or
%X=[Xxx Xyy Xzz]
%
%Dimension of X can be [m n]
%where: m - represents the number of tensors to calculate
%       n = 9, 6 or 3 (see above)
%
%E  - Young's modulus
%ni - Poisson's coefficient
%
%   Copyright (c) 2001 by Aleksander Karolczuk.
%   $Revision: 1.2 $  $Date: 2004/02/13

error(nargchk(4,4,nargin))

[m n]=size(X);

if n==3
   I=[1 1 1];
elseif n==6
   I=[1 1 1 0 0 0];
elseif n==9
   I=[1 1 1 0 0 0 0 0 0];
else
   error('Improper matrix dimension')
end

if lower(td)=='stress_strain'
   Y=((1+ni)/E)*(X-(ni/(1+ni))*(X*I')*I);
   elseif lower(td)=='strain_stress'
      Y=(E/(1+ni))*(X+(ni/(1-2*ni))*(X*I')*I);
   else
      error('Improper name of transform direction')
end

author: Sable

repo: mcbench-benchmarks


text:
function as = make3daxes(h)
% make3daxes creates a dextral set of coordinate axes for animation.
%
% as = make3daxes returns an array (as) continaing handles to three surface
% objects created in figure h representing the unit vectors of a dextral
% reference frame.
%
% Note: This function is intended to be called by animEuler.m
%
% Examples:
%     as = make3daxes(1);

% Written by Dmitry Savransky 29 April 2009

%set figure and remove all children
figure(h);
clf(h,'reset');
hold off

%create an axis
[xrod,yrod,zrod] = cylinder(0.05,12);

%create and place the axes
a3 = surface(xrod,yrod,zrod,'FaceColor','b','AlphaData',0.5);
a2 = surface(xrod,yrod,zrod,'FaceColor','g','AlphaData',0.5);
a1 = surface(xrod,yrod,zrod,'FaceColor','r','AlphaData',0.5);
rotate(a2,[1 0 0],-90,[0,0,0])
rotate(a1,[0 1 0],90,[0,0,0])

%clean up display
axis equal;
axis([-1 1 -1 1 -1 1])
view(100,15)
grid on
xlabel('x');
ylabel('y');
zlabel('z');

%return handles
as = [a1 a2 a3];


text:
function disp (F)
%DISP displays the factorization F
%
% Example
%   F = factorize (A)
%   disp (F)
%
% See also factorize.

% Copyright 2009, Timothy A. Davis, University of Florida

fprintf ('  A:\n') ;
disp (F.A) ;

if (~isempty (F.L))
    fprintf ('  L:\n') ;
    disp (F.L) ;
end

if (~isempty (F.U))
    fprintf ('  U:\n') ;
    disp (F.U) ;
end

if (~isempty (F.Q))
    fprintf ('  Q:\n') ;
    disp (F.Q) ;
end

if (~isempty (F.R))
    fprintf ('  R:\n') ;
    disp (F.R) ;
end

if (~isempty (F.p))
    fprintf ('  p:\n') ;
    disp (F.p) ;
end

if (~isempty (F.q))
    fprintf ('  q:\n') ;
    disp (F.q) ;
end

fprintf ('  is_inverse: %d kind: %d\n', F.is_inverse, F.kind) ;

% print the kind of factorization that F contains

switch F.kind

    case 1

        fprintf ('  Q-less economy sparse QR factorization: ') ;
        fprintf ('(A*q)''*A*q = R''*R\n');

    case 2

        fprintf ('  dense economy QR factorization: A = Q*R\n') ;

    case 3

        fprintf ('  Q-less economy sparse QR factorization: ') ;
        fprintf ('(p*A)*(p*A)'' = R''*R\n') ;

    case 4

        fprintf ('  dense economy QR factorization: A'' = Q*R\n') ;

    case 5

        fprintf ('  sparse Cholesky factorization: ')  ;
        fprintf ('q*A*q'' = L*L''\n');

    case 6

        fprintf ('  dense Cholesky factorization: A = L*L''\n') ;

    case 7

        fprintf ('  sparse LU factorization: p*A*q = L*U\n') ;

    case 8

        fprintf ('  dense LU factorization: p*A = L*U\n') ;

end


text:
function plotpoint(i,x1,y1,Color)

%x1=[source dest]--in x axis
%y1=[source dest]--in y axis

   %Color='r';
  % N=1000;

    if abs(x1(1)-x1(2))>=abs(y1(1)-y1(2))

     if x1(1)<=x1(2)
        %x=x1(1):x1(2);
        x=x1(1)+i;
        p=1;
     else
        %x=x1(1):-1:x1(2);
        x=x1(1)-i;
        p=0;
     end
             if ((p==1 && x<=x1(2)) || (p==0 && x>=x1(2)))

                %%
                Radius=i;%sqrt((abs(i))^2 + (abs(i))^2);
                plotcircle(x1(1),y1(1),Radius,Color)
                %%

                y=(((x-x1(1))./((x1(2)-x1(1)))).*((y1(2)-y1(1))))+y1(1);

                plot(x,y,'.','LineWidth',1,...
                    'MarkerEdgeColor','k',...
                    'MarkerFaceColor','y',...
                    'MarkerSize',8);
                hold on

             %   plot(x1,y1,'^r')

             end
    else

      if y1(1)<=y1(2)
        %x=x1(1):x1(2);
        y=y1(1)+i;
        p=1;
      else
        %x=x1(1):-1:x1(2);
        y=y1(1)-i;
        p=0;
      end

          if ((p==1 && y<=y1(2)) || (p==0 && y>=y1(2)))
                x=(((y-y1(1))./((y1(2)-y1(1)))).*((x1(2)-x1(1))))+x1(1);

                %%
                Radius=i;%sqrt((abs(i))^2 + (abs(i))^2);
                plotcircle(x1(1),y1(1),Radius,Color)

                %%
                plot(x,y,'.','LineWidth',1,...
                    'MarkerEdgeColor','k',...
                    'MarkerFaceColor','y',...
                    'MarkerSize',8);
                hold on

              %  plot(x1,y1,'^r')

          end
    end

end


text:
function ht= axis_yTickLabel(label_list, varargin)
%AXIS_YTICKLABEL - Add y-ticklabels to the current axis
%
%Synopsis:
% H= axis_yTickLabel(LABEL_LIST, <OPT>)
%
%Input:
% LABEL_LIST: string (labels separated by '|') or cell array of strings
% OPT: struct or property/value list of optional properties:
%  .HPos  - horizontal position
%  .YTick - position of yticks. If not provided, get(gca,'YTick') is used.
%  .Color - font color in RGB format. Maybe an [nLabels 3] matrix of color
%           codes
%  .HorizontalAlignment, .VerticalAlignment
%
%Output:
% H: handle to the text object(s)
%
%Note:
% The position of the text object is defined *within* the axis. So you should
% set XLimMode and YLimMode to 'manual' before calling AXIS_TITLE.

% Author(s): Benjamin Blankertz

props= {'YTick'                []       'DOUBLE[-]'
        'HPos'                 -0.03    'DOUBLE[1]'
        'Color'                [0 0 0]  'DOUBLE[- 3]'
        'HorizontalAlignment' 'right'   'CHAR(left center right)'
        'VerticalAlignment'   'middle'  'CHAR(top cap middle baseline bottom)'};

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);

if ~iscell(label_list),
  if ismember('|',label_list,'legacy'),
    label_str= label_list;
    idx= find(label_str=='|');
    idx= [0, idx, length(label_str)+1];
    nLab= length(idx)-1;
    label_list= cell(1,nLab);
    for ii= 1:nLab,
      label_list{ii}= label_str(idx(ii)+1:idx(ii+1)-1);
    end
  else
    label_list= {label_list};
  end
end
nLab= length(label_list);
if isdefault.YTick,
  opt.YTick= get(gca, 'YTick');
end
if length(opt.YTick)~=nLab,
  error('number of labels must match number of xticks');
end
if nLab>1 && size(opt.Color,1)==1,
  opt.Color= repmat(opt.Color, [nLab 1]);
end

opt_fn= fieldnames(opt);
ifp= find(ismember(opt_fn, {'HorizontalAlignment','VerticalAlignment'},'legacy'));
ifp= cat(1, ifp, strmatch('Font', opt_fn));
font_opt= struct_copyFields(opt, opt_fn(ifp));
font_pl= opt_structToProplist(font_opt);

XLim= get(gca, 'XLim');
xx= XLim(1) + opt.HPos*diff(XLim);
yy= opt.YTick;

for tt= 1:nLab,
  ht(tt)= text(xx, yy(tt), label_list{tt});
  set(ht(tt), font_pl{:}, 'Color',opt.Color(tt,:));
end


text:
function [P Fcost Emi Pl]=emield(elddata,emidata,h1,h2,B,Pd)

if nargin ~= 6
     error('Wrong number of input arguments')
end
n=length(elddata(:,1));

     Aeq=ones(1,n);
     a=elddata(:,1);
          b=elddata(:,2);
               c=elddata(:,3);
                    l=elddata(:,4);
                         u=elddata(:,5);
                          a1=emidata(:,1);
          b1=emidata(:,2);
               c1=emidata(:,3);
                           P=l;
                         for i=1:5
                             Pl=P'*B*P;
                             Pd1=Pd+Pl;
                             ll=diag(1-2*B*P);
                             A1=inv(ll)*(h1*a+h2*a1);
                              B1=inv(ll)*(h1*b+h2*b1);
                              H=2*diag(A1);
                              P=quadprog(H,B1,[],[],Aeq,Pd1,l,u);
                              Fcost=a'*(P.*P)+b'*P+sum(c);
                               Emi=a1'*(P.*P)+b1'*P+sum(c1);
                         end


text:
function [output] = mesh2edge(mesh)

% MESH2EDGE finds the edge lines from a triangulated mesh or the edge
% surfaces from a tetrahedral or hexahedral mesh. An edge is defined as an
% element that does not border any other element. This also implies that a
% closed triangulated surface has no edges.
%
% Use as
%   [edge] = mesh2edge(mesh)
%
% See also POLY2TRI

% Copyright (C) 2013-2020, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if isfield(mesh, 'tri')
  % make a list of all edges
  edge1 = mesh.tri(:, [1 2]);
  edge2 = mesh.tri(:, [2 3]);
  edge3 = mesh.tri(:, [3 1]);
  edge = cat(1, edge1, edge2, edge3);

elseif isfield(mesh, 'tet')
  % make a list of all triangles that form the tetraheder
  tri1 = mesh.tet(:, [1 2 3]);
  tri2 = mesh.tet(:, [2 3 4]);
  tri3 = mesh.tet(:, [3 4 1]);
  tri4 = mesh.tet(:, [4 1 2]);
  edge = cat(1, tri1, tri2, tri3, tri4);

elseif isfield(mesh, 'hex')
  % make a list of all "squares" that form the cube/hexaheder
  % FIXME should be checked, this is impossible without a drawing
  square1 = mesh.hex(:, [1 2 3 4]);
  square2 = mesh.hex(:, [5 6 7 8]);
  square3 = mesh.hex(:, [1 2 6 5]);
  square4 = mesh.hex(:, [2 3 7 6]);
  square5 = mesh.hex(:, [3 4 8 7]);
  square6 = mesh.hex(:, [4 1 5 8]);
  edge = cat(1, square1, square2, square3, square4, square5, square6);

end % isfield(mesh)

% soort all polygons in the same direction
% keep the original as "edge" and the sorted one as "sedge"
sedge = sort(edge, 2);

% % find the edges that are not shared -> count the number of occurences
% n = size(sedge,1);
% occurences = ones(n,1);
% for i=1:n
%   for j=(i+1):n
%     if all(sedge(i,:)==sedge(j,:))
%       occurences(i) = occurences(i)+1;
%       occurences(j) = occurences(j)+1;
%     end
%   end
% end
%
% % make the selection in the original, not the sorted version of the edges
% % otherwise the orientation of the edges might get flipped
% edge = edge(occurences==1,:);

% find the edges that are not shared
indx = findsingleoccurringrows(sedge);
edge = edge(indx, :);

% replace pnt by pos
mesh = fixpos(mesh);

% the naming of the edges in the output depends on what they represent
output.pos = mesh.pos;
if isfield(mesh, 'tri')
  % these have two vertices in each edge element
  output.line = edge;
  fprintf('reducing triangles to %d lines\n', size(edge,1));
elseif isfield(mesh, 'tet')
  % these have three vertices in each edge element
  output.tri = edge;
  fprintf('reducing tetraheders to %d triangles\n', size(edge,1));
elseif isfield(mesh, 'hex')
  % these have four vertices in each edge element
  output.poly = edge;
  fprintf('reducing hexaheders to %d polygons\n', size(edge,1));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION, see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1833#c12
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function indx = findsingleoccurringrows(X)
[X, indx] = sortrows(X);
sel  = any(diff([X(1,:)-1; X],1),2) & any(diff([X; X(end,:)+1],1),2);
indx = indx(sel);

function indx = finduniquerows(X)
[X, indx] = sortrows(X);
sel  = any(diff([X(1,:)-1; X],1),2);
indx = indx(sel);


text:
function net = net_deploy(net,regression)
%CNN_IMAGENET_DEPLOY  Deploy a CNN

isDag = isa(net, 'dagnn.DagNN') ;
if isDag
  dagRemoveLayersOfType(net, 'dagnn.Loss') ;
  dagRemoveLayersOfType(net, 'dagnn.DropOut') ;
else
  net = simpleRemoveLayersOfType(net, 'softmaxloss') ;
  net = simpleRemoveLayersOfType(net, 'pdist') ;
  net = simpleRemoveLayersOfType(net, 'dropout') ;
end

if regression == 0
    if isDag
        net.addLayer('prob', dagnn.SoftMax(), 'prediction', 'prob', {}) ;
    else
        net.layers{end+1} = struct('name', 'prob', 'type', 'softmax') ;
    end
end

if isDag
  dagMergeBatchNorm(net) ;
  dagRemoveLayersOfType(net, 'dagnn.BatchNorm') ;
else
  net = simpleMergeBatchNorm(net) ;
  net = simpleRemoveLayersOfType(net, 'bnorm') ;
end

if ~isDag
  net = simpleRemoveMomentum(net) ;
end

% Switch to use MatConvNet default memory limit for CuDNN (512 MB)
if ~isDag
  for l = simpleFindLayersOfType(net, 'conv')
    net.layers{l}.opts = removeCuDNNMemoryLimit(net.layers{l}.opts) ;
  end
else
  for name = dagFindLayersOfType(net, 'dagnn.Conv')
    l = net.getLayerIndex(char(name)) ;
    net.layers(l).block.opts = removeCuDNNMemoryLimit(net.layers(l).block.opts) ;
  end
end

% -------------------------------------------------------------------------
function opts = removeCuDNNMemoryLimit(opts)
% -------------------------------------------------------------------------
remove = false(1, numel(opts)) ;
for i = 1:numel(opts)
  if isstr(opts{i}) && strcmp(lower(opts{i}), 'CudnnWorkspaceLimit')
    remove([i i+1]) = true ;
  end
end
opts = opts(~remove) ;

% -------------------------------------------------------------------------
function net = simpleRemoveMomentum(net)
% -------------------------------------------------------------------------
for l = 1:numel(net.layers)
  if isfield(net.layers{l}, 'momentum')
    net.layers{l} = rmfield(net.layers{l}, 'momentum') ;
  end
end

% -------------------------------------------------------------------------
function layers = simpleFindLayersOfType(net, type)
% -------------------------------------------------------------------------
layers = find(cellfun(@(x)strcmp(x.type, type), net.layers)) ;

% -------------------------------------------------------------------------
function net = simpleRemoveLayersOfType(net, type)
% -------------------------------------------------------------------------
layers = simpleFindLayersOfType(net, type) ;
net.layers(layers) = [] ;

% -------------------------------------------------------------------------
function layers = dagFindLayersWithOutput(net, outVarName)
% -------------------------------------------------------------------------
layers = {} ;
for l = 1:numel(net.layers)
  if any(strcmp(net.layers(l).outputs, outVarName))
    layers{1,end+1} = net.layers(l).name ;
  end
end

% -------------------------------------------------------------------------
function layers = dagFindLayersOfType(net, type)
% -------------------------------------------------------------------------
layers = [] ;
for l = 1:numel(net.layers)
  if isa(net.layers(l).block, type)
    layers{1,end+1} = net.layers(l).name ;
  end
end

% -------------------------------------------------------------------------
function dagRemoveLayersOfType(net, type)
% -------------------------------------------------------------------------
names = dagFindLayersOfType(net, type) ;
for i = 1:numel(names)
  layer = net.layers(net.getLayerIndex(names{i})) ;
  net.removeLayer(names{i}) ;
  net.renameVar(layer.outputs{1}, layer.inputs{1}, 'quiet', true) ;
end

% -------------------------------------------------------------------------
function dagMergeBatchNorm(net)
% -------------------------------------------------------------------------
names = dagFindLayersOfType(net, 'dagnn.BatchNorm') ;
for name = names
  name = char(name) ;
  layer = net.layers(net.getLayerIndex(name)) ;

  % merge into previous conv layer
  player = dagFindLayersWithOutput(net, layer.inputs{1}) ;
  player = net.layers(net.getLayerIndex(player)) ;
  if ~isa(player.block, 'dagnn.Conv')
    error('Batch normalization cannot be merged as it is not preceded by a conv layer.') ;
  end

  filters = net.getParamIndex(player.params{1}) ;
  biases = net.getParamIndex(player.params{2}) ;
  multipliers = net.getParamIndex(layer.params{1}) ;
  offsets = net.getParamIndex(layer.params{2}) ;
  moments = net.getParamIndex(layer.params{3}) ;

  [filtersValue, biasesValue] = mergeBatchNorm(...
    net.params(filters).value, ...
    net.params(biases).value, ...
    net.params(multipliers).value, ...
    net.params(offsets).value, ...
    net.params(moments).value) ;

  net.params(filters).value = filtersValue ;
  net.params(biases).value = biasesValue ;
end

% -------------------------------------------------------------------------
function net = simpleMergeBatchNorm(net)
% -------------------------------------------------------------------------

for l = 1:numel(net.layers)
  if strcmp(net.layers{l}.type, 'bnorm')
    if ~strcmp(net.layers{l-1}.type, 'conv')
      error('Batch normalization cannot be merged as it is not preceded by a conv layer.') ;
    end
    [filters, biases] = mergeBatchNorm(...
      net.layers{l-1}.weights{1}, ...
      net.layers{l-1}.weights{2}, ...
      net.layers{l}.weights{1}, ...
      net.layers{l}.weights{2}, ...
      net.layers{l}.weights{3}) ;
    net.layers{l-1}.weights = {filters, biases} ;
  end
end

% -------------------------------------------------------------------------
function [filters, biases] = mergeBatchNorm(filters, biases, multipliers, offsets, moments)
% -------------------------------------------------------------------------
% wk / sqrt(sigmak^2 + eps)
% bk - wk muk / sqrt(sigmak^2 + eps)
a = multipliers(:) ./ moments(:,2) ;
b = offsets(:) - moments(:,1) .* a ;
biases(:) = biases(:) + b(:) ;
sz = size(filters) ;
numFilters = sz(4) ;
filters = reshape(bsxfun(@times, reshape(filters, [], numFilters), a'), sz) ;

repo: DBCNN


text:
function varargout=zerocros(varargin)
%
% For calling details please see v_zerocros.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_zerocros.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: zerocros.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_zerocros(varargin{:});
else
    v_zerocros(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   example 3 DOF planar robot.
%
%   Author: Arturo Gil. Universidad Miguel Hern�ndez de Elche.
%   email: arturo.gil@umh.es date:   05/03/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function robot = parameters()

robot.name='Example 3DOF planar arm';

%kinematic data DH parameters
robot.DH.theta='[q(1) q(2) q(3)]';
robot.DH.d='[0  0  0]';
robot.DH.a='[1  1  1]';
%robot.DH.a='[3  1  0.5]';
robot.DH.alpha='[0  0  0]';

%number of degrees of freedom
robot.DOF = 3;

%Jacobian matrix
robot.J=[];

robot.J=['[-a(1)*sin(q(1))-a(2)*sin(q(1)+q(2))-a(3)*sin(q(1)+q(2)+q(3))  -a(2)*sin(q(1)+q(2))-a(3)*sin(q(1)+q(2)+q(3)) -a(3)*sin(q(1)+q(2)+q(3));' ...
           'a(1)*cos(q(1))+a(2)*cos(q(1)+q(2))+a(3)*cos(q(1)+q(2)+q(3))  a(2)*cos(q(1)+q(2))+a(3)*cos(q(1)+q(2)+q(3))  a(3)*cos(q(1)+q(2)+q(3));' ...
           '               0                                  0  0;' ...
           '               0                                  0  0;' ...
           '               0                                  0  0;' ...
           '               1                                  1  1]'];

robot.kind=['R' 'R' 'R'];

%Function name to compute inverse kinematic
robot.inversekinematic_fn = 'inversekinematic_3dofplanar(robot, T)';
robot.directkinematic_fn = 'directkinematic(robot, q)';

%minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-180) deg2rad(180); %Axis 1, minimum, maximum
                deg2rad(-180) deg2rad(180);
                deg2rad(-180) deg2rad(180)]; %Axis 2, minimum, maximum

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = []; %empty, not available

robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time
% end effectors maximum velocity
robot.linear_velmax = 0; %m/s, example, not available

%base reference system
robot.T0 = eye(4);

%INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
%position, velocity and acceleration
robot=init_sim_variables(robot);
robot.path = pwd;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GRAPHICS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%read graphics files
robot.graphical.has_graphics=1;
robot.graphical.color = [255 20 40]./255;
%for transparency
robot.graphical.draw_transparent=0;
%draw DH systems
robot.graphical.draw_axes=1;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=1;
%adjust for a default view of the robot
robot.axis=[-3.5 3.5 -3.5 3.5 0 1];
robot = read_graphics(robot);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DYNAMIC PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
robot.has_dynamics=1;

%consider friction in the computations
robot.dynamics.friction=0;

%link masses (kg)
robot.dynamics.masses=[1 1 1];

%COM of each link with respect to own reference system
robot.dynamics.r_com=[-0.5      0         0; %(rx, ry, rz) link 1
                      -0.5      0         0;
                      -0.5      0         0];%(rx, ry, rz) link 2

%link masses
m1 = robot.dynamics.masses(1);
m2 = robot.dynamics.masses(2);
m3 = robot.dynamics.masses(3);

a=eval(robot.DH.a);
L1 = a(1);
L2 = a(2);
L3 = a(3);

%Momentos de inercia de cada eslabon.
% Ixx	Iyy	Izz	Ixy	Iyz	Ixz, por cada fila
robot.dynamics.Inertia=[0   m1*L1^2/3   m1*L1^2/3    0	0	0;
                        0   m2*L2^2/3   m2*L2^2/3    0	0	0;
                        0   m3*L3^2/3   m3*L3^2/3    0	0	0];

%Actuator rotor inertia
robot.motors.Inertia=[0 0 0];
%Reduction ratio motor/joint speed
robot.motors.G=[1 1 1];
%consider friction
robot.friction=0;
%Viscous friction factor, motor referred
%robot.B = [1e-3  1e-3 1e-3];
robot.motors.Viscous = [10  10 10];
%Coulomb friction, motor referred
robot.motors.Coulomb = [0 0;
            0 0;
            0 0];

repo: ARTE


text:
% Fig. 5.8   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami

n6=1;
d6=[1 8 32 0];
pzmap(n6,d6)
title('Fig.5.8 The locus crosses the j\omega axis at *')
axis([-10 6 -6 6] )
hold on
r=roots([1 0 32]);
 plot(r,'*')
z=0:.1:.9;
 wn= 1:6;
 sgrid(z, wn)
 hold off


text:
function out_fwrite_edf(sFile, sfid, SamplesBounds, ChannelsRange, F)
% OUT_FWRITE_EDF: Write a block of recordings from a EDF 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: Martin Cousineau, 2017
%          Francois Tadel, 2019

% ===== PARSE INPUTS =====
[nSignals, nSamples] = size(F);
if isempty(SamplesBounds)
    SamplesBounds = [0, nSamples];
end
if isempty(ChannelsRange)
    ChannelsRange = [1, nSignals];
end

fseek(sfid, 0, 'eof');

% Get the gains of the channels for all the non-Annotation channels
iChanGain = setdiff(1:length(sFile.header.signal), sFile.header.annotchan);
% Apply channel gains before converting to integer
chgain = [sFile.header.signal(iChanGain).unit_gain] ./ ...
            ([sFile.header.signal(iChanGain).physical_max] - [sFile.header.signal(iChanGain).physical_min]) .* ...
            ([sFile.header.signal(iChanGain).digital_max]  - [sFile.header.signal(iChanGain).digital_min]);
F = bst_bsxfun(@times, F, chgain');

% Convert to 2-byte integer in 2's complement
F = int16(F);
negF = F < 0;
F(negF) = bitcmp(abs(F(negF))) + 1;

% Prepare annotations if any.
if sFile.header.annotchan >= 0
    annotations    = 1;
    nAnnots        = numel(sFile.header.annotations);
    nSamplesReal   = round((sFile.prop.times(2) - sFile.prop.times(1)) .* sFile.prop.sfreq);
    nSamplesAnnots = sFile.header.signal(sFile.header.annotchan).nsamples;
    annotThreshold = floor((1:nAnnots) / nAnnots * nSamplesReal);
    annotBounds    = [0, 0];

    % Insert annotation in this record only if it contains the required
    % cutoff sample threshold
    for iThr = 1:nAnnots
        if annotThreshold(iThr) >= SamplesBounds(1) && annotThreshold(iThr) <= SamplesBounds(2)
            if annotBounds(1) < 1
                annotBounds(1) = iThr;
            end
            annotBounds(2) = iThr;
        end
    end

    if annotBounds(2) < 1
        annotsList = [];
    else
        annotsList = sFile.header.annotations(annotBounds(1) : annotBounds(2));
    end

    nAnnots     = numel(annotsList);
    nextAnnot   = 1;
else
    annotations = 0;
end

% Write to file record per record
nSamplesPerRecord = sFile.prop.sfreq * sFile.header.reclen;
nRecords          = ceil((SamplesBounds(2) - SamplesBounds(1)) / nSamplesPerRecord);
ncount            = 0;
bounds            = [1, nSamplesPerRecord];
timeOffset        = SamplesBounds(1) / sFile.prop.sfreq;

for iRec = 1:nRecords
    % Special case when we don't have enough data to fill the last record
    if bounds(2) > nSamples
        if iRec ~= nRecords
            error('Ran out of data before last record.');
        end
        writeZeros = bounds(2) - nSamples;
        bounds(2)  = nSamples;
    else
        writeZeros = 0;
    end

    % Write data
    for iSig = ChannelsRange(1):ChannelsRange(2)
        ncount = ncount + fwrite(sfid, F(iSig, floor(bounds(1)):floor(bounds(2))), 'int16');

        % Fill rest of the record with 0s if required
        if writeZeros
            fwrite(sfid, zeros(writeZeros, 1), 'int16');
        end
    end

    % Write annotations if any, split by records
    if annotations
        bytesLeft = nSamplesAnnots * 2;

        % The first annotation specifies the time offset
        bytesLeft = bytesLeft - fprintf(sfid, '+%f%c%c%c', timeOffset, char(20), char(20), char(0));

        % Write as many annotations as possible in current record
        while nextAnnot <= nAnnots && bytesLeft >= length(annotsList{nextAnnot})
            bytesLeft = bytesLeft - fprintf(sfid, '%s', annotsList{nextAnnot});
            nextAnnot = nextAnnot + 1;
        end

        % Fill remaining of record with 0-bytes.
        fprintf(sfid, '%s', repmat(char(0), 1, bytesLeft));
    end

    % Get ready for next record
    bounds     = bounds + nSamplesPerRecord;
    timeOffset = timeOffset + sFile.header.reclen;
end

% Check number of values written
if (ncount ~= numel(F))
    error('Error writing data to file.');
end

author: brainstorm-tools


text:
function linplus_test152 ( )

%*****************************************************************************80
%
%% TEST152 tests R8BLT_INDICATOR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 6;
  ml = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST152\n' );
  fprintf ( 1, '  R8BLT_INDICATOR sets up a R8BLT indicator matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N     = %d\n', n );
  fprintf ( 1, '  Lower bandwidth ML = %d\n', ml );
%
%  Set the matrix.
%
  a = r8blt_indicator ( n, ml );

  r8blt_print ( n, ml, a, '  The indicator matrix:' );

  return
end