text:
function [rgb,ucs,status,RGB,UCS] = maxdistcolor_view(N,fun,varargin)
% Create a figure for interactive generation and display of MAXDISTCOLOR colors.
%
% (c) 2017-2020 Stephen Cobeldick
%
% This function has exactly the same inputs and outputs as MAXDISTCOLOR.
% See MAXDISTCOLOR for descriptions of the required and optional arguments.
%
% See also MAXDISTCOLOR MAXDISTCOLOR_DEMO SRGB_TO_CAM02UCS SRGB_TO_CIELAB
% SRGB_TO_OSAUCS COLORNAMES COLORMAP RGBPLOT AXES SET LINES PLOT

%% New Figure %%
%
figH = figure('Units','pixels', 'ToolBar','figure', 'NumberTitle','off',...
	'HandleVisibility','off', 'Name',mfilename(), 'Visible','on');
figP = get(figH, 'Position');
figY = get(figH, 'Pointer');
set(figH, 'Pointer','watch')
%
iniH = uicontrol(figH, 'Units','pixels', 'Style','text', 'HitTest','off',...
	'Visible','on',	'String','Initializing the figure... please wait.');
iniX = get(iniH,'Extent');
set(iniH,'Position',[figP(3:4)/2-iniX(3:4)/2,iniX(3:4)])
%
drawnow()
%
% First call to check input arguments and define all output arguments:
[bgr,ucs,status,RGB,UCS] = maxdistcolor(N,fun,varargin{:});
%
opts = status.options;
%
inw = true;
ohm = pow2([opts.bitR,opts.bitG,opts.bitB])-1;
%
if isempty(opts.exc)
	opts.exc = nan(0,3);
end
if isempty(opts.inc)
	opts.inc = nan(0,3);
end
%
delete(iniH)
%
%% Parameters %%
%
gap = 7;
% UI colors:
[R,G,B] = ndgrid(0:1/8:1);
smp = [R(:),G(:),B(:)];
clear R G B
txc = [0.7,0.3,0.2];
bgd = [];
fgd = [];
fBackFore(false)
% Slider steps:
bStp = [1,2]; % bits
cStp = [0.005,0.1]; % chroma
lStp = [0.005,0.1]; % lightness
% Number properties:
nRng = [1,64];
nStp = [1,5];
nFun = @(n) max(nRng(1),min(nRng(2),n));
% Table properties:
tTxt = {'exc'; 'inc'};
% Menu properties:
mTxt = {'plot'; 'sort'; 'path'; 'disp'};
mSrt = {{'none','farthest','hue','zip','lightness','a','b'},{'maxmin','minmax','longest','shortest'}};
mStr = {''; [mSrt{:}]; {'open','closed'}; {'off','time','summary','verbose'}};
% Slider properties:
sTxt = {'Cmax'; 'Cmin'; 'Lmax'; 'Lmin'; 'bitR'; 'bitG'; 'bitB'; 'start'};
sInt = [ false;  false;  false;  false;   true;   true;   true;   false];
sRng = [   0,1;    0,1;    0,1;    0,1;    1,8;    1,8;    1,8;   0,360];
sStp = [  cStp;   cStp;   lStp;   lStp;   bStp;   bStp;   bStp;    1,10];
sTwo = [    +1;     -1;     +2;     -2;      0;      0;      0;       0];
sFun = @(n,k) max(sRng(k,1),min(sRng(k,2),n));
sStr = @(x) sprintf('%.4g',round(1000*x)/1000);
% Colorspace axes:
f2s = upper(regexp(func2str(fun),'(?<=_to_)\w+','match','once'));
tmp = struct('CAM02UCS',{'J''','a''','b'''}, 'CIELAB',{'L*','a*','b*'},...
	'DIN99',{'L_{99}','a_{99}','b_{99}'}, 'OSAUCS',{'L','j','g'});
if isfield(tmp,f2s)
	zyx = {tmp.(f2s)};
else
	zyx = {'L?','a?','b?'};
	f2s = '???';
end
% Maximum chroma:
tmp = fun(smp);
mxc = max(max(abs(tmp(:,2:3))));
%
%% 3D Axes %%
%
az = [];
el = [];
%
%% Interactive Graphics Objects %%
%
% Get text width:
cnc = {'Red','Green','Blue'};
tmp = uicontrol(figH, 'Style','text', 'Units','pixels', 'String',cnc);
txw = get(tmp,'Extent');
txh = txw(4)/numel(cnc);
txw = num2cell(ones(1,3)*(7+txw(3)));
delete(tmp)
% Add drop-down menus:
mTxN = numel(mTxt);
mAxH = axes('Parent',figH, 'Units','pixels', 'Visible','off', 'View',[0,90],...
	'HitTest','off', 'Xlim',[0,1], 'Ylim',[0,2*mTxN]);
mTxH = text(zeros(1,mTxN),2*mTxN-1:-2:1, mTxt, 'parent',mAxH, 'Color',txc,...
	'VerticalAlignment','bottom','HorizontalAlignment','left');
set(mTxH(1),'Color',txc([3,2,1]));
for k = mTxN:-1:1
	if k>1
		idm = find(strcmpi(opts.(mTxt{k}),mStr{k}));
	else
		idm = 1;
	end
	mUiH(k) = uicontrol(figH, 'Style','popupmenu', 'Units','pixels',...
		'String',mStr{k}, 'Callback',{@fOptMenu,k}, 'Value',idm);
end
% Add horizontal sliders:
sTxN = numel(sTxt);
for k = sTxN:-1:1
	val = opts.(sTxt{k});
	sTxH(k) = uicontrol(figH, 'Units','pixels', 'Style','text',...
		'String',sTxt{k}, 'HorizontalAlignment','left', 'ForegroundColor',txc);
	sVaH(k) = uicontrol(figH, 'Units','pixels', 'Style','edit',...
		'String',sStr(val), 'HorizontalAlignment','left', 'Callback',{@fOptEdit,k});
	sUiH(k) = uicontrol(figH, 'Units','pixels', 'Style','slider', 'Value',1,...
		'Min',sRng(k,1), 'Max',sRng(k,2), 'Value',sFun(val,k),...
		'SliderStep',sStp(k,:)./diff(sRng(k,:)), 'Callback',@fUpDtMap);
	addlistener(sUiH(k), 'Value', 'PostSet',@(o,e)fOptSlide(o,e,k));
end
% Add tables:
tTxN = numel(tTxt);
for k = tTxN:-1:1
	tmp = opts.(tTxt{k});
	tmp(end+1:N,:) = NaN;
	tUiH(k) = uitable(figH, 'Units','pixels', 'Data',tmp, 'ColumnWidth',txw,...
		'ColumnName',cnc, 'ColumnEditable',true, 'CellEditCallback',{@fCellEdit,k});
	tTxH(k) = uicontrol(figH, 'Style','text', 'Units','pixels','ForegroundColor',txc,...
		'String',tTxt{k}, 'HorizontalAlignment','left');
	tCbH(k) = uicontrol(figH, 'Style','checkbox', 'Units','pixels','String','X',...
		'Callback',{@fCheckBox,k});
end
% Add colorbar:
bAxH = axes('Parent',figH, 'Units','pixels', 'Visible','off',...
	'HitTest','off', 'Xlim',[0.5,1.4], 'Ylim',[0,N]+0.5, 'View',[0,90],...
	'YDir','normal', 'XTick',[], 'YTick',1:N, 'Box','off');
bImH = image('CData',permute(bgr,[1,3,2]), 'Parent',bAxH);
% Add number slider:
nVaH = uicontrol(figH, 'Units','pixels', 'Style','edit', 'String',sprintf('%d',N),...
	'HorizontalAlignment','center', 'Callback',@fNumEdit);
nUiH = uicontrol(figH, 'Units','pixels', 'Style','slider', 'Value',nFun(N),...
	'Min',nRng(1), 'Max',nRng(2), 'SliderStep',nStp./diff(nRng), 'Callback',@fUpDtMap);
addlistener(nUiH, 'Value', 'PostSet',@fNumSlide);
% Add EVAL button:
eUiH = uicontrol(figH, 'Units','pixels', 'Style','pushbutton',...
	'String','pause', 'ForegroundColor',txc([3,2,1]), 'Callback',@fEvalFun);
eUiX = true;
eUiC = get(eUiH,'BackgroundColor');
%
%% Main Plot Objects %%
%
ndx = 'Colormap Index';
eds = sprintf('Euclidean Distance (%s)',f2s);
%
axp = {'Units','normalized', 'NextPlot','replacechildren', 'Clipping','off',...
	'Color','none', 'XColor',fgd, 'YColor',fgd, 'ZColor',fgd, 'UserData'};
pnp = {figH, 'Units','pixels', 'BorderType','none', 'Title','', 'BackgroundColor',bgd};
%
pAxH(10) = axes('Parent',uipanel(pnp{:}), axp{:},4);
pAxS(10) = struct('title','RGB Gamut (alphashape) [slow]',...
	'X',zyx{3}, 'Y',zyx{2}, 'Z',zyx{1}, 'fun',@fAlphaTry);
%
pAxH(9) = axes('Parent',uipanel(pnp{:}), axp{:},3);
pAxS(9) = struct('title','RGB Gamut (point cloud) [fast]',...
	'X',zyx{3}, 'Y',zyx{2}, 'Z',zyx{1}, 'fun',@fPointCloud);
%
pAxH(8) = axes('Parent',uipanel(pnp{:}), axp{:},3);
pAxS(8) = struct('title','Colors with RGB Cube',...
	'X',zyx{3}, 'Y',zyx{2}, 'Z',zyx{1}, 'fun',@fCubeRGB);
%
pAxH(7) = axes('Parent',uipanel(pnp{:}), axp{:},3);
pAxS(7) = struct('title','Colors with Sort Path',...
	'X',zyx{3}, 'Y',zyx{2}, 'Z',zyx{1}, 'fun',@fSortPath);
%
pAxH(6) = axes('Parent',uipanel(pnp{:}), axp{:},2, 'Visible','off');
pAxS(6) = struct('title','Matrix Scatter Plot',...
	'X','', 'Y','', 'Z','', 'fun',@fScatterMat);
%
pAxH(5) = axes('Parent',uipanel(pnp{:}), axp{:},2);
pAxS(5) = struct('title','Adjacent Color Matrix',...
	'X',ndx, 'Y',ndx, 'Z','', 'fun',@fAdjaNode);
%
pAxH(4) = axes('Parent',uipanel(pnp{:}), axp{:},2, 'XTick',[]);
pAxS(4) = struct('title','Bands with RGB Values',...
	'X','', 'Y',ndx, 'Z','', 'fun',@fRgbValues);
%
pAxH(3) = axes('Parent',uipanel(pnp{:}), axp{:},2, 'XTick',[]);
pAxS(3) = struct('title','Bands with Color Names',...
	'X','', 'Y',ndx, 'Z','', 'fun',@fBandName);
%
pAxH(2) = axes('Parent',uipanel(pnp{:}), axp{:},3);
pAxS(2) = struct('title','Euclidean Distance (3D Plot)',...
	'X',ndx, 'Y',ndx, 'Z',eds, 'fun',@fEuclDist);
%
pAxH(1) = axes('Parent',uipanel(pnp{:}), axp{:},2);
pAxS(1) = struct('title','Euclidean Distance (2D Plot)',...
	'X',ndx, 'Y',eds, 'Z','', 'fun',@fEuclDist);
%
pPnH = get(pAxH,'Parent');
pPnH = [pPnH{:}];
pScH = [];
%
for k = 1:numel(pAxH)
	xlabel(pAxH(k), pAxS(k).X, 'Color',fgd)
	ylabel(pAxH(k), pAxS(k).Y, 'Color',fgd)
	zlabel(pAxH(k), pAxS(k).Z, 'Color',fgd)
	title( pAxH(k), pAxS(k).title, 'Color',fgd, 'Visible','on')
	if get(pAxH(k), 'UserData')>2
		% Axes 3D view:
		view(pAxH(k),3)
		% Zlabel orientation for short strings:
		if numel(pAxS(k).Z)<7
			set(get(pAxH(k), 'ZLabel'), 'Rotation',0, 'HorizontalAlignment','right')
		end
	end
end
%
% Optional linking of 3D axes:
xud = cellfun(@(m)isequal(m,zyx{3}),{pAxS.X});
linkprop(pAxH(xud), 'View');
%
%% Initialize GUI %%
%
rgb2str = @(f,d,m)cellfun(@(v)sprintf(f,v),reshape(num2cell(round(m*d)/d,2),1,[]),'uni',0);
arrayfun(@(s,a) s.fun(a), pAxS, pAxH)
set(pPnH(:), 'Visible','off', 'HitTest','off')
set(pPnH(1), 'Visible','on', 'HitTest','on')
set(mUiH(1), 'String',{pAxS.title})
set(figH, 'Pointer',figY, 'ResizeFcn',@fSizeObj)
fNumTick()
fSizeObj()
%
%% Main Plot Functions %%
%
	function fOrient3D(axh,az,el)
		axis(axh,'equal')
		grid(axh,'on')
		view(axh,az,el)
		mat = fun([0,0,0;1,1,1]);
		set(axh, 'XGrid','on', 'YGrid','on', 'ZGrid','on',...
			'XLim',[-mxc,mxc], 'YLim',[-mxc,mxc], 'ZLim',[mat(1),mat(2)])
	end
%
	function fAlphaTry(~) % ALPHASHAPE is slow, run only on demand.
		plv = get(mUiH(strcmpi('plot',mTxt)), 'Value');
		if eUiX && inw && 4==get(pAxH(plv), 'UserData')
			set(figH, 'Pointer','watch')
			drawnow()
			fAlphaRGB(pAxH(plv))
			set(figH, 'Pointer',figY)
			inw = false;
		end
	end
	function fAlphaRGB(axh) % Show the RGB gamut using alphashape.
		delete(get(axh, 'Children'))
		[az,el] = view(axh);
		% Get gamut boundary:
		nrw = size(UCS,1);
		mxn = 1e4; % max elements for ALPHASHAPE: more elements -> slower runtime.
		stp = ceil(nrw/mxn);
		vec = [];
		for idk = 1:stp
			idu = idk:stp:nrw;
			try
				bnd = fGetBound(idu);
			catch %#ok<CTCH>
				text(0.5, 0.5, 0.5,...
					{'ALPHASHAPE not found','(requires R2014b or later)'},...
					'Parent',axh, 'HorizontalAlignment','center',...
					'FontWeight','bold', 'Color',fgd)
				return
			end
			vec = union(vec,idu(bnd(:)));
		end
		bnd = fGetBound(vec);
		% Show RGB gamut:
		patch('Faces',bnd, 'Vertices',UCS(vec,[3,2,1]), 'Parent',axh,...
			'EdgeColor','none', 'FaceColor','interp', 'FaceVertexCData',RGB(vec,:));
		% Show color nodes:
		hold(axh,'on')
		scatter3(ucs(:,3),ucs(:,2),ucs(:,1), 13, fgd, 'filled', 'Parent',axh)
		% Orient axes:
		fOrient3D(axh,az,el)
	end
	function bnd = fGetBound(idu) % Get boundary of the node cloud.
		shp = alphaShape(UCS(idu,[3,2,1]),Inf); % requires R2014b or later.
		shp.Alpha = shp.criticalAlpha('all-points');
		bnd = shp.boundaryFacets();
	end
%
	function fPointCloud(axh) % Show the RGB gamut using a point cloud.
		delete(get(axh, 'Children'))
		[az,el] = view(axh);
		% Plot point cloud:
		npt = 1e4;
		I = unique(round(linspace(1,size(UCS,1),npt)));
		scatter3(UCS(I,3),UCS(I,2),UCS(I,1), 3, RGB(I,:), 'filled', 'Parent', axh)
		hold(axh,'on')
		scatter3(ucs(:,3),ucs(:,2),ucs(:,1), 13, fgd, 'filled', 'Parent',axh)
		% Orient axes:
		fOrient3D(axh,az,el)
	end
%
	function fSortPath(axh) % Colors shown in UCS space, with path.
		delete(get(axh, 'Children'))
		[az,el] = view(axh);
		% Show color nodes:
		n2s = cellstr(strjust(num2str((1:N).'),'left'));
		scatter3(ucs(:,3),ucs(:,2),ucs(:,1), 256, bgr, 'filled', 'Parent',axh)
		text(ucs(:,3),ucs(:,2),ucs(:,1), n2s, 'Parent',axh,...
			'Color',fgd, 'HorizontalAlignment','center')
		hold(axh,'on')
		% Show path:
		ipc = strcmpi(opts.path,'closed');
		mat = ucs([1:N,1:+ipc],:);
		plot3(mat(:,3),mat(:,2),mat(:,1),'-', 'Parent',axh, 'Color',fgd)
		% Show distances:
		vec = sqrt(sum(diff(mat,1,1).^2,2));
		[mxv,mxi] = max(vec);
		[mnv,mni] = min(vec);
		mxp = mean(mat(mxi:mxi+1,:),1);
		mnp = mean(mat(mni:mni+1,:),1);
		text(mxp(3),mxp(2),mxp(1), sprintf('max:%.5g',mxv), 'Parent',axh, 'Color',fgd)
		text(mnp(3),mnp(2),mnp(1), sprintf('min:%.5g',mnv), 'Parent',axh, 'Color',fgd)
		% Orient axes:
% 		axis(axh,'equal')
% 		view(axh,az,el)
		fOrient3D(axh,az,el)
	end
%
	function fCubeRGB(axh) % Colors shown in UCS space, with RGB cube.
		delete(get(axh, 'Children'))
		[az,el] = view(axh);
		% Show color nodes:
		n2s = cellstr(strjust(num2str((1:N).'),'left'));
		scatter3(ucs(:,3),ucs(:,2),ucs(:,1), 256, bgr, 'filled', 'Parent',axh)
		text(ucs(:,3),ucs(:,2),ucs(:,1), n2s, 'Parent',axh,...
			'Color',fgd, 'HorizontalAlignment','center')
		% Show excluded nodes:
		if isinteger(opts.exc)
			exc = bsxfun(@rdivide,double(opts.exc),ohm);
		else
			exc = opts.exc;
		end
		cxe = fun(exc);
		hold(axh,'on')
		scatter3(cxe(:,3),cxe(:,2),cxe(:,1), 256, exc, 'filled', 'Parent',axh)
		text(cxe(:,3),cxe(:,2),cxe(:,1), 'X', 'Parent',axh,...
			'Color',fgd, 'HorizontalAlignment','center')
		% Show outline of RGB cube:
		M = 23;
		[X,Y,Z] = ndgrid(linspace(0,1,M),0:1,0:1);
		mat = fun([X(:),Y(:),Z(:);Y(:),Z(:),X(:);Z(:),X(:),Y(:)]);
		X = reshape(mat(:,3),M,[]);
		Y = reshape(mat(:,2),M,[]);
		Z = reshape(mat(:,1),M,[]);
		line(X,Y,Z, 'Color',fgd, 'Parent',axh)
		% Orient axes:
% 		axis(axh,'equal')
% 		view(axh,az,el)
		fOrient3D(axh,az,el)
	end
%
	function fScatterMat(axh) % Matrix scatter plot.
		uih = get(axh, 'Parent');
		ise = isempty(pScH);
		for rr = 3:-1:1
			for cc = 3:-1:1
				if ise
					pos = [rr-1,(3-cc)*0.93,1,0.93]./3;
					pScH(rr,cc) = axes('Parent',uih, 'Visible','on',...
						'Units','normalized', 'OuterPosition',pos,...
						'NextPlot','replacechildren',...
						'Color','none', 'XColor',fgd, 'YColor',fgd); %#ok<LAXES>
				else
					delete(get(pScH(rr,cc), 'Children'))
				end
				if rr==cc
					text(0.5,0.5, zyx{rr}, 'Parent',pScH(rr,cc),...
						'FontWeight','bold', 'Color',fgd)
				else
					scatter(pScH(rr,cc), ucs(:,rr),ucs(:,cc), 32, bgr, 'filled');%, 'Parent',axh)
				end
			end
		end
		%set(pScH, 'XLimMode','auto', 'YLimMode','auto')
		set(pScH([1,5,9]), 'XLim',0:1, 'YLim',0:1, 'Visible','off')
		set(pScH, 'XColor',fgd, 'YColor',fgd, 'ZColor',fgd)
	end
%
	function fAdjaNode(axh) % Show matrix with all colors adjacent.
		delete(get(axh, 'Children'))
		[idr,idc] = ndgrid(0:N);
		F = (1:N*(N+1)).';
		F(N+1:N+1:end) = [];
		image('Parent',axh, 'CData',repmat(permute(bgr,[1,3,2]),[1,N,1]));
		patch('Parent',axh, 'Faces',[F,F+N+1,F+N+2], 'Vertices',0.5+[idr(:),idc(:)],...
			'FaceColor','flat', 'FaceVertexCData',bgr(mod(F,N+1),:), 'EdgeColor','none');
	end
%
	function fBandName(axh) % Show bands with colornames (if COLORNAMES is available).
		delete(get(axh, 'Children'))
		image(permute(bgr,[1,3,2]), 'Parent',axh)
		try
			cnm = colornames('CSS',bgr);
		catch %#ok<CTCH>
			cnm = repmat({'COLORNAMES not found: download FEX #48155'},1,N);
		end
		text(ones(1,N), 1:N, cnm, 'Parent',axh, 'Color',fgd,...
			'BackgroundColor',bgd, 'HorizontalAlignment','center')
	end
%
	function fRgbValues(axh) % Show bands with RGB values (floating point and integer).
		delete(get(axh, 'Children'))
		image(permute(bgr,[1,3,2]), 'Parent',axh)
		% Integer:
		cnm = rgb2str('[%d,%d,%d] =',1,bsxfun(@times,bgr,ohm));
		text(ones(1,N), 1:N, cnm, 'Parent',axh, 'Color',fgd,...
			'BackgroundColor',bgd, 'HorizontalAlignment','right')
		% Normalized:
		cnm = rgb2str(' [%.5g,%.5g,%.5g]',1e5,bgr);
		text(ones(1,N), 1:N, cnm, 'Parent',axh, 'Color',fgd,...
			'BackgroundColor',bgd, 'HorizontalAlignment','left')
	end
%
	function fEuclDist(axh) % Show Euclidean distances (2D or 3D).
		delete(get(axh, 'Children'))
		aud = get(axh, 'UserData');
		sgc = {'Parent',axh, 'LineWidth',2.8, 'Marker','o'};
		for idk = 1:N
			dst = sqrt(sum(bsxfun(@minus,ucs,ucs(idk,:)).^2,2));
			switch aud
				case 2
					scatter(1:N, dst, 123, bgr,...
						sgc{:}, 'MarkerFaceColor',bgr(idk,:), 'Parent',axh);
				case 3
					scatter3(idk*ones(1,N), 1:N, dst, 123, bgr,...
						sgc{:}, 'MarkerFaceColor',bgr(idk,:), 'Parent',axh);
				otherwise
					error('SC:maxdistcolor_view:TooManyDimensions',...
						'Sorry, I don''t know how to plot %d dimensions.',aud)
			end
			hold(axh,'on')
		end
		%
		if aud==2
			dst = status.minDistOutput;
			text(N+0.5,dst,sprintf(' %#.5g',dst), 'Parent',axh,...
				'Color',fgd, 'BackgroundColor','none')
		end
	end
%
%% Callback Functions %%
%
	function fEvalFun(obj,~) % Turn evaluation on and off.
		eUiX = ~eUiX;
		if eUiX % eval
			set(obj,'String','pause', 'BackgroundColor',eUiC)
		else % paused
			set(obj,'String','eval', 'BackgroundColor',[1,1,0])
		end
		fUpDtMap()
	end
	function fUpDtMap(~,~) % Update colormap using options structure.
		if eUiX
			% Generate colormap:
			set(figH, 'Pointer','watch')
			drawnow()
			[bgr,ucs,status,RGB,UCS] = maxdistcolor(N,fun,opts);
			set(figH, 'Pointer',figY)
			inw = true;
			ohm = pow2([opts.bitR,opts.bitG,opts.bitB])-1;
			% Update colorbar:
			set(bImH, 'CData',permute(bgr,[1,3,2]))
			% Update main plots:
			arrayfun(@(s,h) s.fun(h), pAxS, pAxH)
		end
	end
%
	function fNumEdit(obj,~) % Number edit box callback.
		str = get(obj, 'String');
		if all(isstrprop(str,'digit')) && sscanf(str,'%d')
			N = sscanf(str,'%d');
			set(nUiH, 'Value',nFun(N))
			fNumTick()
			fUpDtMap()
		else
			set(obj,'String',sprintf('%d',N))
		end
	end
	function fNumSlide(~,evt) % Number slider listener callback.
		try
			N = round(get(evt,'NewValue'));
		catch %#ok<CTCH>
			N = round(evt.AffectedObject.Value);
		end
		set(nVaH, 'String',sprintf('%d',N))
		fNumTick()
	end
	function fNumTick() % Number adjust limits and tickmarks.
		fPermSort()
		% Colorbar limits:
		set(bAxH, 'Ylim',[0,N]+0.5, 'YTick',1:N)
		% Table sizes:
		for idk = 1:tTxN
			new = get(tUiH(idk), 'Data');
			new(end+1:N,:) = NaN;
			new(N+1:end,:) = [];
			set(tUiH(idk), 'Data',new)
		end
		% Main plot axes limits:
		for idk = 1:numel(pAxH)
			for idc = 'XYZ'
				if strcmpi(pAxS(idk).(idc),ndx)
					set(pAxH(idk), [idc,'Lim'],[0,N]+0.5, [idc,'Tick'],1:N);
				end
			end
		end
	end
%
	function fOptEdit(obj,~,idk) % Options edit box callback.
		str = get(obj, 'String');
		new = str2double(str);
		if ~isreal(new) || (sInt(idk) && fix(new)~=new) || new<sRng(idk,1) || new>sRng(idk,2)
			new = NaN;
		end
		if isnan(new)
			set(obj, 'String',sStr(opts.(sTxt{idk})))
		else
			opts.(sTxt{idk}) = new;
			set(sUiH(idk), 'Value',new)
			fUpDtMap()
		end
	end
	function fOptSlide(~,evt,idk) % Options slider listener callback.
		try
			new = get(evt,'NewValue');
		catch %#ok<CTCH>
			new = evt.AffectedObject.Value;
		end
		if sInt(idk)
			new = round(new);
			set(sUiH(idk), 'Value',new)
		end
		opts.(sTxt{idk}) = new;
		set(sVaH(idk), 'String',sStr(new))
		% Move paired slider:
		if sTwo(idk)
			idt = setdiff(find(abs(sTwo)==abs(sTwo(idk))),idk);
			old = get(sUiH(idt), 'Value');
			if ((new-old) * diff(sign(sTwo([idt,idk])))) < 0
				set(sUiH(idt), 'Value',new)
			end
		end
	end
%
	function fOptMenu(obj,~,idk) % Options menu callback.
		idv = get(obj,'Value');
		if idk>1
			opts.(mTxt{idk}) = mStr{idk}{idv};
			fUpDtMap()
		else
			set(pPnH,     'Visible','off', 'HitTest','off')
			set(pPnH(idv),'Visible','on',  'HitTest','on')
		end
		fAlphaTry()
	end
%
	function fCellEdit(obj,~,idk) % Options table cell callback.
		new = get(obj,'Data');
		isn = sum(isnan(new),2);
		idz = isn==0;
		if all(idz|(isn==3))
			new = new(idz,:);
			if get(tCbH(idk),'Value') % uint
				pwr = pow2(max(3,ceil(log2(max(log2(ohm+1))))));
				new = cast(new,sprintf('uint%d',pwr));
			end
			opts.(tTxt{idk}) = new;
			fBackFore(true)
			fUpDtMap()
		end
	end
	function fCheckBox(obj,~,idk) % Options table checkbox callback.
		old = get(tUiH(idk),'Data');
		if get(obj,'Value')
			% float->uint
			new = round(bsxfun(@times,old,ohm));
		else
			% uint->float
			new = bsxfun(@rdivide,old,ohm);
		end
		set(tUiH(idk),'Data',new);
	end
%
	function fPermSort() % No permutation sorting if N>9.
		ids = strcmpi('sort',mTxt);
		idn = get(mUiH(ids), 'Value');
		if N>9
			if idn > numel(mSrt{1})
				idn = 1;
			end
			opts.sort = mSrt{1}{idn};
			set(mUiH(ids), 'String',mSrt{1}, 'Value',idn)
		else
			set(mUiH(ids), 'String',mStr{ids})
		end
	end
%
	function fBackFore(chg) % Background and Foreground colors.
		if isempty(opts.exc)
			bgd = [1,1,1];
			fgd = [0,0,0];
		else
			new = opts.exc(1,:);
			if isfloat(new)
				bgd = double(new);
			else
				bgd = double(new) ./ ohm;
			end
			% Define forground as farthest color from background:
			[~,idf] = max(sum(bsxfun(@minus,fun(bgd),fun(smp)).^2,2));
			fgd = smp(idf,:);
		end
		if chg
			uih = get(pAxH, 'Parent');
			set([uih{:}], 'BackgroundColor',bgd)
			uih = get(pAxH, 'Title');
			set([uih{:}], 'Color',fgd)
			set(pAxH, 'XColor',fgd, 'YColor',fgd, 'ZColor',fgd)
		end
	end
%
	function fSizeObj(~,~) % Resize the figure contents.
		drawnow()
		try
			figP = get(figH, 'Position');
		catch %#ok<CTCH>
			return
		end
		% Ensure minimum virtual figure size:
		adj = max(figP(3:4),[425,254]);
		pFg = [figP(1:2)+min(0,figP(3:4)-adj),adj];
		% Get object sizes:
		mUiX = cell2mat(get(mUiH, 'Extent'));
		sTxX = cell2mat(get(sTxH, 'Extent'));
		tCbX = cell2mat(get(tCbH, 'Extent'));
		tUiX = cell2mat(get(tUiH, 'Extent'));
		tTxX = cell2mat(get(tTxH, 'Extent'));
		% Group widths and heights:
		tWd = max(tUiX(:,3))+21; % table width
		bWd = 36; % colorbar axes width
		aWd = pFg(3)-tWd-bWd-gap*4; % main axes width
		mHt = max(mUiX(:,4)); % menu height
		% Menu UI positions:
		mUiP = mUiX;
		mUiP(:,1) = aWd+bWd+3*gap;
		mUiP(:,2) = gap+2*mHt*(mTxN-1:-1:0)+3;
		mUiP(:,3) = tWd;
		mUiP(:,4) = mHt;
		set(mUiH,{'Position'},num2cell(mUiP,2))
		% Menu text positions:
		mAxP = mUiP(end,:);
		mAxP(:,4) = 2*mTxN*mHt;
		set(mAxH,'Position',mAxP)
		% Horizontal slider text positions:
		sHt = (2*mTxN*mHt)/sTxN;
		sTxP = sTxX;
		sTxP(:,1) = gap;
		sTxP(:,2) = gap+sHt*(sTxN-1:-1:0);
		sTxP(:,3) = max(sTxX(:,3));
		sTxP(:,4) = sHt;
		set(sTxH,{'Position'},num2cell(sTxP,2))
		% Horizontal slider value positions:
		sVaP = sTxP;
		sVaP(:,1) = sum(sTxP(:,[1,3]),2);
		set(sVaH,{'Position'},num2cell(sVaP,2))
		% Horizontal slider UI positions:
		sUiP = sVaP;
		sUiP(:,1) = sum(sVaP(:,[1,3]),2);
		sUiP(:,3) = aWd-sTxP(:,3)-sVaP(:,3);
		set(sUiH,{'Position'},num2cell(sUiP,2))
		% Table UI positions:
		bHt = pFg(4)-2*mTxN*mHt-3*gap;
		tUiP = tUiX;
		tUiP(:,1) = aWd+bWd+3*gap;
		tUiP(:,2) = 2*mTxN*mHt+2*gap+[bHt/2;0];
		tUiP(:,3) = tWd;
		tUiP(:,4) = bHt/2;
		set(tUiH,{'Position'},num2cell(tUiP,2))
		% Table text positions:
		tTxP = tTxX;
		tTxP(:,1) = tUiP(:,1)+3;
		tTxP(:,2) = tUiP(:,2)+bHt/2-tTxX(:,4);
		tTxP(:,4) = tTxX(:,4)-3;
		set(tTxH,{'Position'},num2cell(tTxP,2))
		% Table checkbox UI positions:
		tCbP = tCbX(:,[1,2,4,4]);
		tCbP(:,1) = tUiP(:,1)+tUiP(:,3)-21;
		tCbP(:,2) = tUiP(:,2)+tUiP(:,4)-21;
		set(tCbH,{'Position'},num2cell(tCbP,2))
		% Colorbar axes position:
		bAxP = tUiP(end,:);
		bAxP(1) = aWd+gap*2;
		bAxP(3) = bWd;
		bAxP(4) = bHt;
		set(bAxH,'Position',bAxP)
		% Number slider UI position:
		nUiP = bAxP;
		nUiP(2) = gap;
		nUiP(4) = 2*mTxN*mHt-txh-gap-mHt;
		set(nUiH,'Position',nUiP)
		% Number slider value position:
		nVaP = nUiP;
		nVaP(2) = gap+nUiP(4);
		nVaP(4) = txh;
		set(nVaH,'Position',nVaP)
		% Eval button UI position:
		eUiP = nVaP;
		eUiP(2) = 2*mTxN*mHt-mHt+gap;
		eUiP(4) = mHt;
		set(eUiH,'Position',eUiP)
		% Main plots uipanel position:
		pPnP = bAxP;
		pPnP(1) = gap;
		pPnP(3) = aWd;
		set(pPnH,'Position',pPnP)
	end
%
if nargout
	waitfor(figH)
	rgb = bgr;
end
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%maxdistcolor_view

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function R = RiemannND(x)
% Riemann's non-differentiable function, R(x) for rational x = p/q

[p,q] = rat(x);
R = 0;
for k = 1:q-1
    R = R + sin(k^2*p*pi/q)/(sin(k*pi/2/q))^2;
end
R = R*pi/4/q/q;
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function varargout = cosh(varargin)
%COSH (overloaded)

switch class(varargin{1})

    case 'double'
        error('Overloaded SDPVAR/COSH CALLED WITH DOUBLE. Report error')

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = struct('convexity','none','monotonicity','none','definiteness','positive','model','callback');
        operator.convexhull = [];
        operator.bounds = @bounds;
        operator.derivative = @(x)(sinh(x));

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error('SDPVAR/COSH called with CHAR argument?');
end

function [L,U] = bounds(xL,xU)
if xL<0 & xU>0
    L = 0;
    U = max([cosh(xL) cosh(xU)]);
elseif xL<0
    L = cosh(xU);
    U = cosh(xL);
else
    L = cosh(xL);
    U = cosh(xU);
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function test_libltfat_fifo
gl = 20;
fifoLen = 120; % Must be at least as big as gl + max expected bufLen
M = gl;
a = 7;
g = firwin('hann',gl);
gd = gabdual(g,a,M);
gg = fftshift(g.*gd)*M;

fifoPtr = calllib('libltfat','rtdgtreal_fifo_init_d',fifoLen,gl,a,1);
fifoPtr.Value.buf.setdatatype('doublePtr',fifoLen+1)
ififoPtr = calllib('libltfat','rtidgtreal_fifo_init_d',fifoLen,gl,a,1);
ififoPtr.Value.buf.setdatatype('doublePtr',fifoLen+gl+1)

bufIn = (1:1000)';
%bufIn = ones(1,1000);
bufOut = zeros(size(bufIn));
bufLen = 100;
bufOutPtr = libpointer('doublePtr',zeros(gl,1));

for ii=1:length(bufIn)/bufLen
slice = (ii-1)*bufLen + 1 : ii*bufLen;
buf = bufIn(slice);
bufInPtr = libpointer('doublePtr',buf);
bufInPtrTmp = libpointer('doublePtr',zeros(size(buf)));

written = calllib('libltfat','rtdgtreal_fifo_write_d',fifoPtr,bufLen,bufInPtr);

while calllib('libltfat','rtdgtreal_fifo_read_d',fifoPtr,bufOutPtr) > 0
    bufOutPtr.Value = bufOutPtr.Value.*gg;
    written = calllib('libltfat','rtidgtreal_fifo_write_d',ififoPtr,bufOutPtr)
end

read = calllib('libltfat','rtidgtreal_fifo_read_d',ififoPtr,bufLen,bufInPtrTmp)

bufOut(slice) = bufInPtrTmp.Value;

end

inshift = circshift(bufIn,(gl-1));
inshift(1:(gl-1)) = 0;
stem([bufOut, inshift]);shg;

calllib('libltfat','rtdgtreal_fifo_done_d',fifoPtr);
calllib('libltfat','rtidgtreal_fifo_done_d',ififoPtr);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function bsdsAnalyzeResult(pb, imfn, outdir)

for f = 1:numel(pb)
    tmp = load(['~/data/eccv08/bsds/pb/' strtok(imfn{f}, '.') '_pb.mat']);
    pb{f} = bsdsNonMaxSuppression(pb{f}, tmp.pb);
    imwrite(pb{f}, [outdir strtok(imfn{f}, '.') '.bmp']);
end

boundaryBench(outdir,'color', 25, 1)

author: Cloud-CV

repo: object-proposals


text:
function [ net,res,opts ] = qrnn_ff( net,opts )
%NET_FF Summary of this function goes here
%   Detailed explanation goes here

    n_frames=opts.parameters.n_frames;
    n_hidden_nodes=opts.parameters.n_hidden_nodes;
    n_input_nodes=size(opts.input_data,1);
    batch_size=opts.parameters.batch_size;

    %classification
    if isfield(opts,'input_labels')
        res.Fit(1).class=opts.input_labels(:);%vectorize all outputs
    end

    %Process inputs, clarity is more important here
    res.Gates(1).x=reshape(opts.input_data,n_input_nodes,[]);%inputs
    res.Input(1).x=res.Gates(1).x;%inputs

    [ net{1},res.Gates,opts ] = net_ff( net{1},res.Gates,opts );
    [ net{2},res.Input,opts ] = net_ff( net{2},res.Input,opts );

    res.Gates(end).x=reshape(res.Gates(end).x,n_hidden_nodes,batch_size,n_frames);
    res.Input(end).x=reshape(res.Input(end).x,n_hidden_nodes,batch_size,n_frames);
    res.Hidden=zeros(n_hidden_nodes,batch_size,n_frames,'like',res.Input(end).x);
    %Hidden_diff is used in bp.
    res.Hidden_diff=zeros(n_hidden_nodes,batch_size,n_frames,'like',res.Input(end).x);

    res.Hidden_diff(:,:,1)=res.Input(end).x(:,:,1);
    res.Hidden(:,:,1)=res.Gates(end).x(:,:,1).*res.Hidden_diff(:,:,1);
    for f=2:n_frames
        %Update hidden nodes;
        res.Hidden_diff(:,:,f)=res.Input(end).x(:,:,f)-res.Hidden(:,:,f-1);
        res.Hidden(:,:,f)=res.Hidden(:,:,f-1) + res.Gates(end).x(:,:,f).*res.Hidden_diff(:,:,f);
    end

    res.Fit(1).x=res.Hidden;

    %Data fitting transform
    res.Fit(1).x=reshape(res.Fit(1).x,n_hidden_nodes,[]);
    [ net{3},res.Fit,opts ] = net_ff( net{3},res.Fit,opts );

    %%%summarize
    if isfield(opts,'input_labels')
        opts.err=error_multiclass(res.Fit(1).class,res.Fit)./numel(res.Fit(1).class);
    end
    opts.loss=mean(res.Fit(end).x(:));
end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function [V] = GNMFtest(U, X)
% Graph regularized Non-negative Matrix Factorization (GNMF) handling the test data
%
% Notation:
% X ... (mFea x nSmp) data matrix
%       mFea  ... number of words (vocabulary size)
%       nSmp  ... number of documents
% U ... (mFea x K) basis matrix learned by GNMF.
%
% X = U*V'
%
% References:
% [1] Deng Cai, Xiaofei He, Xiaoyun Wu, and Jiawei Han. "Non-negative
% Matrix Factorization on Manifold", Proc. 2008 Int. Conf. on Data Mining
% (ICDM'08), Pisa, Italy, Dec. 2008.
%
% [2] Deng Cai, Xiaofei He, Jiawei Han, Thomas Huang. "Graph Regularized
% Non-negative Matrix Factorization for Data Representation", IEEE
% Transactions on Pattern Analysis and Machine Intelligence, , Vol. 33, No.
% 8, pp. 1548-1560, 2011.
%
%
%   version 1.0 --Jan./2012
%
%   Written by Deng Cai (dengcai AT gmail.com)
%

if min(min(X)) < 0
    error('Input should be nonnegative!');
end

UTU = U'*U;
UTU = max(UTU,UTU');
UTX = U'*X;
V = max(0,UTU\UTX);
V = V';

author: ZJULearning

repo: MatlabFunc


text:
function  exponent_m = separable(exponent_m,exponent_p,options);
%SEPARABLE Internal function, not used

% %exponent_m(sum((exponent_m>0),2)>2,:)=[];
%
% card = max(sum((exponent_p>0),2));
%
% n_less = exponent_m(sum((exponent_m>0),2)<card,:);
% n_equ  = exponent_m(sum((exponent_m>0),2)==card,:);
% n_larg = exponent_m(sum((exponent_m>0),2)>card,:);
%
% A = minksum(n_less,n_less);
% B = minksum(n_less,n_equ);
% C = minksum(n_less,n_larg);
% D = minksum(n_equ,n_equ);
% E = minksum(n_equ,n_larg);
% F = minksum(n_larg,n_larg);

disconnected = [];
for i = 1:size(exponent_p,2)
    for j = i+1:size(exponent_p,2)
        if ~any(exponent_p(:,i) & exponent_p(:,j))
            disconnected = [disconnected;i j];
        end
    end
end

for i = 1:size(disconnected,1)
    j = disconnected(i,1);
    k = disconnected(i,2);
    n0 = find(~exponent_m(:,j) & ~exponent_m(:,k));
    nx = find(exponent_m(:,j) & ~exponent_m(:,k));
    nz = find(~exponent_m(:,j) & exponent_m(:,k));
    nxz = find(exponent_m(:,j) & exponent_m(:,k));

%     m0 = exponent_m(n0,:);
%     mx = exponent_m(nx,:);
%     mz = exponent_m(nz,:);
%     mxz = exponent_m(nxz,:);
%
%     from_E = minksum(mx,mz);
%     from_B = minksum([m0;mx;mz],mxz);
%     from_C = minksum(mxz,mxz);
%     m_e = exponent_m(union(nx,nz),:)
%     m_cb = exponent_m(union(nx,nz),:)
    exponent_m = exponent_m([n0;nx;nz],:);
end

function msum = minksum(a,b);
msum = [];
for i = 1:size(a,1)
    for j = i:size(b,1)
        msum = [msum;a(i,:)+b(j,:)];
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function oldSettings = fillPage(h, varargin)
%FILLPAGE set figure h to fill a page when printed/exported
%   FILLPAGE will adjust printing/exporting related properties so the
%   figure (h) will take up a full page when printed/exported. Optional
%   arguments (varargin) allow you to adjust the margins and page size used.
%
%	OLDSETTINGS = FILLPAGE(fig) returns a structure containing the names of
%	page-related properties that **may** have been affected by FILLPAGE, and
%	the original values of those properties (before any FILLPAGE changes)
%   The figure settings can be restored by
%      set(fig, oldsettings);
%
%   FILLPAGE also accepts up to 2 optional parameter-value pairs:
%       parameter      value/meaning
%       ----------     --------------------------------------
%        'margins'     1x4 numeric array specifying the
%                      margins to use on the left, right, top,
%                      and bottom, respectively. Units used are the
%                      figure's current PaperUnits.
%
%                      If not specified, a margin of [.25 .25 .5 .5]
%                      is used.
%
%        'papersize'   a 1x2 numeric array specifying the
%                      width and height of the page, respectively.
%                      Units used are the figure's current PaperUnits.
%                                  OR
%                      a string specifying the paper type to use
%                        (see the figure's PaperType property for
%                        valid values)
%
%                      If not specified, the figure's current PaperSize
%                      is used.
%
%   Example usage:
%     % fills page based on current settings
%        figure;
%        surf(peaks);
%        oldSettings = fillpage(gcf);
%     % fills page leaving a margin of .5 around all sides
%        figure;
%        surf(peaks);
%        oldSettings = fillpage(gcf, 'margins', [.5 .5 .5 .5]);
%     % fills page using A4 paper size
%        figure;
%        surf(peaks);
%        oldSettings = fillpage(gcf, 'papersize', 'A4');
%     % fills page using a paper size of 4x6
%        figure;
%        surf(peaks);
%        oldSettings = fillpage(gcf, 'papersize', [4 6]);
%     % fills page using USLetter paper size with margin of 1 all around
%        figure;
%        surf(peaks);
%        oldSettings = fillpage(gcf, 'papersize', 'usletter', 'margins', [1 1 1 1]);
%
% See also PRINT

%   Copyright 2007 The MathWorks, Inc.

% set up some meaningful aliases for index values
IDX_PAPERTYPE = 1;
IDX_PAPERSIZE = 2;
IDX_PAPERPOSITION = 3;
IDX_PAPERPOSITIONMODE = 4;

if ishandle(h)
    setprops = {'PaperType', 'PaperSize', 'PaperPosition', 'PaperPositionMode'};
    try
        settings = get(h, setprops);
    catch
        error('fillPage unable to get figure properties');
    end

    oldSettings = cell2struct(settings, setprops, 2); % to return to caller
    paramSettings = processArgs(varargin{:});
    if isstruct(paramSettings)
        % process papersize/type
        if isfield(paramSettings, 'papertype')
            if ~isempty(paramSettings.papertype)
                settings{IDX_PAPERTYPE} = paramSettings.papertype;
                % surround with try/catch in case setting is invalid
                try
                   set(h, 'PaperType', paramSettings.papertype);
                catch
                    error('fillPage: invalid PaperType ''%s'' requested ', ...
                        paramSettings.papertype );
                end
                if isfield(paramSettings, 'papersize') && ...
                    ~isempty(paramSettings.papersize)
                    settings{IDX_PAPERSIZE} = paramSettings.papersize;
                else
                    % since changing papertype might change papersize...
                    % re-read value
                    settings{IDX_PAPERSIZE} = get(h, 'PaperSize');
                end
            end
        end

        % process margins
        if isfield(paramSettings, 'margins')
            margins = paramSettings.margins;
            papersize = settings{IDX_PAPERSIZE};
            settings{IDX_PAPERPOSITIONMODE} = 'manual';
            settings{IDX_PAPERPOSITION} = [ ...
                 margins.left ... % X = leftMargin
                 margins.bottom ... % Y = bottomMargin
                 papersize(1) - (margins.right + margins.left) ... % W =  paperwidth - rightMargin - leftMargin
                 papersize(2) - (margins.top + margins.bottom) ... % H = paperheight - topMargin - bottomMargin
                 ];
        end
    end
    set(h, setprops, settings);
else
    error('fillPage requires a handle to a figure');
end
end

% helper function to deal with parsing command line arguments
function param_settings = processArgs(varargin)
    margins.left   = .25;
    margins.right  = .25;
    margins.top    = .5;
    margins.bottom = .5;
    papertype = [];
    papersize = [];

    for i = 1 : 2 : length(varargin)-1
        param_arg = varargin{i};

        % if caller specified margins to use
        if ischar(param_arg) && strcmpi(param_arg, 'margins')
           val_arg = varargin{i+1};
           if isnumeric(val_arg) && length(val_arg) == 4
               margins.left   = val_arg(1);
               margins.right  = val_arg(2);
               margins.top    = val_arg(3);
               margins.bottom = val_arg(4);
           else
               warning('fillpage:InvalidMargin', ...
                   'fillpage ignoring invalid margin setting; using default');
           end
        end

        % if caller specified papersize to use
        if ischar(param_arg) && strcmpi(param_arg, 'papersize')
           val_arg = varargin{i+1};
           if ischar(val_arg)
              papertype = val_arg;
            elseif isnumeric(val_arg) && length(val_arg) == 2
              papertype = '<custom>';
              papersize = val_arg;
           else
               warning('fillpage:InvalidPapersize', ...
                   'fillpage ignoring invalid papersize setting; using figure''s current papersize');
           end
        end
    end

    param_settings.margins = margins;
    param_settings.papertype = papertype;
    param_settings.papersize = papersize;
end

author: Sable

repo: mcbench-benchmarks


text:
function meshfaces_test ( )

%*****************************************************************************80
%
%% MAIN tests the MESHFACES library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 August 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MESHFACES_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the MESHFACES library.\n' );

  meshfaces_test01 ( );
  meshfaces_test02 ( );
  meshfaces_test03 ( );
  meshfaces_test04 ( );
  meshfaces_test05 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MESHFACES_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clear all; clc;

newlineInAscii1 = [13 10];
spaceInInAscii = 32;
% for printing, newline causes much confusion in matlab and is provided here as an alternative
newline = char(newlineInAscii1);
spaceChar = char(spaceInInAscii);

targetFolder = 'images';


text:
function option_num = p08_option_num ( )

%*****************************************************************************80
%
%% P08_OPTION_NUM returns the number of options for problem 8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer OPTION_NUM, the number of options.
%
  option_num = 1;

  return
end


text:
classdef LIRCMOP6 < PROBLEM
% <multi> <real> <large/none> <constrained>
% Constrained benchmark MOP with large infeasible regions

%------------------------------- Reference --------------------------------
% Z. Fan, W. Li, X. Cai, H. Huang, Y. Fang, Y. You, J. Mo, C. Wei, and E.
% Goodman, An improved epsilon constraint-handling method in MOEA/D for
% CMOPs with large infeasible regions, Soft Computing, 2019, 23:
% 12491-12510.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Wenji Li

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values and constraint violations
        function Population = Evaluation(obj,varargin)
            X = varargin{1};
            X = max(min(X,repmat(obj.upper,size(X,1),1)),repmat(obj.lower,size(X,1),1));
            [popsize,variable_length] = size(X);
            sum1 = zeros(popsize,1);
            sum2 = zeros(popsize,1);
            for j = 2 : variable_length
                if mod(j,2) == 1
                    sum1 = sum1+(X(:,j)-sin((0.5*j/variable_length*pi)*X(:,1))).^2;
                else
                    sum2 = sum2+(X(:,j)-cos((0.5*j/variable_length*pi)*X(:,1))).^2;
                end
            end
            gx          = 0.7057;
            PopObj(:,1) = X(:,1)+10*sum1+gx;
            PopObj(:,2) = 1-X(:,1).^2+10.*sum2+gx;
            Population  = SOLUTION(X,PopObj,Constraint(PopObj),varargin{2:end});
            obj.FE      = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^2;
            R      = R + 0.7057;
            R(any(Constraint(R)>0,2),:) = [];
        end
        %% Generate the feasible region
        function R = GetPF(obj)
            [x,y] = meshgrid(linspace(0.7057,5,400));
            fes   = all(Constraint([x(:),y(:)])<=0,2);
            z     = nan(size(x));
            z(reshape(fes,size(z)) & (x-0.7057).^2+y>=1.7057) = 0;
            R = {x,y,z};
        end
    end
end

function PopCon = Constraint(PopObj)
    p     = [1.8,2.8];
    q     = [1.8,2.8];
    a     = [2,2];
    b     = [8,8];
    r     = 0.1;
    theta = -0.25 * pi;
    for k = 1 : 2
        PopCon(:,k) = r - ((PopObj(:,1)-p(k))*cos(theta)-(PopObj(:,2)-q(k))*sin(theta)).^2/(a(k)^2) -...
                      ((PopObj(:,1)-p(k))*sin(theta)+(PopObj(:,2)-q(k))*cos(theta)).^2/(b(k)^2);
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function val = dtiGetImageCoordsFromFibers(dt, fg)
%Obsolete
%
%  Returns the image coordinates of the fibers in fg
%
% WARNING:  Doesn't work.  See notes below.  Should be obsolete soon.
%
%  val = dtiGetImageCoordsFromFibers(dt, fg)
%
% The data in val are Nx3
%
% (c) Stanford VISTA Team

% Written by Sherbondy and Dougherty.  It doesn't seem to do what it says.
% It is not called by any other functions. We need a function that does
% just this properly.  It should be inside of dtiGet(), I think.

error('Obsolete %s\n',mfilename)

return

% Here are the coords
coords = horzcat(fg.fibers{:})';
coords = unique(round(mrAnatXformCoords(inv(dt.xformToAcpc),coords)),'rows');

% Here are the tensors as 6-vectors
inds = sub2ind(size(dt.dt6(:,:,:,1)),coords(:,1),coords(:,2),coords(:,3));
vecDt6 = zeros(size(coords,1),6);
for ii=1:6
    temp = dt.dt6(:,:,:,ii);
    vecDt6(:,ii) = temp(inds);
end

% These statistics should be computed based on nargout.
[eigVec, eigVal] = dtiEig(vecDt6);
l2 = eigVal
[fa,md,rd] = dtiComputeFA(eigVal);
[cl, cp, cs] = dtiComputeWestinShapes(eigVal);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_bug2440

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_freqsimulation ft_freqanalysis ft_singleplotTFR

cfg = [];
cfg.method = 'broadband';
cfg.numtrl = 5;
cfg.trllen = 2;

data = ft_freqsimulation(cfg);

cfg = [];
cfg.method = 'mtmconvol';
cfg.taper = 'hanning';
cfg.foi = 8:12;
cfg.t_ftimwin = 4./cfg.foi;
cfg.toi = 0:0.1:2;
cfg.keeptrials = 'yes';
freq = ft_freqanalysis(cfg, data);

freq.logspctrm = log(freq.powspctrm);

% plot
cfg = [];
cfg.parameter = 'logspctrm';
cfg.trials = 2;
ft_singleplotTFR(cfg,freq);

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [select] = parameterselection(param, data)

% PARAMETERSELECTION selects the parameters that are present as a volume in the data
% add that have a dimension that is compatible with the specified dimensions of the
% volume, i.e. either as a vector or as a 3D volume.
%
% Use as
%   [select] = parameterselection(param, data)
% where
%   param    cell-array, or single string, can be 'all'
%   data     structure with anatomical or functional data
%   select   returns the selected parameters as a cell-array

% Copyright (C) 2005-2008, 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 ischar(param)
  param = {param};   % it should be a cell-array
elseif isempty(param)
  param = {};        % even being empty, it should be a cell-array
end

sel = find(strcmp(param, 'all'));
if ~isempty(sel)
  % the old default was a list of all known volume parameters
  % the new default is to try all fields present in the data
  allparam = fieldnames(data);

  % Around 1 May 2021, seg has been renamed in tissue in the automatic
  % conversion of probablistic to indexed segmentations/parcellations.
  if ismember('seg', param) && ~ismember('seg', allparam) && ismember('tissue', allparam)
    ft_warning('the field ''seg'' has been renamed to ''trial'', please update your code')
    param{strcmp(param, 'seg')} = 'tissue';
  end

  % fields can be nested in source.avg
  if isfield(data, 'avg') && isstruct(data.avg)
    tmp = fieldnames(data.avg);
    for i=1:length(tmp)
      tmp{i} = ['avg.' tmp{i}];
    end
    allparam = cat(1, allparam, tmp);
  end

  % fields can be nested in source.trial
  if isfield(data, 'trial') && isstruct(data.trial)
    tmp = fieldnames(data.trial);
    for i=1:length(tmp)
      tmp{i} = ['trial.' tmp{i}];
    end
    allparam = cat(1, allparam, tmp);
  end
  param(sel) = [];                          % remove the 'all'
  param      = [param(:)' allparam(:)'];    % add the list of all possible parameters, these will be tested later

else
  % check all specified parameters and give support for some parameters like 'pow' and 'coh'
  % which most often will indicate 'avg.pow' and 'avg.coh'
  for i=1:length(param)
    if ~issubfield(data, param{i}) && issubfield(data, ['avg.' param{i}])
      % replace the parameter xxx by avg.xxx
      param{i} = ['avg.' param{i}];
    end
  end
end

% remove empty fields
param(cellfun('isempty', param)) = [];

% ensure that there are no double entries
param = unique(param);

select = {};
for i=1:length(param)
  if issubfield(data, param{i})
    % the field is present, check whether the dimension is correct
    dim = size(getsubfield(data, param{i}));
    if isfield(data, 'dim') && isequal(dim(:), data.dim(:))
      select{end+1} = param{i};
    elseif isfield(data, 'dim') && prod(dim)==prod(data.dim)
      select{end+1} = param{i};
    elseif isfield(data, 'dim') && numel(dim)==3 && isequal(dim(1:3)', data.dim(:))
      select{end+1} = param{i};
    elseif isfield(data, 'pos') && (prod(dim)==size(data.pos, 1) || dim(1)==size(data.pos,1))
      select{end+1} = param{i};
    elseif isfield(data, 'dimord') && (isfield(data, 'pos') || isfield(data, 'transform'))
      dimtok = tokenize(data.dimord, '_');
      nels   = 1;
      for k=1:numel(dimtok)
        if strcmp(dimtok{k}, 'rpt') || strcmp(dimtok{k}, 'rpttap')
          nels = nels*dim(k);
        elseif strcmp(dimtok{k}, 'pos') && isfield(data, 'pos')
          nels = nels*size(data.pos,1);
        elseif strcmp(dimtok{k}, '{pos}') && isfield(data, 'pos')
          nels = nels*size(data.pos,1);
        elseif isfield(data, dimtok{k})
          nels = nels*numel(getfield(data, dimtok{k}));
        end
      end
      if nels==prod(dim)
        select{end+1} = param{i};
      end
    elseif isfield(data, 'dim') && numel(dim)>3 && isequal(dim(1:3), data.dim(1:3))
      select{end+1} = param{i};
    end
  end
end


text:
function a = band(a,p,q)
%BAND         Extract band from matrix a, lower bandwidth p, upper bandwidth q
%   if parameter q is omitted, q:=p
%
%   c = band(a,p,q)
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if nargin<3
    q = p;
  end

  index = ( tril(triu(ones(size(a.x)),-p),q) == 0 );
  a.x(index) = 0;
  a.dx(:,index) = 0;
  a.hx(:,index) = 0;

repo: OpenMAS


text:
function tf = isTestCaseSubclass(name)
%isTestCaseSubclass True for name of a TestCase subclass
%   tf = isTestCaseSubclass(name) returns true if the string name is the name of
%   a TestCase subclass on the MATLAB path.

%   Steven L. Eddins
%   Copyright 2008-2009 The MathWorks, Inc.

tf = false;

class_meta = meta.class.fromName(name);
if isempty(class_meta)
    % Not the name of a class
    return;
end

if strcmp(class_meta.Name, 'TestCase')
    tf = true;
else
    tf = isMetaTestCaseSubclass(class_meta);
end

function tf = isMetaTestCaseSubclass(class_meta)

tf = false;

if strcmp(class_meta.Name, 'TestCase')
    tf = true;
else
    % Invoke function recursively on parent classes.
    super_classes = class_meta.SuperClasses;
    for k = 1:numel(super_classes)
        if isMetaTestCaseSubclass(super_classes{k})
            tf = true;
            break;
        end
    end
end

author: ricedsp


text:
function [xvec, yvec]=generatepointsregular(npoints,diamcirc,origin,phaseoffset)
% [xvec, yvec]=generatepointsregular(npoints,diamcirc,origin,pahseoffset)
% generates points qually distributed around the circle with 2r=diamcirc
% with center at origin
% phaseoffset: offster of the phase of the first point (default == 0);
if ~exist('phaseoffset','var') %to make oblique configuration
    phaseoffset = 0;
end
angstep=2*pi/npoints;
angvectmp=(0:angstep:2*pi)+phaseoffset;
angvec=angvectmp(1:end-1); %last opint = first point
xvec=(origin(1)+0.5*diamcirc*cos(angvec))';
yvec=(origin(2)+0.5*diamcirc*sin(angvec))';

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [P, p] = prob_node(CPD, self_ev, pev)
% PROB_NODE Compute prod_m P(x(i,m)| x(pi_i,m), theta_i) for node i (discrete)
% [P, p] = prob_node(CPD, self_ev, pev)
%
% self_ev(m) is the evidence on this node in case m.
% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents).
% (These may also be cell arrays.)
%
% p(m) = P(x(i,m)| x(pi_i,m), theta_i)
% P = prod p(m)

if iscell(self_ev), usecell = 1; else usecell = 0; end

ncases = length(self_ev);
sz = dom_sizes(CPD);

nparents = length(sz)-1;
if nparents == 0
  assert(isempty(pev));
else
  assert(isequal(size(pev), [nparents ncases]));
end

n = length(sz);
dom = 1:n;
p = zeros(1, ncases);
if nparents == 0
  for m=1:ncases
    if usecell
      evidence = {self_ev{m}};
    else
      evidence = num2cell(self_ev(m));
    end
    T = convert_to_table(CPD, dom, evidence);
    p(m) = T;
  end
else
  for m=1:ncases
    if usecell
      evidence = cell(1,n);
      evidence(1:n-1) = pev(:,m);
      evidence(n) = self_ev(m);
    else
      evidence = num2cell([pev(:,m)', self_ev(m)]);
    end
    T = convert_to_table(CPD, dom, evidence);
    p(m) = T;
  end
end
P = prod(p);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function p = predict(Theta1, Theta2, X)
%PREDICT Predict the label of an input given a trained neural network
%   p = PREDICT(Theta1, Theta2, X) outputs the predicted label of X given the
%   trained weights of a neural network (Theta1, Theta2)

% Useful values
m = size(X, 1);
num_labels = size(Theta2, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned neural network. You should set p to a
%               vector containing labels between 1 to num_labels.
%
% Hint: The max function might come in useful. In particular, the max
%       function can also return the index of the max element, for more
%       information see 'help max'. If your examples are in rows, then, you
%       can use max(A, [], 2) to obtain the max for each row.
%

a1 = [ones(m, 1) X]; % add 1

z2 = a1 * Theta1';

a2 = [ones(size(sigmoid(z2), 1), 1) sigmoid(z2)]; % add 1

z3 = a2 * Theta2';
a3 = sigmoid(z3); % H_theta(x)

[x, ix] = max(a3, [], 2);

p = ix;

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

end

repo: machine-learning-coursera


text:
function hb_io_test11 ( )

%*****************************************************************************80
%
%% HB_IO_TEST11 tests HB_MATVEC_A_MEM;
%
%  Modified:
%
%    04 February 2005
%
  ncol = 32;
  neltvl = 0;
  nnzero = 126;
  nrhs = 2;
  nrhsix = 0;
  nrow = 32;

  colptr = [ ...
      1,   7,  12,  18,  22,  26,  29,  34,  39,  46, ...
     53,  58,  61,  63,  65,  68,  71,  74,  79,  82, ...
     85,  88,  90,  94,  97, 102, 106, 110, 112, 117, ...
    121, 124, 127 ]';
  exact = [ ...
    0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   1.0, ...
    0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0, ...
    0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0, ...
    0.0,   0.0;
    1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0, ...
    1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0, ...
    1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0, ...
    1.0,   1.0 ]';
  guess = [ ...
    1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0, ...
    1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0, ...
    1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0, ...
    1.0,   1.0;
    1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0, ...
    1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0, ...
    1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0, ...
    1.0,   1.0 ]';
  indcrd = 8;
  indfmt = '(16I5)';
  key = 'RUA_32';
  mxtype = 'RUA';
  output_file = 'rua_32_ax.txt';
  ptrcrd = 3;
  ptrfmt = '(16I5)';
  rhscrd = 12;
  rhsfmt = '(10F7.1)';
  rhsind = 0;
  rhsptr = 0;
  rhstyp = 'FGX';
  rhsvec = 0;
  rowind = [ ...
    1,    2,    3,    4,    7,   26,    1,    2,    9,   21, ...
   28,    2,    3,    6,    8,    9,   29,    3,    4,    5, ...
   12,    3,    5,   23,   27,    1,    6,   16,    3,    7, ...
   14,   21,   31,    1,    8,   12,   17,   27,    7,    9, ...
   10,   13,   19,   23,   27,    1,   10,   11,   21,   23, ...
   25,   27,    2,   11,   15,   18,   29,    6,   12,   24, ...
   11,   13,    3,   14,    2,   15,   20,    4,   16,   22, ...
    4,   16,   17,    6,   10,   18,   20,   30,    1,   19, ...
   26,    8,   16,   20,    3,   21,   32,   11,   22,    2, ...
   17,   21,   23,   12,   24,   26,    6,   15,   18,   24, ...
   25,   13,   18,   22,   26,    5,   24,   26,   27,    9, ...
   28,    3,    5,   27,   29,   32,   12,   17,   23,   30, ...
   13,   14,   31,   24,   28,   32 ]';
  title = '1Real unsymmetric assembled matrix based on IBM32';
  totcrd = 36;
  valcrd = 13;
  valfmt = '(10F7.1)';
  values = [ ...
  101.0,  102.0,  103.0,  104.0,  107.0, ...
  126.0,  201.0,  202.0,  209.0,  221.0, ...
  228.0,  302.0,  303.0,  306.0,  308.0, ...
  309.0,  329.0,  403.0,  404.0,  405.0, ...
  412.0,  503.0,  505.0,  523.0,  527.0, ...
  601.0,  606.0,  616.0,  703.0,  707.0, ...
  714.0,  721.0,  731.0,  801.0,  808.0, ...
  812.0,  817.0,  827.0,  907.0,  909.0, ...
  910.0,  913.0,  919.0,  923.0,  927.0, ...
 1001.0, 1010.0, 1011.0, 1021.0, 1023.0, ...
 1025.0, 1027.0, 1102.0, 1111.0, 1115.0, ...
 1118.0, 1129.0, 1206.0, 1212.0, 1224.0, ...
 1311.0, 1313.0, 1403.0, 1414.0, 1502.0, ...
 1515.0, 1520.0, 1604.0, 1616.0, 1622.0, ...
 1704.0, 1716.0, 1717.0, 1806.0, 1810.0, ...
 1818.0, 1820.0, 1830.0, 1901.0, 1919.0, ...
 1926.0, 2008.0, 2016.0, 2020.0, 2103.0, ...
 2121.0, 2132.0, 2211.0, 2222.0, 2302.0, ...
 2317.0, 2321.0, 2323.0, 2412.0, 2424.0, ...
 2426.0, 2506.0, 2515.0, 2518.0, 2524.0, ...
 2525.0, 2613.0, 2618.0, 2622.0, 2626.0, ...
 2705.0, 2724.0, 2726.0, 2727.0, 2809.0, ...
 2828.0, 2903.0, 2905.0, 2927.0, 2929.0, ...
 2932.0, 3012.0, 3017.0, 3023.0, 3030.0, ...
 3113.0, 3114.0, 3131.0, 3224.0, 3228.0, ...
 3232.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HB_IO_TEST11\n' );
  fprintf ( 1, '  HB_MATVEC_A_MEM multiplies a matrix times a vector.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This particular version assumes:\n' );
  fprintf ( 1, '  * the matrix is in "A" format (assembled),\n' );
  fprintf ( 1, '  * the matrix and vectors can fit in memory,\n' );
  fprintf ( 1, '  * the matrix and multiplicand have been read into\n' );
  fprintf ( 1, '    memory before the routine is called.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  For this example, the first vector X is zero except\n' );
  fprintf ( 1, '  for a 1 in row 10.  This means A*X should return\n' );
  fprintf ( 1, '  column 10 of A.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The second vector X is all 1''s.  A*X should be\n' );
  fprintf ( 1, '  the sum of the entries of each row.\n' );

  rhsval = hb_matvec_a_mem ( nrow, ncol, nnzero, nrhs, colptr, rowind, ...
    values, exact );

  r8mat_print ( nrow, nrhs, rhsval,  '  The product vectors A*X' );

  output_unit = fopen ( output_file, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST11 - Fatal error!\n' );
    fprintf ( 1, '  Error opening the file.\n' );
    return;
  end

  hb_file_write ( output_unit, title, key, totcrd, ptrcrd, indcrd, ...
    valcrd, rhscrd, mxtype, nrow, ncol, nnzero, neltvl, ptrfmt, indfmt, ...
    valfmt, rhsfmt, rhstyp, nrhs, nrhsix, colptr, rowind, values, ...
    rhsval, rhsptr, rhsind, rhsvec, guess, exact );

  fclose ( output_unit );

  return
end


text:
function varargout = subsref(f,s)
%overloads subsref

switch s(1).type
  case '()'

    f.h = subsref(f.h,s(1));
    f.o1 = subsref(f.o1,s(1));
    f.o2 = subsref(f.o2,s(1));

    if numel(s)>1
      [varargout{1:nargout}] = builtin('subsref',f,s(2:end));
    else
      varargout{1} = f;
    end
  otherwise

    [varargout{1:nargout}] = builtin('subsref',f,s);

end
end

repo: mtex


text:
function value = get_field_value(structure, name, default)

if ~isstruct(structure) || ~isfield(structure, name)
    value = default;
    return;
end;

value = structure.(name);

author: votchallenge

repo: toolkit-legacy


text:
function r = msmoothboxrnd(a,b,sigma,n)
%MSMOOTHBOXRND Random arrays from the multivariate smooth-box distribution.
%   R = MSMOOTHBOXRND(A,B,SIGMA) returns an N-by-D matrix R of random
%   vectors chosen from the multivariate smooth-box distribution
%   with pivots A and B and scale SIGMA. A, B and SIGMA are N-by-D matrices,
%   and MSMOOTHBOXRND generates each row of R using the corresponding row
%   of A, B and SIGMA.
%
%   R = MSMOOTHBOXRND(A,B,SIGMA,N) returns a N-by-D matrix R of random
%   vectors chosen from the multivariate smooth-box distribution
%   with pivots A and B and scale SIGMA.
%
%   See also MSMOOTHBOXPDF.

% Luigi Acerbi 2022

[Na,Da] = size(a);
[Nb,Db] = size(b);
[Nsigma,Dsigma] = size(sigma);

if any(sigma(:) <= 0)
    error('msmoothboxrnd:NonPositiveSigma', ...
        'All elements of SIGMA should be positive.');
end

if nargin < 4 || isempty(n)
    n = max([Na,Nb,Nsigma]);
else
    if (Na ~= 1 && Na ~= n) || (Nb ~= 1 && Nb ~= n) || ...
            (Nsigma ~= 1 && Nsigma ~= n)
        error('msmoothboxrnd:SizeError', ...
            'A, B, SIGMA should be 1-by-D or N-by-D arrays.');
    end
end
if Na ~= Nb || Da ~= Db || Na ~= Nsigma || Da ~= Dsigma
    error('msmoothboxrnd:SizeError', ...
        'A, B, SIGMA should be arrays of the same size.');
end

D = Da;

if size(a,1) == 1; a = repmat(a,[n,1]); end
if size(b,1) == 1; b = repmat(b,[n,1]); end
if size(sigma,1) == 1; sigma = repmat(sigma,[n,1]); end

r = zeros(n,D);

nf = 1 + 1/sqrt(2*pi)./sigma.*(b - a);

% Sample one dimension at a time
for d = 1:D
    % Draw component (left/right tails or plateau)
    u = nf(:,d) .* rand(n,1);

    % Left Gaussian tails
    idx = u < 0.5;
    if any(idx)
        z1 = abs(randn(sum(idx),1).*sigma(idx,d));
        r(idx,d) = a(idx) - z1;
    end

    % Right Gaussian tails
    idx = (u >= 0.5 & u < 1);
    if any(idx)
        z1 = abs(randn(sum(idx),1).*sigma(idx,d));
        r(idx,d) = b(idx) + z1;
    end

    % Plateau
    idx = u >= 1;
    if any(idx)
        r(idx,d) = a(idx,d) + (b(idx,d) - a(idx,d)).*rand(sum(idx),1);
    end
end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function bpsksys(bin,f)

disp('========================================');
disp(' HAM DIEU CHE DICH 2 PHA: BPSK');
disp(' VI DU:bpsksys([0 1 0 1 1 0 1 1 0],3)');
disp('Written by Nguyen Hoang Minh DHCNTPHCM. he..he..');
disp('========================================');

bin=[0 1 0 1 1 0 1 1 1 0];f=3;k=1000;

t=0:2*pi/(k-1):2*pi;

L = length(bin);sig0=cos(f*t);sig1=cos(f*t+pi);
bit1=ones(1,k);bit0=zeros(1,k);
mbit=[];mcw=[];

for n=1:L;
    if bin(n)==0;
       cw=sig0;bit=bit0;
    else
       cw=sig1;bit=bit1;
    end
   mbit=[mbit bit];
   mcw=[mcw cw];
end
psk=mcw;
%================================================
s=length(mcw);mrec=[];

for m=0:k:s
    if mcw(m)<0
       rec=bit1;
    elseif mcw(m)>0
       rec=bit0;
    end
    mrec=[mrec rec];
end
depsk=mrec;
subplot(3,1,1);plot(mbit,'r','linewidth',2);axis([0  k*L -0.5 1.5]);grid on;title('Data in');
subplot(3,1,2);plot(psk,'m','linewidth',1.5);axis([0  k*L -1.5 1.5]);grid on;title('PSK modulation');
subplot(3,1,3);plot(depsk,'g','linewidth',2);axis([0  k*L -.5 1.5]);grid on;title('PSK demodulation,Data out');


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Example code for estimating the HRF using the Inverse-Logit Model, a
% Finte Impluse Response Model and the Canonical HRF with 2 derivatives.
% Also the code illustrates our code for detecting model misspecification.
%
% This example illustrates the multi-stimulus version
%
% By Martin Lindquist
% Created: 05/26/10
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Load time course
%

mypath = which('ilogit');
if isempty(mypath), error('Cannot find directory with ilogit.m and other functions. Not on path?'); end
[mydir] = fileparts(mypath)

h = spm_hrf(1);
h = h./max(h);

RunA = zeros(60,10);
RunA(1,:) = 1;
RunA = reshape(RunA,600,1);
RunB = zeros(60,10);
RunB(31,:) = 1;
RunB = reshape(RunB,600,1);

tc = 2*conv(RunA,h) + conv(RunB,h);
tc = tc(1:600);
tc = tc+normrnd(0,1,600,1);

Run =[];
Run{1} = RunA;
Run{2} = RunB;

len = length(tc);

figure; subplot(3,1,1); han = plot(tc);
title('Sample time course'); drawnow

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Settings
%

TR = 1;
T = round(30/TR);
t = 1:T;                        % samples at which to get Logit HRF Estimate
FWHM = 4;                       % FWHM for residual scan
pval = 0.01;
df = 600;
alpha = 0.001;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Create stick function (sample event onsets)
% Variable R contains onset times
% Variable Run contains stick (a.k.a. delta or indicator) function

RA = [1 61 121 181 241 301 361 421 481 541];
RB = [31 91 151 211 271 331 391 451 511 571];

try
    hold on;
    hh = plot_onsets(RA,'k',-3,1);
    hh = plot_onsets(RB,'r',-3,0.5);
    drawnow
catch
    disp('Couldn''t find function to add onset sticks to plot. Skipping.')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fit HRF using IL-function

% Choose mode (deterministic/stochastic)

mode = 0;   % 0 - deterministic aproach
            % 1 - simulated annealing approach
            % Please note that when using simulated annealing approach you
            % may need to perform some tuning before use.

[h1, fit1, e1, param] = Fit_Logit(tc,Run,t,mode);
[pv sres sres_ns1] = ResidScan(e1, FWHM);
[PowLoss1] = PowerLoss(e1, fit1, (len-7) , tc, TR, Run, alpha);

hold on; han(2) = plot(fit1,'r');

disp('Summary: IL_function');

disp('HRF - Event A');
disp('Amplitude:'); disp(param(1,1));
disp('Time-to-peak:'); disp(param(2,1));
disp('Width:'); disp(param(3,1));

disp('HRF - Event B');
disp('Amplitude:'); disp(param(1,2));
disp('Time-to-peak:'); disp(param(2,2));
disp('Width:'); disp(param(3,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fit HRF using FIR-model

% Choose mode (FIR/sFIR)

mode = 1;   % 0 - FIR
            % 1 - smooth FIR

[h2, fit2, e2, param] = Fit_sFIR(tc,TR,Run,T,mode);
[pv sres sres_ns2] = ResidScan(e2, FWHM);
[PowLoss2] = PowerLoss(e2, fit2, (len-T) , tc, TR, Run, alpha);

hold on; han(3) = plot(fit2,'g');

disp('Summary: FIR');

disp('HRF - Event A');
disp('Amplitude:'); disp(param(1,1));
disp('Time-to-peak:'); disp(param(2,1));
disp('Width:'); disp(param(3,1));

disp('HRF - Event B');
disp('Amplitude:'); disp(param(1,2));
disp('Time-to-peak:'); disp(param(2,2));
disp('Width:'); disp(param(3,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fit HRF using Canonical HRF + 2 derivatives

p=1;

[h3, fit3, e3, param, info] = Fit_Canonical_HRF(tc,TR,Run,T,p);
[pv sres sres_ns3] = ResidScan(e3, FWHM);
[PowLoss3] = PowerLoss(e3, fit3, (len-p) , tc, TR, Run, alpha);

hold on; han(4) = plot(fit3,'m');

legend(han,{'Data' 'IL' 'sFIR' 'DD'})

disp('Summary: Canonical + 2 derivatives');

disp('HRF - Event A');
disp('Amplitude:'); disp(param(1,1));
disp('Time-to-peak:'); disp(param(2,1));
disp('Width:'); disp(param(3,1));

disp('HRF - Event B');
disp('Amplitude:'); disp(param(1,2));
disp('Time-to-peak:'); disp(param(2,2));
disp('Width:'); disp(param(3,2));

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

%figure;

subplot(3,2,5);
han2 = plot(h1(:,1),'r');
hold on; han2(2) = plot(h2(:,1),'g');
hold on; han2(3) = plot(h3(:,1),'m');
legend(han2,{'IL' 'sFIR' 'DD'})
title('Estimated HRF - Event A');

subplot(3,1,2); hold on;
hh = plot_onsets(RA,'k',-3,1);
hh = plot_onsets(RB,'r',-3,0.5);
drawnow

han3 = plot(sres_ns1,'r');
hold on; han3(2) = plot(sres_ns2,'g');
hold on; han3(3) = plot(sres_ns3,'m');
hold on; plot((1:len),zeros(len,1),'--k');
legend(han3,{'IL' 'sFIR' 'DD'})
title('Mis-modeling (time course)');

subplot(3,2,6); hold on;

han4 = plot(h2(:,1),'r');
hold on; han4(2) = plot(h2(:,2),'g');
hold on; han4(3) = plot(h3(:,2),'m');
legend(han4,{'IL' 'sFIR' 'DD'})
title('Estimated HRF - Event A');

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% helpforexe() - Write help files for exe version
%
% Usage:
%   histtoexe( mfile, folder)
%
% Inputs:
%   mfile  - [cell of string] Matlab files with help message
%   folder - [string] Output folder
%
% Output:
%   text files name help_"mfile".m
%
% Author: Arnaud Delorme, 2006
%
% See also: eeglab()

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

function helpforexe( funct, fold );

if nargin <1
	help histforexe;
	return;
end;
nonmatlab = 0;

% write all files
% ---------------
for index = 1:length(funct)
    doc1 = readfunc(funct{index}, nonmatlab);

    fid = fopen( fullfile(fold, [ 'help_' funct{index} ]), 'w');
    for ind2 = 1:length(doc1)
        if isempty(doc1{ind2}) fprintf(fid, [ 'disp('' '');\n' ]);
        else                   fprintf(fid, [ 'disp(' vararg2str({ doc1{ind2} }) ');\n' ]);
        end;
    end;
    fclose(fid);
    %fprintf(fid, 'for tmpind = 1:length(tmptxt), if isempty(tmptxt{tmpind}), disp('' ''); else disp(tmptxt{tmpind}); end; end; clear tmpind tmptxt;\n' );
end;

% try all functions
% -----------------
tmppath = pwd;
cd(fold);
for index = 1:length(funct)
    evalc([ 'help_' funct{index}(1:end-2) ]);
end;
cd(tmppath);
return;

%-------------------------------------
function [doc] = readfunc(funct, nonmatlab)

doc = {};
if nonmatlab
	fid = fopen( funct, 'r');
else
	if findstr( funct, '.m')
		fid = fopen( funct, 'r');
	else
		fid = fopen( [funct '.m'], 'r');
	end;
end;

if fid == -1
	error('File not found');
end;

sub = 1;
try,
    if ~isunix, sub = 0; end;
catch, end;

if nonmatlab
	str = fgets( fid );
	while ~feof(fid)
		str = deblank(str(1:end-sub));
        doc = { doc{:} str(1:end) };
        str = fgets( fid );
	end;
else
	str = fgets( fid );
	while (str(1) == '%')
		str = deblank(str(1:end-sub));
        doc = { doc{:} str(2:end) };
		str = fgets( fid );
	end;
end;
fclose(fid);

author: PatternRecognition

repo: OpenBMI


text:
% marshalling of parameters and available warping functions
function varargout = warpLog(y,hyp,varargin)
  varargout = cell(nargout, 1);  % allocate the right number of output arguments

  % return number of parameters
  m = 2;
  if nargin<1 && nargout>0, varargout{1} = m; return, end

  if nargin <= 2
      [varargout{:}] = g(y,hyp);
  elseif nargin == 3
      if strcmpi(varargin{1},'inv')
          [varargout{:}] = ig(y,hyp);
      elseif isnumeric(varargin{1})
          [varargout{:}] = g(y,hyp,varargin{1});
      end
  else
      error('Derivative of inverse not supported yet.');
  end

% Logarithmic warping function g(y) and log of the derivative log(g'(y))>0
% or derivatives of the latter w.r.t. ith hyperparameter
function [gy,lgpy] = g(y,hyp,i)

  mu = hyp(1);
  delta = exp(hyp(2));

  idx = y > mu;

  if nargin==2                                                 % function values
    gy = y;
    gy(idx) = delta*log((y(idx) - mu)/delta +1)+mu;

    lgpy = zeros(size(y));
    lgpy(idx) = log(delta) - log(y(idx)-mu+delta);
  else                                                          % derivatives

    if i == 1
        gy = zeros(size(y));
        gy(idx) = (y(idx) - mu)./(y(idx) - mu + delta);
        lgpy = zeros(size(y));
        lgpy(idx) = 1./(y(idx)-mu+delta);
    elseif i == 2
        gy = zeros(size(y));
        gy(idx) = delta*(log(1 + (y(idx) - mu)/delta)+ (mu - y(idx))./(y(idx) - mu + delta));
        lgpy = zeros(size(y));
        lgpy(idx) = (y(idx)-mu)./(y(idx)-mu+delta);
    end
  end

% invert g(y)
function [y,n,d] = ig(z,hyp)
  y = z;
  mu = hyp(1);
  delta = exp(hyp(2));
  idx = z > mu;
  y(idx) = delta*(exp((z(idx)-mu)./delta)-1) + mu;
  n = 0;
  d = 0;

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


% Creates a ZC sequence that is mapped onto an OFDM symbol and FFT'd
%
% The results of this function can be used to cross correlate for the specified OFDM
% symbol (numbers 4 or 6).  There is no cyclic prefix added by this function!
%
% 601 samples are created by the ZC, but the middle sample is zeroed out.  This is done
% so that the DC carrier of the FFT is not populated
%
% @param fft_size Size of the OFDM FFT window (must be power of 2)
% @param symbol_index Which ZC sequence symbol should be created (must be 4 or 6)
%
function [samples] = create_zc(fft_size, symbol_index)
    % Validate inputs
    assert(symbol_index == 4 || symbol_index == 6, "Invalid symbol index (must be 4 or 6)");
    assert(log2(fft_size) == round(log2(fft_size)), "Invalid FFT size.  Must be power of 2");

    % Pick the correct root for the ZC sequence
    if (symbol_index == 4)
        root = 600;
    else
        root = 147;
    end

    % Would use MATLAB's zadoffChuSeq function, but Octave doesn't have that
    % The logic below was tested against the MATLAB function
    zc = reshape(exp(-1j * pi * root * (0:600) .* (1:601) / 601), [], 1);

    % Remove the middle value (this would be DC in the FFT)
    zc(301) = [];

    % Create a buffer to hold the freq domain carriers
    samples_freq = zeros(fft_size, 1);

    % Get which FFT bins should be used for data carriers
    data_carrier_indices = get_data_carrier_indices(fft_size * 15e3);

    % Assign just the data carrier bins (left to right) the ZC sequence values
    samples_freq(data_carrier_indices) = zc;

    % Convert to time domain making sure to flip the spectrum left to right first
    samples = ifft(fftshift(samples_freq));
end


text:
function  p = updateboundsfromupper(p,upper)
if nargin == 1
    upper = p.upper;
end
if ~isinf(upper)
    LU = [p.lb p.ub];
    if nnz(p.c.*(p.ub-p.lb)) == 1 & nnz(p.Q)==0
        i = find(p.c.*(p.ub-p.lb));
        if p.c(i) > 0
            % We are minimizing x(i), since an upper bound is UPPER
            % this means c(i)*x(i) has to be < UPPER in optimal solution
            p.ub(i) = min([p.ub(i) upper/p.c(i)]);
        elseif p.c(i) < 0
            % We are maximizing x(i), since an lower bound is -UPPER
            % this means x(i) has to be > -UPPER in optimal solution
            p.lb(i) = max([p.lb(i) -upper/abs(p.c(i))]);
        end

    end
    if ~isempty(p.bilinears) & nnz(p.Q)==0
        quad_v = find(p.bilinears(:,2) == p.bilinears(:,3));
        quad_x = p.bilinears(quad_v,2);
        quad_v = p.bilinears(quad_v,1);
        i = find(p.c(quad_v)>0);
        quad_v = quad_v(i);
        quad_x = quad_x(i);
        if ~isempty(quad_v) & ~any(isinf(p.lb(quad_x))) & ~any(isinf(p.ub(quad_x)))
            % x'*D*x + c'*x + f + h(x) < U
            D = diag(p.c(quad_v));
            if all(diag(D)>0)
                c = p.c(quad_x);
                h = p.c;
                h(quad_x) = 0;
                h(quad_v) = 0;
                % Complete
                % (x + D^-0.5*c/2)'*D*(x + D^-0.5*c/2) - c'*D*c/4 + f + h(x) < U
                xc = -(D^-1)*c/2;
                rhs = upper - p.f - (h(h<0)'*p.ub(h<0) + h(h>=0)'*p.lb(h>=0)) + xc'*D*xc;
                if rhs>0
                    D = D/rhs;
                    % (x-xc)'*D*(x-xc) <= 1
                    for i = 1:length(quad_x)
                        p.lb(quad_x(i)) = max(p.lb(quad_x(i)),xc(i) - 1/sqrt(D(i,i)));
                        p.ub(quad_x(i)) = min(p.ub(quad_x(i)),xc(i) + 1/sqrt(D(i,i)));
                    end
                end
            end
        end
    end
    if ~isequal(LU,[p.lb p.ub])
        p.changedbounds = 1;
    end

    % Some initial code for using inverse objective to derive bounds.
    if nnz(p.Q)==0 && nnz(p.c)==1
        [pos,~,val] = find(p.c);
        if val == -1
            fi = find(p.evalVariables == pos);
            if ~isempty(fi)
                % We're maximizing f(x)
                if ~isempty(p.evalMap{fi}.properties.inverse)
                    if strcmp(p.evalMap{fi}.properties.definiteness,'positive')
                        if strcmp(p.evalMap{fi}.properties.monotonicity,'increasing')
                            if length(p.evalMap{fi}.arg)==2
                                lower = -upper;
                                p.lb(p.evalMap{fi}.variableIndex) = max([p.lb(p.evalMap{fi}.variableIndex) p.evalMap{fi}.properties.inverse(lower)]);
                            end
                        end
                    end
                end
            end
        end
        if val == 1
            fi = find(p.evalVariables == pos);
            if ~isempty(fi)
                % We're minimizing f(x)
                if ~isempty(p.evalMap{fi}.properties.inverse)
                    if strcmp(p.evalMap{fi}.properties.definiteness,'positive')
                        if strcmp(p.evalMap{fi}.properties.monotonicity,'increasing')
                            if length(p.evalMap{fi}.arg)==2
                                p.ub(p.evalMap{fi}.variableIndex) = min([p.ub(p.evalMap{fi}.variableIndex) p.evalMap{fi}.properties.inverse(upper)]);
                            end
                        end
                    end
                end
            end
        end
    end
    if any(p.lb > p.ub + 1e-7)
        p.feasible = 0;
    end
end

author: zarathustr

repo: LibQPEP


text:
function [bar] = Torr2bar(Torr)
% Convert pressure from torr (same as mmHg) to bar
% Chad Greene 2012
bar = Torr*0.00133322;


text:
function [varargout] = routlm(varargin)

% ROUTLM computes the projection of a point from its la/mu parameters
% these equal the "Barycentric" coordinates
%
% Use as
%   [proj] = routlm(v1, v2, v3, la, mu)
% where v1, v2 and v3 are three vertices of the triangle

% Copyright (C) 2002-2009, 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$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is the native MATLAB implementation.
% The mex file is many times faster and is therefore preferred.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [proj] = routlm(v1, v2, v3, la, mu);
% proj = (1-la-mu)*v1 + la*v2 + mu*v3;

% compile the missing mex file on the fly
% remember the original working directory
pwdir = pwd;

% determine the name and full path of this function
funname = mfilename('fullpath');
mexsrc  = [funname '.c'];
[mexdir, mexname] = fileparts(funname);

try
  % try to compile the mex file on the fly
  warning('trying to compile MEX file from %s', mexsrc);
  cd(mexdir);

  if ispc
    mex -I. -c geometry.c
    mex -I. -c routlm.c ; mex routlm.c routlm.obj geometry.obj
  else
    mex -I. -c geometry.c
    mex -I. -c routlm.c ; mex -o routlm routlm.o geometry.o
  end

  cd(pwdir);
  success = true;

catch
  % compilation failed
  disp(lasterr);
  error('could not locate MEX file for %s', mexname);
  cd(pwdir);
  success = false;
end

if success
  % execute the mex file that was just created
  funname   = mfilename;
  funhandle = str2func(funname);
  [varargout{1:nargout}] = funhandle(varargin{:});
end


text:
%% Example of using KST class for interfacing with KUKA iiwa robots
% Controlling KUKA iiwa from V-rep simulation using MATLAB and KST as a
% middle-ware.
% This script works together with V-rep scene: iiwaFromVrep.ttt

%% To run this example do the following:
% 1- Run the Vrep simulation: iiwaFromVrep.ttt
% 2- Run the (MarlabToolboxServer) application on the robot
% 3- Run this script from MATLAB.

% Copyright: Mohammad SAFEEA, 25th of June 2018

%% Initiation part of the code
close all;clear;clc;
warning('off')
disp('Program started');
%% Add path of KST class to work space
cDir = pwd;
cDir=getTheKSTDirectory(cDir);
addpath(cDir);

%% Declare V-rep objects
vrep=remApi('remoteApi');
vrep.simxFinish(-1);
clientID=vrep.simxStart('127.0.0.1',19999,true,true,5000,5);

%% Instantiate the KST object
ip='172.31.1.147'; % The IP of the controller
arg1=KST.LBR7R800; % choose the robot iiwa7R800 or iiwa14R820
arg2=KST.Medien_Flansch_elektrisch; % choose the type of flange
Tef_flange=eye(4); % transofrm matrix of EEF with respect to flange
iiwa=KST(ip,arg1,arg2,Tef_flange); % create the object

if (clientID>-1)
    %% Start a connection with the server
    flag=iiwa.net_establishConnection();
    if flag==0
      return;
    end

    jHandles=zeros(7,1);
    for i=1:7
            s=['LBR_iiwa_7_R800_joint',num2str(i)];
            [res, daHandle]=vrep.simxGetObjectHandle(clientID,s,vrep.simx_opmode_oneshot_wait);
            jHandles(i)=daHandle;
    end
else
    return;
end
jPos={0., pi / 180 * 20., 0., -pi / 180 * 70., 0.,...
                        pi / 180 * 90., 0.};
relVel=0.25;
iiwa.movePTPJointSpace(jPos, relVel);
%% Control loop
iiwa.realTime_startDirectServoJoints();
pause(10);
counter=0;
tic;
totalTimeSecs=60; % 60 seconds of execution time
vrep.simxStartSimulation(clientID,vrep.simx_opmode_oneshot_wait);
for i=1:7
        [res,tempPos]=vrep.simxGetJointPosition(clientID,jHandles(i),...
            vrep.simx_opmode_streaming);
end
vrep.simxWriteStringStream(clientID,'path_velocity','0.07',vrep.simx_opmode_oneshot_wait);
jpos=zeros(7,25000);
t_0=toc;
while(t_0<totalTimeSecs)
    for i=1:7
            [res,tempPos]=vrep.simxGetJointPosition(clientID,jHandles(i),...
                vrep.simx_opmode_buffer);
            jPos{i}=tempPos;
    end

    if(toc-t_0>0.003)
        counter=counter+1;
        iiwa.sendJointsPositionsf(jPos);
        for i=1:7
            jpos(i,counter)=jPos{i};
        end
        t_0=toc;
	end
end
%% Turn off KST server
iiwa.net_turnOffServer();
%% Stop simulation
vrep.simxStopSimulation(clientID,vrep.simx_opmode_oneshot_wait);
%% End simulation
vrep.simxFinish(clientID);
vrep.delete();

repo: KST-Kuka-Sunrise-Toolbox


text:
function varargout = max(varargin)
%MAX   Maximum value of a CHEBFUN2 in one direction.
%   MAX(f) returns a chebfun representing the maximum of the CHEBFUN2 along the
%   y direction, i.e, MAX(f) = @(x) max( f ( x, : ) )
%
%   MAX(f, [], dim) returns a CHEBFUN representing the maximum of f along the
%   DIM direction. If DIM = 1 is along the y-direction and DIM = 2 is along the
%   x-direction.
%
%   WARNING: This function is not always accurate to the expected precision.
%
%   For the global maximum use MAX2.
%
% See also MIN, MAX2, MIN2, MINANDMAX2.

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

[varargout{1:nargout}] = max@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


function [ L b info ] = ldml_learn( X, Y, k, it, verbose, A0 )
% [ L b info ] = ldml_learn( X, Y, k, it, verbose, A0 )
%
%  * Mandatory inputs:
% X : (m x d) data matrix (m data points with d dimensions)
% Y : (m x 1) class labels
%
%  * To learn a kernelized metric, provide the kernel Gram matrix instead of the data matrix X.
%
%  * Optional inputs:
% k : number of dimensions of projection space (default=d)
% it : number of iterations (default=100)
% verbose : boolean for verbosity (default=false)
% A0 : parameter initialization vector (default=random)
%
%  * Outputs:
% L : (d x k) projection matrix
% b : (scalar) distance threshold as learnt by LDML (may or may not be optimal depending on application)
% info : struct with information about the optimization, with following fields
%     .fA : sequence of function values during gradient descent
%     .it : effective number of iterations
%     .A  : parameter vector that you can re-use for initialization (cf "A0" in input section)
%
% Matthieu Guillaumin, INRIA.

[ n d ] = size(X);
m = numel(Y);
I=(1:n);

if n ~= m,
   error('ldml','Second input must have a number of elements equal to the number of rows of the first input');
end

YY=sparse(1:n,Y(:),true);
c=size(YY,2);

% Set default values
opt=2;
if nargin<opt+1,    k=d;                    end
if nargin<opt+2,    it=100;                 end
if nargin<opt+3,    verbose=false;          end

if verbose, fprintf('Learning (%d x %d) LDML with %d instances in %d bags with %d labels\n',d,k,n,m,c); end

% Set initial projection as a random orthogonal matrix. TODO: several random initialization.
if nargin<opt+4,
    A0 = orth(rand(d,k));
    A0 = A0(:);
    A0(end+1)=rand;
end

H=vec(hist_count(I));
if m>size(H,1),
    H(size(H,1)+1:m)=0;
end

% Perform gradient descent
[ A fA it ] = minimize( A0, 'mildml_fg', it, verbose, double(X), double(I)-1, double(sparse(YY)), double(H), n, d, k, m, c);
info.A=A;info.fA=fA;info.it=it;

% Get projection matrix and bias from learned parameter vector
b = A(end);
L = reshape(A(1:end-1),d,k);


text:
% STD_READSPEC - load spectrum measures for data channels or
%                  for all components of a specified cluster.
%                  Called by plotting functions
%                  STD_ENVTOPO, STD_ERPPLOT, STD_ERSPPLOT, ...
% Usage:
%         >> [STUDY, specdata, allfreqs, setinds, cinds] = ...
%                   std_readspec(STUDY, ALLEEG, varargin);
% Inputs:
%       STUDY - studyset structure containing some or all files in ALLEEG
%      ALLEEG - vector of loaded EEG datasets
%
% Optional inputs:
%  'design'    - [integer] read files from a specific STUDY design. Default
%                is empty (use current design in STUDY.currentdesign).
%  'channels'  - [cell] list of channels to import {default: none}
%  'clusters'  - [integer] list of clusters to import {[]|default: all but
%                the parent cluster (1) and any 'NotClust' clusters}
%  'singletrials' - ['on'|'off'] load single trials spectral data (if
%                available). Default is 'off'.
%  'subject'   - [string] select a specific subject {default:all}
%  'component' - [integer] select a specific component in a cluster
%                 {default:all}
%
% Spectrum specific inputs:
%  'freqrange' - [min max] frequency range {default: whole measure range}
%  'rmsubjmean' - ['on'|'off'] remove mean subject spectrum from every
%                 channel spectrum, making them easier to compare
%                 { default: 'off' }
% Output:
%  STUDY    - updated studyset structure
%  specdata - [cell array] spectral data (the cell array size is
%             condition x groups)
%  freqs    - [float array] array of frequencies
%  setinds  - [cell array] datasets indices
%  cinds    - [cell array] channel or component indices
%
% Example:
%  std_precomp(STUDY, ALLEEG, { ALLEEG(1).chanlocs.labels }, 'spec', 'on');
%  [spec freqs] = std_readspec(STUDY, ALLEEG, 'channels', { ALLEEG(1).chanlocs(1).labels });
%
% Author: Arnaud Delorme, CERCO, 2006-

% Copyright (C) Arnaud Delorme, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [STUDY, specdata, allfreqs] = std_readspec(STUDY, ALLEEG, varargin)

if nargin < 2
    help std_readspec;
    return;
end

[STUDY, specdata, allfreqs] = std_readerp(STUDY, ALLEEG, 'datatype', 'spec', varargin{:});

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,slice,n,localglobal,cst)
% gamma index calculation
% according to http://www.ncbi.nlm.nih.gov/pubmed/9608475
%
% call
%   [gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,cst)
%   [gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,slice,cst)
%   [gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,n,cst)
%   [gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,localglobal,cst)
%   [gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,slice,n,cst)
%               ...
%   [gammaCube,gammaPassRateCell] = matRad_gammaIndex(cube1,cube2,resolution,criteria,slice,n,localglobal,cst)
%
% input
%   cube1:         dose cube as an M x N x O array
%   cube2:         dose cube as an M x N x O array
%   resolution:    resolution of the cubes [mm/voxel]
%   criteria:      [1x2] vector specifying the distance to agreement
%                  criterion; first element is percentage difference,
%                  second element is distance [mm]
%   slice:         (optional) slice in cube1/2 that will be visualized
%   n:             (optional) number of interpolations. there will be 2^n-1
%                  interpolation points. The maximum suggested value is 3.
%   localglobal:   (optional) parameter to choose between 'global' and 'local'
%                  normalization
%   cst:           list of interessing volumes inside the patient
%
% output
%
%   gammaCube:          result of gamma index calculation
%   gammaPassRateCell:  rate of voxels passing the specified gamma criterion
%                  evaluated for every structure listed in 'cst'.
%                  note that only voxels exceeding the dose threshold are
%                  considered.
%
% References
%   [1]  http://www.ncbi.nlm.nih.gov/pubmed/9608475
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[env, ~] = matRad_getEnvironment();

% set parameters for gamma index calculation
if exist('criteria','var')
    relDoseThreshold = criteria(1); % in [%]
    dist2AgreeMm     = criteria(2); % in [mm]
else
    dist2AgreeMm     = 3; % in [mm]
    relDoseThreshold = 3; % in [%]
end

% set parameters for gamma index calculation
if ~exist('n','var')
    n = 0;
end

% set global or local gamma evaluation
if ~exist('localglobal','var')
    localglobal = 'global';
end

% check if cubes consistent
if ~isequal(size(cube1),size(cube2))
   error('dose cubes must be the same size\n');
end

% define search nighborhood and resolution
if round(n) ~= n
    error('n must be an integer value');
else
    resolution = resolution./(2.^n);
end

neighborhood = 1.5* dist2AgreeMm; % [mm]
searchX = ceil(neighborhood./resolution(1));
searchY = ceil(neighborhood./resolution(2));
searchZ = ceil(neighborhood./resolution(3));

% init cube
% cut all the zeros surrounding the interesting part
cut1 = any(any(cube1,2),3);
cut2 = any(any(cube1,1),3);
cut3 = any(any(cube1,1),2);

% avoids that "zero-slides" are deleted between two interest regions
k = find(cut1); cut1(k(1):k(end)) = 1;
k = find(cut2); cut2(k(1):k(end)) = 1;
k = find(cut2); cut2(k(1):k(end)) = 1;

% copy cubes
cubex1c = cube1;
cubex2c = cube2;

% cut cubes
cubex1c( ~cut1, :, :) = []; % rows
cubex1c( :, ~cut2, :) = []; % columns
cubex1c( :, :, ~cut3) = []; % slices
cubex2c( ~cut1, :, :) = []; % rows
cubex2c( :, ~cut2, :) = []; % columns
cubex2c( :, :, ~cut3) = []; % slices

% pad cubes
cubex2 = zeros([size(cubex2c,1)+2*searchX size(cubex2c,2)+2*searchY size(cubex2c,3)+2*searchZ]);
cubex2((1+searchX):(end-searchX), ...
       (1+searchY):(end-searchY), ...
       (1+searchZ):(end-searchZ)) = cubex2c;

% interpolate if necessary
if n > 0
    switch env
        case 'MATLAB'
            cubex2 = interp3(cubex2,n,'cubic');
        case 'OCTAVE'
            cubex2 = interp3(cubex2,n,'linear');
    end
end

% set up temporary cubes required for calculation
gammaCubeSq = inf*ones(size(cubex1c));

% adjust dose threshold
if strcmp(localglobal,'local')
    doseThreshold = cubex1c .* relDoseThreshold/100;
elseif strcmp(localglobal,'global')
    doseThreshold = relDoseThreshold/100 * max(cube1(:));
end

% search for min
for i = -searchX:searchX
    for j = -searchY:searchY
        for k = -searchZ:searchZ

            delta_sq = ((i*resolution(1))^2 + ...
                        (j*resolution(2))^2 + ...
                        (k*resolution(3))^2) / dist2AgreeMm^2;

            tmpCube = cubex1c - cubex2((1+((2^n)*searchX)+i) : 2^n : (end-((2^n)*searchX)+i), ...
                                       (1+((2^n)*searchY)+j) : 2^n : (end-((2^n)*searchY)+j), ...
                                       (1+((2^n)*searchZ)+k) : 2^n : (end-((2^n)*searchZ)+k));

            tmpCube = tmpCube.^2 ./ doseThreshold.^2 + delta_sq;

            gammaCubeSq = min(gammaCubeSq,tmpCube);

        end
    end

%     display '.';

end

% evaluate gamma cube and set to zero all the voxel that contain an
% infinite value
gammaCubeSqx                 = zeros(size(cube1));
gammaCubeSqx(cut1,cut2,cut3) = gammaCubeSq;
gammaCube                    = sqrt(gammaCubeSqx);

% set values where we did not compute gamma keeping inf in the cube to zero
gammaCube(cube1<=0 & cube2<=0) = 0;

% compute gamma pass rate
doseIx          = cube1 > relDoseThreshold/100*max(cube1(:)) | cube2 > relDoseThreshold/100*max(cube2(:));
numOfPassGamma  = sum(gammaCube(doseIx) < 1);
gammaPassRate   = 100 * numOfPassGamma / sum(doseIx(:));

% compute stats for all segmented sturtures
if exist('cst','var')
    gammaPassRateCell = cell(1,2);
    gammaPassRateCell{1,1} = 'Whole CT';
    gammaPassRateCell{1,2} = gammaPassRate;

    for i = 1:size(cst, 1)
        volume = cst{i,4}{1,1}; % indices of voxels of the interesting volume
        doseIxVol = false(size(doseIx));
        doseIxVol(volume) = doseIx(volume);
        numOfPassGammaVol  = sum(gammaCube(doseIxVol) < 1);
        gammaPassRateVol   = 100 * numOfPassGammaVol / sum(doseIxVol(:));
        if isnan(gammaPassRateVol) % remove organs not receiving any dose (resulting in NaN)
            continue;
        end
        gammaPassRateCell{end+1, 1} = cst{i,2};
        gammaPassRateCell{end, 2} = gammaPassRateVol;

    end
end

% visualize if applicable
if exist('slice','var') && ~isempty(slice)
    figure
    set(gcf,'Color',[1 1 1]);
    imagesc(gammaCube(:,:,slice),[0 2])
    myColormap = matRad_getColormap('gammaIndex');

    colormap(gca,myColormap);
    colorbar

    title({[num2str(gammaPassRate,5) '% of points > ' num2str(relDoseThreshold) ...
            '% pass gamma criterion (' num2str(relDoseThreshold) '% / ' ...
            num2str(dist2AgreeMm) 'mm)']; ['with ' num2str(2^n-1) ' interpolation points']});
end

repo: matRad


text:
function [varargout] = likGauss(hyp, y, mu, s2, inf, i)

% likGauss - Gaussian likelihood function for regression. The expression for the
% likelihood is
%   likGauss(t) = exp(-(t-y)^2/2*sn^2) / sqrt(2*pi*sn^2),
% where y is the mean and sn is the standard deviation.
%
% The hyperparameters are:
%
% hyp = [  log(sn)  ]
%
% Several modes are provided, for computing likelihoods, derivatives and moments
% respectively, see likFunctions.m for the details. In general, care is taken
% to avoid numerical issues when the arguments are extreme.
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2015-07-13.
%                                      File automatically generated using noweb.
%
% See also LIKFUNCTIONS.M.

if nargin<3, varargout = {'1'}; return; end   % report number of hyperparameters

sn2 = exp(2*hyp);

if nargin<5                              % prediction mode if inf is not present
  if isempty(y),  y = zeros(size(mu)); end
  s2zero = 1; if nargin>3&&numel(s2)>0&&norm(s2)>eps, s2zero = 0; end  % s2==0 ?
  if s2zero                                                    % log probability
    lp = -(y-mu).^2./sn2/2-log(2*pi*sn2)/2; s2 = 0;
  else
    lp = likGauss(hyp, y, mu, s2, 'infEP');                         % prediction
  end
  ymu = {}; ys2 = {};
  if nargout>1
    ymu = mu;                                                   % first y moment
    if nargout>2
      ys2 = s2 + sn2;                                          % second y moment
    end
  end
  varargout = {lp,ymu,ys2};
else
  switch inf
  case 'infLaplace'
    if nargin<6                                             % no derivative mode
      if isempty(y), y=0; end
      ymmu = y-mu; dlp = {}; d2lp = {}; d3lp = {};
      lp = -ymmu.^2/(2*sn2) - log(2*pi*sn2)/2;
      if nargout>1
        dlp = ymmu/sn2;                      % dlp, derivative of log likelihood
        if nargout>2                    % d2lp, 2nd derivative of log likelihood
          d2lp = -ones(size(ymmu))/sn2;
          if nargout>3                  % d3lp, 3rd derivative of log likelihood
            d3lp = zeros(size(ymmu));
          end
        end
      end
      varargout = {lp,dlp,d2lp,d3lp};
    else                                                       % derivative mode
      lp_dhyp = (y-mu).^2/sn2 - 1;  % derivative of log likelihood w.r.t. hypers
      dlp_dhyp = 2*(mu-y)/sn2;                               % first derivative,
      d2lp_dhyp = 2*ones(size(mu))/sn2;   % and also of the second mu derivative
      varargout = {lp_dhyp,dlp_dhyp,d2lp_dhyp};
    end

  case 'infEP'
    if nargin<6                                             % no derivative mode
      lZ = -(y-mu).^2./(sn2+s2)/2 - log(2*pi*(sn2+s2))/2;    % log part function
      dlZ = {}; d2lZ = {};
      if nargout>1
        dlZ  = (y-mu)./(sn2+s2);                    % 1st derivative w.r.t. mean
        if nargout>2
          d2lZ = -1./(sn2+s2);                      % 2nd derivative w.r.t. mean
        end
      end
      varargout = {lZ,dlZ,d2lZ};
    else                                                       % derivative mode
      dlZhyp = ((y-mu).^2./(sn2+s2)-1) ./ (1+s2./sn2);   % deriv. w.r.t. hyp.lik
      varargout = {dlZhyp};
    end

  case 'infVB'
    % variational lower site bound
    % t(s) = exp(-(y-s)^2/2sn2)/sqrt(2*pi*sn2)
    % the bound has the form: (b+z/ga)*f - f.^2/(2*ga) - h(ga)/2
    n = numel(s2); b = zeros(n,1); y = y.*ones(n,1); z = y;
    varargout = {b,z};
  end
end


text:
function y = pow_p( x, p )

%POW_P   Internal cvx version.

error(nargchk(2,2,nargin));
y = pow_cvx( x, p, 'pow_p' );

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

repo: REST