text:
function SaveFig(fname, varargin)
% OLD now use export_fig
% Usage: savefig(filename, fighdl, options)
%
% Saves a pdf, eps, png, jpeg, and/or tiff of the contents of the fighandle's (or current) figure.
% It saves an eps of the figure and the uses Ghostscript to convert to the other formats.
% The result is a cropped, clean picture. There are options for using rgb or cmyk colours,
% or grayscale. You can also choose the resolution.
%
% The advantage of savefig is that there is very little empty space around the figure in the
% resulting files, you can export to more than one format at once, and Ghostscript generates
% trouble-free files.
%
% If you find any errors, please let me know! (peder at axensten dot se)
%
% filename: File name without suffix.
%
% fighdl:  (default: gcf) Integer handle to figure.
%
% options: (default: '-r300', '-lossless', '-rgb') You can define your own
%          defaults in a global variable savefig_defaults, if you want to, i.e.
%          savefig_defaults= {'-r200','-gray'};.
% 'eps':   Output in Encapsulated Post Script (no preview yet).
% 'pdf':   Output in (Adobe) Portable Document Format.
% 'png':   Output in Portable Network Graphics.
% 'jpeg':  Output in Joint Photographic Experts Group format.
% 'tiff':  Output in Tagged Image File Format (no compression: huge files!).
% '-rgb':  Output in rgb colours.
% '-cmyk': Output in cmyk colours (not yet 'png' or 'jpeg' -- '-rgb' is used).
% '-gray': Output in grayscale (not yet 'eps' -- '-rgb' is used).
% '-fonts':  Include fonts in eps or pdf. Includes only the subset needed.
% '-lossless':  Use lossless compression, works on most formats. same as '-c0', below.
% '-c<float>':  Set compression for non-indexed bitmaps in PDFs -
%               0: lossless; 0.1: high quality; 0.5: medium; 1: high compression.
% '-r<integer>':  Set resolution.
% '-crop': Removes points and line segments outside the viewing area -- permanently.
%          Only use this on figures where many points and/or line segments are outside
%          the area zoomed in to. This option will result in smaller vector files (has no
%          effect on pixel files).
% '-dbg':  Displays gs command line(s).
%
% EXAMPLE:
% savefig('nicefig', 'pdf', 'jpeg', '-cmyk', '-c0.1', '-r250');
% Saves the current figure to nicefig.pdf and nicefig.png, both in cmyk and at 250 dpi,
%          with high quality lossy compression.
%
% REQUIREMENT: Ghostscript. Version 8.57 works, probably older versions too, but '-dEPSCrop'
%          must be supported. I think version 7.32 or newer is ok.
%
% HISTORY:
% Version 1.0, 2006-04-20.
% Version 1.1, 2006-04-27:
% - No 'epstopdf' stuff anymore! Using '-dEPSCrop' option in gs instead!
% Version 1.2, 2006-05-02:
% - Added a '-dbg' option (see options, above).
% - Now looks for a global variable 'savefig_defaults' (see options, above).
% - More detailed Ghostscript options (user will not really notice).
% - Warns when there is no device for a file-type/color-model combination.
% Version 1.3, 2006-06-06:
% - Added a check to see if there actually is a figure handle.
% - Now works in Matlab 6.5.1 (R13SP1) (maybe in 6.5 too).
% - Now compatible with Ghostscript 8.54, released 2006-06-01.
% Version 1.4, 2006-07-20:
% - Added an option '-soft' that enables anti-aliasing on pixel graphics (on by default).
% - Added an option '-hard' that don't do anti-aliasing on pixel graphics.
% Version 1.5, 2006-07-27:
% - Fixed a bug when calling with a figure handle argument.
% Version 1.6, 2006-07-28:
% - Added a crop option, see above.
% Version 1.7, 2007-03-31:
% - Fixed bug: calling print with invalid renderer value '-none'.
% - Removed GhostScript argument '-dUseCIEColor' as it sometimes discoloured things.
% Version 1.8, 2008-01-03:
% - Added MacIntel: 'MACI'.
% - Added 64bit PC (I think, can't test it myself).
% - Added option '-nointerpolate' (use it to prevent blurring of pixelated).
% - Removed '-hard' and '-soft'. Use '-nointerpolate' for '-hard', default for '-soft'.
% - Fixed the gs 8.57 warning on UseCIEColor (it's now set).
% - Added '-gray' for pdf, but gs 8.56 or newer is needed.
% - Added '-gray' and '-cmyk' for eps, but you a fairly recent gs might be needed.
% Version 1.9, 2008-07-27:
% - Added lossless compression, see option '-lossless', above. Works on most formats.
% - Added lossy compression, see options '-c<float>...', above. Works on 'pdf'.
%   Thanks to Olly Woodford for idea and implementation!
% - Removed option '-nointerpolate' -- now savefig never interpolates.
% - Fixed a few small bugs and removed some mlint comments.
% Version 2.0, 2008-11-07:
% - Added the possibility to include fonts into eps or pdf.
%
% TO DO: (Need Ghostscript support for these, so don't expect anything soon...)
% - svg output.
% - '-cmyk' for 'jpeg' and 'png'.
% - Preview in 'eps'.
% - Embedded vector fonts, not bitmap, in 'eps'.
%
% Copyright (C) Peder Axensten (peder at axensten dot se), 2006.

% KEYWORDS:     eps, pdf, jpg, jpeg, png, tiff, eps2pdf, epstopdf, ghostscript
%
% INSPIRATION:  eps2pdf (5782), eps2xxx (6858)
%
% REQUIREMENTS: Works in Matlab 6.5.1 (R13SP1) (maybe in 6.5 too).

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

	op_dbg=		false;												% Default value.

	% Compression
	compr=			[' -dUseFlateCompression=true -dLZWEncodePages=true -dCompatibilityLevel=1.6' ...
					 ' -dAutoFilterColorImages=false -dAutoFilterGrayImages=false ' ...
					 ' -dColorImageFilter=%s -dGrayImageFilter=%s'];	% Compression.
	lossless=		sprintf (compr, '/FlateEncode', '/FlateEncode');
	lossy=			sprintf (compr, '/DCTEncode',   '/DCTEncode'  );
	lossy=			[lossy ' -c ".setpdfwrite << /ColorImageDict << /QFactor %g ' ...
					 '/Blend 1 /HSample [%s] /VSample [%s] >> >> setdistillerparams"'];

	% Create gs command.
	cmdEnd=			' -sDEVICE=%s -sOutputFile="%s"';					% Essential.
	epsCmd=			'';
	epsCmd=	[epsCmd ' -dSubsetFonts=true -dNOPLATFONTS'];				% Future support?
	epsCmd=	[epsCmd ' -dUseCIEColor=true -dColorConversionStrategy=/UseDeviceIndependentColor'];
	epsCmd=	[epsCmd ' -dProcessColorModel=/%s'];						% Color conversion.
	pdfCmd=	[epsCmd ' -dAntiAliasColorImages=false' cmdEnd];
	epsCmd=	[epsCmd cmdEnd];

	% Get file name.
	if((nargin < 1) || isempty(fname) || ~ischar(fname))				% Check file name.
		error('No file name specified.');
	end
	[pathstr, namestr] = fileparts(fname);
	if(isempty(pathstr)), fname= fullfile(cd, namestr);	end

	% Get handle.
	fighdl=		get(0, 'CurrentFigure'); % See gcf.						% Get figure handle.
	if((nargin >= 2) && (numel(varargin{1}) == 1) && isnumeric(varargin{1}))
		fighdl=		varargin{1};
		varargin=	{varargin{2:end}};
	end
	if(isempty(fighdl)), error('There is no figure to save!?');	end
	set(fighdl, 'Units', 'centimeters')									% Set paper stuff.
	sz=			get(fighdl, 'Position');
	sz(1:2)=	0;
	set(fighdl, 'PaperUnits', 'centimeters', 'PaperSize', sz(3:4), 'PaperPosition', sz);

	% Set up the various devices.
	% Those commented out are not yet supported by gs (nor by savefig).
	% pdf-cmyk works due to the Matlab '-cmyk' export being carried over from eps to pdf.
	device.eps.rgb=		sprintf(epsCmd, 'DeviceRGB',	'epswrite', [fname '.eps']);
	device.jpeg.rgb=	sprintf(cmdEnd,	'jpeg', 					[fname '.jpeg']);
%	device.jpeg.cmyk=	sprintf(cmdEnd,	'jpegcmyk', 				[fname '.jpeg']);
	device.jpeg.gray=	sprintf(cmdEnd,	'jpeggray',					[fname '.jpeg']);
	device.pdf.rgb=		sprintf(pdfCmd, 'DeviceRGB',	'pdfwrite', [fname '.pdf']);
	device.pdf.cmyk=	sprintf(pdfCmd, 'DeviceCMYK',	'pdfwrite', [fname '.pdf']);
	device.pdf.gray=	sprintf(pdfCmd, 'DeviceGray',	'pdfwrite', [fname '.pdf']);
	device.png.rgb=		sprintf(cmdEnd,	'png16m', 					[fname '.png']);
%	device.png.cmyk=	sprintf(cmdEnd,	'png???', 					[fname '.png']);
	device.png.gray=	sprintf(cmdEnd,	'pnggray', 					[fname '.png']);
	device.tiff.rgb=	sprintf(cmdEnd,	'tiff24nc',					[fname '.tiff']);
	device.tiff.cmyk=	sprintf(cmdEnd,	'tiff32nc', 				[fname '.tiff']);
	device.tiff.gray=	sprintf(cmdEnd,	'tiffgray', 				[fname '.tiff']);

	% Get options.
	global savefig_defaults;											% Add global defaults.
 	if( iscellstr(savefig_defaults)), varargin=	{savefig_defaults{:}, varargin{:}};
	elseif(ischar(savefig_defaults)), varargin=	{savefig_defaults, varargin{:}};
	end
	varargin=	{'-r300', '-lossless', '-rgb', varargin{:}};			% Add defaults.
	res=		'';
	types=		{};
	fonts=		'false';
	crop=		false;
	for n= 1:length(varargin)											% Read options.
		if(ischar(varargin{n}))
			switch(lower(varargin{n}))
			case {'eps','jpeg','pdf','png','tiff'},		types{end+1}=	lower(varargin{n});
			case '-rgb',				color=	'rgb';	deps= {'-depsc2'};
			case '-cmyk',				color=	'cmyk';	deps= {'-depsc2', '-cmyk'};
			case '-gray',				color=	'gray';	deps= {'-deps2'};
			case '-fonts',				fonts=			'true';
			case '-lossless',			comp=			0;
			case '-crop',				crop=			true;
			case '-dbg',				op_dbg=			true;
			otherwise
				if(regexp(varargin{n}, '^\-r[0-9]+$')), 	 res=  varargin{n};
				elseif(regexp(varargin{n}, '^\-c[0-9.]+$')), comp= str2double(varargin{n}(3:end));
				else	warning('pax:savefig:inputError', 'Unknown option in argument: ''%s''.', varargin{n});
				end
			end
		else
			warning('pax:savefig:inputError', 'Wrong type of argument: ''%s''.', class(varargin{n}));
		end
	end
	types=		unique(types);
	if(isempty(types)), error('No output format given.');	end

	if (comp == 0)														% Lossless compression
		gsCompr=		lossless;
	elseif (comp <= 0.1)												% High quality lossy
		gsCompr=		sprintf(lossy, comp, '1 1 1 1', '1 1 1 1');
	else																% Normal lossy
		gsCompr=		sprintf(lossy, comp, '2 1 1 2', '2 1 1 2');
	end

	% Generate the gs command.
	switch(computer)													% Get gs command.
		case {'MAC','MACI'},			gs= '/usr/local/bin/gs';
		case {'PCWIN','PCWIN64'},		gs= 'gswin32c.exe';
		otherwise,						gs= 'gs';
	end
	gs=		[gs		' -q -dNOPAUSE -dBATCH -dEPSCrop'];					% Essential.
	gs=		[gs     ' -dPDFSETTINGS=/prepress -dEmbedAllFonts=' fonts];	% Must be first?
	gs=		[gs		' -dUseFlateCompression=true'];						% Useful stuff.
	gs=		[gs		' -dAutoRotatePages=/None'];						% Probably good.
	gs=		[gs		' -dHaveTrueTypes'];								% Probably good.
	gs=		[gs		' ' res];											% Add resolution to cmd.

	if(crop && ismember(types, {'eps', 'pdf'}))							% Crop the figure.
		fighdl= do_crop(fighdl);
	end

	% Output eps from Matlab.
	renderer=	['-' lower(get(fighdl, 'Renderer'))];					% Use same as in figure.
	if(strcmpi(renderer, '-none')), renderer=	'-painters';	end		% We need a valid renderer.
	print(fighdl, deps{:}, '-noui', renderer, res, [fname '-temp']);	% Output the eps.

	% Convert to other formats.
	for n= 1:length(types)												% Output them.
		if(isfield(device.(types{n}), color))
			cmd=		device.(types{n}).(color);						% Colour model exists.
		else
			cmd=		device.(types{n}).rgb;							% Use alternative.
			if(~strcmp(types{n}, 'eps'))	% It works anyways for eps (VERY SHAKY!).
				warning('pax:savefig:deviceError', ...
						'No device for %s using %s. Using rgb instead.', types{n}, color);
			end
		end
		cmp=	lossless;
		if (strcmp(types{n}, 'pdf')),	cmp= gsCompr;		end			% Lossy compr only for pdf.
		if (strcmp(types{n}, 'eps')),	cmp= '';			end			% eps can't use lossless.
		cmd=	sprintf('%s %s %s -f "%s-temp.eps"', gs, cmd, cmp, fname);% Add up.
			status= system(cmd);										% Run Ghostscript.
			if (op_dbg || status), 		display (cmd),		end
	end
	delete([fname '-temp.eps']);										% Clean up.
end

function fig= do_crop(fig)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%	Remove line segments that are outside the view.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	haxes=	findobj(fig, 'Type', 'axes', '-and', 'Tag', '');
	for n=1:length(haxes)
		xl=		get(haxes(n), 'XLim');
		yl=		get(haxes(n), 'YLim');
		lines=	findobj(haxes(n), 'Type', 'line');
		for m=1:length(lines)
			x=				get(lines(m), 'XData');
			y=				get(lines(m), 'YData');

			inx=			(xl(1) <= x) & (x <= xl(2));	% Within the x borders.
			iny=			(yl(1) <= y) & (y <= yl(2));	% Within the y borders.
			keep=			inx & iny;						% Within the box.

			if(~strcmp(get(lines(m), 'LineStyle'), 'none'))
				crossx=		((x(1:end-1) < xl(1)) & (xl(1) < x(2:end))) ...	% Crossing border x1.
						|	((x(1:end-1) < xl(2)) & (xl(2) < x(2:end))) ...	% Crossing border x2.
						|	((x(1:end-1) > xl(1)) & (xl(1) > x(2:end))) ...	% Crossing border x1.
						|	((x(1:end-1) > xl(2)) & (xl(2) > x(2:end)));	% Crossing border x2.
				crossy=		((y(1:end-1) < yl(1)) & (yl(1) < y(2:end))) ...	% Crossing border y1.
						|	((y(1:end-1) < yl(2)) & (yl(2) < y(2:end))) ...	% Crossing border y2.
						|	((y(1:end-1) > yl(1)) & (yl(1) > y(2:end))) ...	% Crossing border y1.
						|	((y(1:end-1) > yl(2)) & (yl(2) > y(2:end)));	% Crossing border y2.
				crossp=	[(	(crossx & iny(1:end-1) & iny(2:end)) ...	% Crossing a x border within y limits.
						|	(crossy & inx(1:end-1) & inx(2:end)) ...	% Crossing a y border within x limits.
						|	crossx & crossy ...							% Crossing a x and a y border (corner).
						 ),	false ...
						];
				crossp(2:end)=	crossp(2:end) | crossp(1:end-1);		% Add line segment's secont end point.

				keep=			keep | crossp;
			end
			set(lines(m), 'XData', x(keep))
			set(lines(m), 'YData', y(keep))
		end
	end
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
% Correlation map post processing in PIVlab
function [u_out,v_out] = PIVlab_correlation_filter (u,v,corr2_thresh,corr2_value)
u(corr2_value<corr2_thresh)=nan;
v(corr2_value<corr2_thresh)=nan;
u_out=u;
v_out=v;

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function cortMag = buildCortMagBins(cortMag, flatView, grayView)
%
%  cortmag = buildCortMagBins(cortMag, flatView, grayView)
%
%  Loads the bins for all the ROIs. This is a necessary and computationally
%  expensive
%   step torward computing the cortical magnification function. We should
%   avoid doing it
%   more often than necessary.
%
% HISTORY:
%   2002.03.20 RFD (bob@white.stanford.edu) wrote it, based
%   on code by Wandell, Baseler and Brewer.
%   2002.12.11 RFD: added flatDataFlag check so that we can
%   either use the data pulled from the flat map, or get the
%   data from the volume.

% The basic steps that we need to do:
% * find a start point
%    - ROI coords are already sorted, so we just take the first coord
% * find a gray node for each point
% * find distance from start node to every other node
% * grab the co and ph for each node.

bins = {};

% if(strcmp(cortMag.hemisphere,'left'))
%     nodes = grayView.allLeftNodes;
%     edges = grayView.allLeftEdges;
% else
%     nodes = grayView.allRightNodes;
%     edges = grayView.allRightEdges;
% end
nodes = grayView.nodes;
edges = grayView.edges;

for(roiNum=1:length(cortMag.ROIs))
    disp(['Binning nodes for ROI ',num2str(roiNum),' of ',num2str(length(cortMag.ROIs)),'...']);

    % We loop until we've dealt with all the nodes. The basic idea is to
    % pick a binCenterNode, bin all nodes within binDist of that
    % binCenterNode,
    % compute the phase of all those points, then find the next
    % binCenterNode and
    % do it again.
    %
    binNum = 1;
    thisNodeIndex = cortMag.nodeIndices{roiNum};
    % Here we use a sparse matix to create a look-up table of phase and co values
    % for each nodeIndex. These data have been derived from the flat map
    % coranal. The corMag struct has a flag (flatDataFlag) that tells us if we should use
    % the flat data or actually go and get the data from the volume.
    nodePhTable = sparse(cortMag.nodeIndices{roiNum}, 1, cortMag.data{roiNum}.ph);
    nodeCoTable = sparse(cortMag.nodeIndices{roiNum}, 1, cortMag.data{roiNum}.co);
    while(~isempty(thisNodeIndex))
        % Select binCenterNode
        bins{roiNum}(binNum).binCenterNode = thisNodeIndex(1);
        if(length(thisNodeIndex)>1)
            thisNodeIndex = thisNodeIndex(2:end);

            % Find all points within binDist of binCenterNode.
            % We actually don't care about the distances right now, so we
            % just grab the
            % indices of the non-NaN entries.
            nodesInBin = find(~isnan(mrManDist(nodes, edges, bins{roiNum}(binNum).binCenterNode, ...
                                            cortMag.mmPerPix, NaN, cortMag.binDist)));
            [c, ia] = intersect(thisNodeIndex, nodesInBin);
            bins{roiNum}(binNum).allNodes = [bins{roiNum}(binNum).binCenterNode, c];
            % now we can remove all the nodes that intersected (ia) from
            % thisNodeIndex.
            % BUT- if ia is not consecutive- that's a problem. These nodes
            % are supposed to be sorted for
            % us! We do a hack for now.
            if(~isempty(ia))
                thisNodeIndex = thisNodeIndex(max(ia):end);
            end

            if(isfield(cortMag,'flatDataFlag') & cortMag.flatDataFlag==1)
                bins{roiNum}(binNum).allPh = full(nodePhTable(bins{roiNum}(binNum).allNodes));
                bins{roiNum}(binNum).allCo = full(nodeCoTable(bins{roiNum}(binNum).allNodes));
            else

                bins{roiNum}(binNum).allPh = grayView.ph{cortMag.expNumber}(bins{roiNum}(binNum).allNodes);
                bins{roiNum}(binNum).allCo = grayView.co{cortMag.expNumber}(bins{roiNum}(binNum).allNodes);
            end

            % But, the clipped nodes may contain some unconnected nodes.
            % So, we may need to go back to
            % the complete set of nodes. And, since those have different
            % indices, we need to extract the
            % volume coords and then look up the corresponding data.
            % Unfortunately, this makes this loop an order of magnitude
            % slower. So, we cheat with the above shortcut.
            % The following code is a stream-lined version of
            % 'getCurDataROI' to try to speed things up a bit.
%             coords = nodes([2,1,3], bins{roiNum}(binNum).allNodes);
%             [inter,ROIIndices,viewIndices] = intersectCols(coords,
%             grayView.coords);
%             bins{roiNum}(binNum).allPh = NaN*ones([1,size(coords,2)]);
%             bins{roiNum}(binNum).allPh(ROIIndices) =
%             grayView.co{cortMag.expNumber}(viewIndices);
%             bins{roiNum}(binNum).allCo = NaN*ones([1,size(coords,2)]);
%             bins{roiNum}(binNum).allCo(ROIIndices) =
%             grayView.ph{cortMag.expNumber}(viewIndices);

            binNum = binNum+1;
        else
            thisNodeIndex = [];
        end
    end
    bins{roiNum}(1).distToPrev = 0;
    for(binNum=2:length(bins{roiNum}))
        allDist = mrManDist(nodes, edges, bins{roiNum}(binNum).binCenterNode, ...
                            cortMag.mmPerPix, NaN, 0);
        bins{roiNum}(binNum).distToPrev = allDist(bins{roiNum}(binNum-1).binCenterNode);
    end
    if(any(isnan(cat(1,bins{roiNum}.distToPrev))))
        % Oops! this shouldn't happen with a properly connected mesh!
        bins{roiNum} = [];
        disp(['   ROI ',num2str(roiNum),' DROPPED BECAUSE OF UNCONNECTED NODES.']);
    else
        disp(['   ',num2str(binNum-1),' bins created from ',num2str(length(cortMag.nodeIndices{roiNum})),' nodes.']);
    end
end
cortMag.bins = bins;

return;

% Debugging
dist={}; meanPh={};
for(roiNum=1:length(cortMag.ROIs))
    for(binNum=1:length(bins{roiNum}))
        dist{roiNum}(binNum) = bins{roiNum}(binNum).distToPrev;
        meanPh{roiNum}(binNum) = mean(unwrap(bins{roiNum}(binNum).allPh));
    end
    figure;plot(cumsum(dist{roiNum}), meanPh{roiNum});
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% LFDefaultField - Convenience function to set up structs with default field values
%
% Usage:
%
%   ParentStruct = LFDefaultField( ParentStruct, FieldName, DefaultVal )
%
% This provides an elegant way to establish default field values in a struct. See LFDefaultValue for
% setting up non-struct variables with default values.
%
% Inputs:
%
%   ParentStruct: string giving the name of the struct; the struct need not already exist
%   FieldName: string giving the name of the field
%   DefaultVal: default value for the field
%
% Outputs:
%
%   ParentStruct: if the named struct and field already existed, the output matches the struct's
%                 original value; otherwise the output is a struct with an additional field taking
%                 on the specified default value
%
% Example:
%
%   clearvars
%   ExistingStruct.ExistingField = 42;
%   ExistingStruct = LFDefaultField( 'ExistingStruct', 'ExistingField', 3 )
%   ExistingStruct = LFDefaultField( 'ExistingStruct', 'NewField', 36 )
%   OtherStruct = LFDefaultField( 'OtherStruct', 'Cheese', 'Indeed' )
%
%
%   Results in :
%       ExistingStruct =
%           ExistingField: 42
%                NewField: 36
%       OtherStruct =
%                  Cheese: 'Indeed'
%
% Usage for setting up default function arguments is demonstrated in most of the LF Toolbox
% functions.
%
% User guide: LFToolbox.pdf
% See also: LFDefaultVal

% Copyright (c) 2013-2020 Donald G. Dansereau

function ParentStruct = LFDefaultField( ParentStruct, FieldName, DefaultVal )

%---First make sure the struct exists---
CheckIfExists = sprintf('exist(''%s'', ''var'') && ~isempty(%s)', ParentStruct, ParentStruct);
VarExists = evalin( 'caller', CheckIfExists );
if( ~VarExists )
    ParentStruct = [];
else
    ParentStruct = evalin( 'caller', ParentStruct );
end

%---Now make sure the field exists---
if( ~isfield( ParentStruct, FieldName ) )
    ParentStruct.(FieldName) = DefaultVal;
end

end

repo: LFToolbox


text:
%% Axis
% Enumeration class for the Cartesian axes x, y, z.

%%% Description
% The only three instances |Axis.x|, |Axis.y|, and |Axis.z| of |Axis| are used
% extensively in MaxwellFDFD to represent and specify the x-, y-, and z-axis of
% the Cartesian coordinate system.
%
% |Axis| is a subclass of |Enumerated|, so it inherits all the methods of
% |Enumerated|.  See <Enumerated.html |Enumerated|> for more details.

%%% Instances
% * |Axis.x|: instance representing the x-axis
% * |Axis.y|: instance representing the y-axis
% * |Axis.z|: instance representing the z-axis

%%% Methods
% Below, |axis| is an instance of |Axis|.
%
% * |[p, q, r] = cycle(axis)|: |[p, q, r]| is a cyclic permutation of |[Axis.x,
% Axis.y, Axis.z]| satisfying |r == axis|.

%%% Example
%   % Test conversion to integers and strings.
%   fprintf('# of instances of Axis: %d\n', Axis.count);
%   for w = Axis.elems
%       fprintf('The integer value of Axis.%s is %d\n', char(w), int(w));
%   end
%
%   % Test cyclic permutation.
%   [p q r] = cycle(Axis.y);
%   fprintf('The cyclic permutation of [x, y, z] beginning with y is [%s, %s, %s]\n', char(r), char(p), char(q));

%%% See Also
% <Enumerated.html |Enumerated|>

classdef Axis < Enumerated
	enumeration
		x('x')
		y('y')
		z('z')
	end

	methods (Static)
		function elems = elems(ind)
			elems = [Axis.x, Axis.y, Axis.z];
			if nargin > 0  % ind
				elems = elems(ind);
			end
		end

		function count = count()
			count = length(Axis.elems);
		end
	end

	methods
		function [p q r] = cycle(this)
			r = this;
			p = Axis.elems(mod(int(this), Axis.count) + 1);
			q = Axis.elems(mod(int(this)+1, Axis.count) + 1);
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function rasterSegs = structFillin(rasterSegs1, scanNum, planC)
%"structFillin"
%   Take the rasterSegments from one structure and fill in all slices with
%   no defined structure.  Return the rasterSegments of the new filled in
%   structure.  Uses nearest neighbor interpolation to duplicate the
%   contours.
%
%   By JRA 1/7/05
%
%   rasterSegs1    : rasterSegments of a structure
%   planC          : CERR planC
%
%   rasterSegs     : rasterSegments of filled in structure
%
%Usage:
%   rasterSegs = structFillin(rasterSegs1, scanNum, planC)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};
rasterSegs = [];
[xV, yV, zV] = getScanXYZVals(planC{indexS.scan}(scanNum));

%sort input rasterSegments by CTSliceValue
rasterSegs1 = sortrows(rasterSegs1, 6);

%get list of CTSlices to iterate over.
slices1 = unique(rasterSegs1(:,6));
minSlice = min(slices1);
maxSlice = max(slices1);

%Iterate over all slices.  If no rasterSegs exist on a slice, find the
%nearest slice, take its rasterSegs, modify them to fit the new
%slice/zvalue, and add them to the rasterSegs list.
for sliceNum = minSlice:maxSlice
    rasterIndices = find(rasterSegs1(:,6) == sliceNum);
    if isempty(rasterIndices)
        sliceZValue = zV(sliceNum);
        nearestDefinedSliceNum = interp1(zV(slices1), slices1, sliceZValue, 'nearest');
        tmpIndV = find(rasterSegs1(:,6) == nearestDefinedSliceNum);
        tmpSegsV = rasterSegs1(tmpIndV,:);
        tmpSegsV(:,1) = zV(sliceNum);
        tmpSegsV(:,6) = sliceNum;
        rasterSegs = [rasterSegs;tmpSegsV];
    else
        rasterSegs = [rasterSegs;rasterSegs1(rasterIndices,:)];
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function pathstr = getPathStrDialog(startPathStr,title,filterspec)
%
% fullPath = getPathStrDialog(startPathStr,title,[filterspec])
%
% Calls uigetfile to display a dialog box that allows the user to pick a file.
%
% startPathStr: starting point for the dialogbox
% title: window title of dialog box, passed on to uigetfile.
% filterspec: passed onto uigetfile. default: '*'
%
% returns full path string
%
% djh, 2/5/2001

if ~exist('filterspec','var')
    filterspec = '*';
end
if ~exist(startPathStr,'dir')
    startPathStr = pwd;
end

% save current file position
cDir=pwd;
% move to appropriate directory
chdir(startPathStr);
% call uigetfile
[filename, pathStr]=uigetfile(filterspec,title);
% move back to original directory
chdir(cDir);

if (filename==0) % cancel pressed
    pathstr = [];
else
    pathstr=fullfile(pathStr,filename);
end


text:
% Change to your downloaded location
clear
addpath('C:\liblinear\matlab')
addpath('../training_code/');
addpath('../utilities/');
addpath('../../data extraction/');
%% load shared definitions and AU data
shared_defs;

% Set up the hyperparameters to be validated
hyperparams.c = 10.^(-7:0.5:1);
hyperparams.e = 10.^(-3);

hyperparams.validate_params = {'c', 'e'};

% Set the training function
svm_train = @svm_train_linear;

% Set the test function (the first output will be used for validation)
svm_test = @svm_test_linear;

pca_loc = '../../pca_generation/generic_face_rigid.mat';

hog_data_dir_BP4D = hog_data_dir;

aus = [1, 2, 4, 6, 7, 10, 12, 14, 15, 17, 23];
%%
for a=1:numel(aus)

    au = aus(a);

    rest_aus = setdiff(all_aus, au);

    % load the training and testing data for the current fold
    [train_samples, train_labels, valid_samples, valid_labels, ~, PC, means, scaling] = Prepare_HOG_AU_data_generic(train_recs, devel_recs, au, BP4D_dir, hog_data_dir_BP4D);

    train_samples = sparse(train_samples);
    valid_samples = sparse(valid_samples);

    %% Cross-validate here
    [ best_params, ~ ] = validate_grid_search_no_par(svm_train, svm_test, false, train_samples, train_labels, valid_samples, valid_labels, hyperparams);
    model = svm_train(train_labels, train_samples, best_params);

    [~, predictions_all] = svm_test(valid_labels, valid_samples, model);

    name = sprintf('results_BP4D_devel/AU_%d_static.mat', au);

    [ accuracies, F1s, corrs, ccc, rms, classes ] = evaluate_regression_results( predictions_all, valid_labels );

    save(name, 'model', 'F1s', 'accuracies', 'predictions_all', 'valid_labels');

    % Write out the model
    name = sprintf('models/AU_%d_static.dat', au);

    pos_lbl = model.Label(1);
    neg_lbl = model.Label(2);

    w = model.w(1:end-1)';
    b = model.w(end);

    svs = bsxfun(@times, PC, 1./scaling') * w;

    write_lin_svm(name, means, svs, b, pos_lbl, neg_lbl);
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function [tstop,res] = testProblem(optObj,solvers)
%TESTPROBLEM  Test an OPTI Problem against all available solvers
%
%   [times,results] = testProblem(optObj) runs the supplied OPTI problem
%   against all available solvers for the given problem type.
%
%   [times,results] = testProblems(optObj,solvers) runs only the selected
%   solvers across the supplied problem.

%   Copyright (C) 2011 Jonathan Currie (I2C2)

%Setup default arguments
if(~exist('solvers','var') || isempty(solvers))
    solvers = checkSolver(['all_' optObj.prob.type]);
end
if(~iscell(solvers))
    solvers = {solvers};
end

%Remove troublesome solvers
ind = strcmpi(solvers,'lipsol');
if(any(ind))
    solvers(ind) = [];
end

%Setup problem generation function (for warm up)
switch(lower(optObj.prob.type))
    case 'lp'
        testSet = @lp_prob;
    case 'milp'
        testSet = @milp_prob;
    case 'bilp'
        testSet = @bilp_prob;
    case 'qp'
        testSet = @qp_prob;
    case 'qcqp'
        testSet = @qcqp_prob;
    case 'sdp'
        testSet = @sdp_prob;
    case 'miqp'
        testSet = @miqp_prob;
    case 'nlp'
        testSet = @nlp_prob;
        if(isempty(optObj.prob.x0))
            error('Please include x0 in optiprob when using this function on NLPs');
        end
    case 'minlp'
        testSet = @minlp_prob;
        if(isempty(optObj.prob.x0))
            error('Please include x0 in optiprob when using this function on MINLPs');
        end
    otherwise
        error('Unknown problem type: %s',optObj.prob.type);
end

%Remove troublesome solvers
if(strcmpi(optObj.prob.type,'lp'))
    ind = strcmpi(solvers,'SEDUMI'); solvers(ind) = []; %no good with neq = nvar
end

%Setup no runs
len = length(solvers);
%confirm solver exists!
for i = 1:len
    rub = checkSolver(solvers{i}); %#ok<NASGU>
end

%Setup wait bar
h = waitbar(0,['Solver: ' upper(solvers{1})],'name','OPTI Problem Test','CreateCancelBtn','setappdata(gcbf,''canceling'',1)');

%Preallocate
tstop = zeros(len,1);
res = struct('sol',{cell(len,1)},'fval',zeros(len,1),'ok',zeros(len,1),'msg',{cell(len,1)});

%Clear display warnings etc
opts = optiset(optObj.prob,'display','off','warnings','off');

try
    for i = 1:len
        %Check Cancel
        if getappdata(h,'canceling')
            break;
        end
        %Assign the new solver
        solver = solvers{i};
        opts = optiset(opts,'solver',solver);
        %Initialize with two solves of problem 1
        prob = testSet(1);
        Opt = opti(prob,opts);
        solve(Opt); solve(Opt);
        %Now run our problem
        Opt = opti(getProb(optObj),opts);
        [x,sval,~,stat] = solve(Opt);
        res.sol{i} = x;
        res.fval(i) = sval;
        res.ok(i) = checkSol(Opt);
        res.msg{i} = stat.Status;
        tstop(i) = stat.Time;
        %Update Waitbar
        waitbar(i/len,h,['Solver: ' upper(solvers{min(len,i+1)})]);
    end
catch ME
    delete(h);
    rethrow(ME);
end

delete(h);

%Draw Table
fprintf('\nTest Problem [%s] Result\n',Opt.prob.type);
fprintf('--------------------------------------------------\n');
fprintf('  Solver    Status        Fval         Time\n');
for i = 1:len
    fprintf('%8s    %4s    %+12.4f     %1.4fs\n',solvers{i},status(res.ok(i)),res.fval(i),tstop(i));
end
fprintf('--------------------------------------------------\n');

function msg = status(ok)

if(ok)
    msg = 'OK';
else
    msg = 'FAIL';
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [K, dK_logtheta, dK_x2] = gpcovPeriodic(logtheta, x1, x2, distfun)

if nargin == 0
  % Return the number of hyperparameters
  K = 2;
  return
end

if nargin < 3
  x2 = [];
end
if nargin < 4
%  distfun = @sqdistEuclidean;
end

% 2-norm distances
if ~isempty(x2)
  D = bsxfun(@minus, x1(:), x2(:)');
  %dD_dx2 = reshape(-D, [1,length(x1),length(x2)]);
%  [D, dD_dx2] = feval(distfun, x1, x2, false);
else
  % Distances for variances
  n = cols(x1);
  D = zeros(n,1);
end

% TODO:
% helper, inverse squared length scale:
invl2 = exp(-2*logtheta(1));
wavel = exp(logtheta(2));

% Covariance matrix

angle = pi*D/wavel;
exponent = -2*sin(angle).^2*invl2;
K = exp(exponent);

% Gradient for hyperparameters
if nargout >= 2
  dK_logtheta = zeros([size(D), 2]);
  dK_logtheta(:,:,1) = K .* exponent .* (-2);
  dK_logtheta(:,:,2) = K.*(-2*invl2*2*sin(angle)).*cos(angle).*(-pi*D/wavel);
end

% Gradients for inputs x2
if nargout >= 3
  if isempty(x2)
    error('Can''t calculate gradient: x2 not given');
  end
  d = rows(x2); % dimensionality of inputs
  m = cols(x1); % number of other inputs
  n = cols(x2); % number of inputs
  dK_x2 = zeros([d,m,n]);
  dK_x2 = bsxfun(@times, ...
                 reshape(K,[1,m,n]), ...
                 reshape((-2*invl2*2*sin(angle)) .* ...
                 cos(angle) .* pi/wavel .* (-1), [1,m,n]));
% $$$   for j=1:n
% $$$     dK_x2(:,:,j) = bsxfun(@times, K(:,j)', (-0.5*invl2) * dD2_dx2(:,:,j));
% $$$     %2*bsxfun(@minus, x2(:,j),x1));
% $$$   end
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function r = spm_invIcdf(F,n,p)
% Inverse Cumulative Distribution Function (CDF) of Binomial distribution
% FORMAT r = spm_invIcdf(F,n,p)
%
% F - CDF (lower tail p-value)
% n - Binomial n
% p - Binomial p [Defaults to 0.5]
% r - ordinate
%__________________________________________________________________________
%
% spm_invIcdf returns the inverse Cumulative Distribution Function for
% the Binomial family of distributions.
%
% Definition:
%--------------------------------------------------------------------------
% The Bin(n,p) distribution is the distribution of the number of
% successes from n identical independent Bernoulli trials each with
% success probability p. If random variable R is the number of
% successes from such a set of Bernoulli trials, then the CDF F(r) is
% the probability of r or less sucesses.
%
% The Binomial distribution is discrete, defined for p in [0,1] and r
% in {0,1,...,n}, so F(r) is a discrete function. This inverse CDF
% function returns the smallest Whole r such that the F(r) equals or
% exceeds the given CDF probability F. I.e. F(r) is treated as a step
% function.
%
% Algorithm:
%--------------------------------------------------------------------------
% r is found by direct summation of the Binomial PDFs until F is exceeded.
%
% References:
%--------------------------------------------------------------------------
% Evans M, Hastings N, Peacock B (1993)
%       "Statistical Distributions"
%        2nd Ed. Wiley, New York
%
% Abramowitz M, Stegun IA, (1964)
%       "Handbook of Mathematical Functions"
%        US Government Printing Office
%
% Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
%       "Numerical Recipes in C"
%        Cambridge
%
%__________________________________________________________________________
% Copyright (C) 1999-2011 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_invIcdf.m 4182 2011-02-01 12:29:09Z guillaume $

%-Format arguments, note & check sizes
%--------------------------------------------------------------------------
if nargin<3, p=0.5; end
if nargin<2, error('Insufficient arguments'), end
ad = [ndims(F);ndims(n);ndims(p)];
rd = max(ad);
as = [[size(F),ones(1,rd-ad(1))];...
      [size(n),ones(1,rd-ad(2))];...
      [size(p),ones(1,rd-ad(3))]];
rs = max(as);
xa = prod(as,2)>1;
if sum(xa)>1 && any(any(diff(as(xa,:)),1))
    error('non-scalar args must match in size');
end

%-Computation
%--------------------------------------------------------------------------
%-Initialise result to zeros
r = zeros(rs);

%-Only defined for whole n, p&F in [0,1]. Return NaN if undefined.
md = ( F>=0  &  F<=1  &  n==floor(n)  &  n>=0  &  p>=0  &  p<=1 );
if any(~md(:))
    r(~md) = NaN;
    warning('Returning NaN for out of range arguments');
end

%-Non-zero only where defined, & F>0
Q  = find( md  &  F>0 );
if isempty(Q), return, end
if xa(1), QF=Q; else QF=1; end
if xa(2), Qn=Q; else Qn=1; end
if xa(3), Qp=Q; else Qp=1; end

%-Compute by directly summing Bin PDF's for successive r & comparing with F
tr  = 0;
Ftr = spm_Ipdf(tr,n(Qn),p(Qp));
while any(F(QF)>Ftr) && any(n(Qn)>tr)
    tr      = tr+1;
    i       = find(Ftr<F(QF));
    r(Q(i)) = r(Q(i)) + 1;
    Ftr     = Ftr + spm_Ipdf(tr,n(Qn),p(Qp));
end

author: spm


text:
%% MULTIGRID OF FOR THE STOKES EQNS IN 2D
%
% This example is to show the convergence of multigrid methods for various
% finite element approximation of the Stokes equation on the unit square:
%
%       -div(mu*grad u) + grad p = f in \Omega,
%                        - div u = 0  in \Omega,
%                              u = g_D  on \Gamma.
%
% with the pure Dirichlet boundary condition.
%

%% Setting
clear all; close all;
[node,elem] = squaremesh([0,1,0,1],0.5);
% [node,elem] = circlemesh(0,0,1,0.25);
[node,elem] = uniformrefine(node,elem);
showmesh(node,elem);
pde = Stokesdata1;
%  pde = StokesZulehnerdata;
bdFlag = setboundary(node,elem,'Dirichlet');
option.L0 = 1;
option.maxIt = 4;
option.printlevel = 1;
option.plotflag = 0;
option.rateflag = 0;

%% MG options
option.solver = 'mg';
option.smoothingstep = 2;
option.smootherbarSp = 'SGS';

%% RT0-P0
display('RT0-P0')
option.elemType = 'RT0-P0';
option.refType = 'bisect';
femStokesHdiv(node,elem,pde,bdFlag,option);

% %% RT0-P0
% display('BDM1B-P0')
% option.elemType = 'BDM1B-P0';
% femStokesHdiv(node,elem,pde,bdFlag,option);

%% CR-P0 element
display('CR-P0')
option.elemType = 'CRP0';
femStokes(node,elem,pde,bdFlag,option);

%% P2-P0 element
display('P2-P0')
option.elemType = 'P2P0';
femStokes(node,elem,pde,bdFlag,option);

%% isoP2-P0 element
display('isoP2-P0')
option.elemType = 'isoP2P0';
femStokes(node,elem,pde,bdFlag,option);

%% isoP2-P1 element
display('isoP2-P1')
option.elemType = 'isoP2P1';
femStokes(node,elem,pde,bdFlag,option);

% %% P1b-P1 element
% display('P1b-P0')
% option.elemType = 'P1bP1';
% femStokes(node,elem,pde,bdFlag,option);

%% P2-P1 element
display('P2-P1')
option.elemType = 'P2P1';
% option.smoothingStep = 3;
% option.smootherbarSp   = 'VCYCLE';
% option.smootherbarSpPara = 0.75;
femStokes(node,elem,pde,bdFlag,option);

repo: Aero-matlab


text:
function [I_b_deep,I_d_deep, show_matrix, count_all, I_vector] = vector_deep_latent(img, de_level, L, unit, is_overlap, stride, w_num)

[t1,t2] = size(img);
show_matrix = [];
% I_d = zeros(t1,t2,de_level);
I_b_deep = zeros(t1,t2);
for i=1:de_level
    if i==1
        [temp_b, temp_d, temp_d_v, count_m, I_vector] = vector_decomposition(img, L, unit, is_overlap, stride, w_num);
    else
        [temp_b, temp_d, temp_d_v, count_m, I_vector] = vector_decomposition(temp_b, L, unit, is_overlap, stride, w_num);
    end
    I_d_deep(:,:,i) = temp_d_v(:,:);
    show_matrix = cat(2,show_matrix, temp_d);

end
% show_matrix = cat(2,show_matrix, temp_b);
% figure;imshow(show_matrix);
% imwrite(temp_b,['./fused_vector/features/base_level_',num2str(i),'.png'],'png');

I_b_deep(:,:) = temp_b(:,:);
count_all = count_m;
% I_d_deep = I_d;

end

author: thfylsty


text:
function Test_symYY()
% function test98()
% Test for symfixedrankYYquotientfactory geometry (low-rank PSD matrix completion)
%
% Paper link: http://www.di.ens.fr/~fbach/journee2010_sdp.pdf
%

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, July 11, 2013.
% Contributors:
% Change log:

% We know about this warning, so it is safe to turn it off for now.
warning('off', 'manopt:symfixedrankYYfactory:exp');

clc; close all;

% Problem data
n = 1000;
r = 5;
Y_org = randn(n, r);
A = Y_org*Y_org';

% Create the problem structure
% quotient YY geometry
problem.M = symfixedrankYYfactory(n, r);

df = problem.M.dim();
p = 3.5*df/(n*n);
mask = spones(sprandsym(n, p));
prop_known = sum(sum(mask == 1))/(n*n);
fprintf('Fraction of entries given: %f \n', full(prop_known));

problem.cost = @(Y) norm(mask.*(Y*Y' - A),'fro')^2;

egrad = @(Y) 4*((mask.^2).*(Y*Y' - A))*Y;
ehess = @(Y, U) 4*((mask.^2) .* (Y*Y' - A))*U + 4*((mask.^2) .* (Y*U' + U*Y'))*Y;

problem.grad = @(Y) problem.M.egrad2rgrad(Y, egrad(Y));
problem.hess = @(Y, U) problem.M.ehess2rhess(Y, egrad(Y), ehess(Y, U), U);

% % Check numerically whether gradient and Hessian are correct
%     checkgradient(problem);
%     drawnow;
%     pause;
%     checkhessian(problem);
%     drawnow;
%     pause;

% Initialization
[U0, S0, ~ ] = svds(mask.*A, r);
Y0 = U0*(S0.^0.5);

% Options (not mandatory)
options.maxiter = inf;
options.maxinner = 30;
options.maxtime = 120;
options.tolgradnorm = 1e-5;
options.Delta_bar = n * r;
options.Delta0 = options.Delta_bar / 8;

% Pick an algorithm to solve the problem
[Yopt costopt info] = trustregions(problem, Y0, options);
% [Yopt costopt info] = steepestdescent(problem, Y0, options);
% [Yopt costopt info] = conjugategradient(problem, Y0, options);

end

author: NicolasBoumal

repo: manopt


text:
function CPD = root_gaussian_CPD(bnet, self, mu, Sigma, mu0, n0, alpha0, beta0)
% ROOT_GAUSSIAN_CPD Make an unconditional Gaussian distrib.
%
% CPD = root_gaussian_CPD(bnet, self, mu, Sigma)
% This defines the distribution Y ~ N(mu, Sigma),
% Pass in [] to generate a default random value for a parameter.
%
% CPD = root_gaussian_CPD(bnet, self, [], [], mu0, n0, alpha0, beta0)
% defines a Normal-Wishart prior over the parameters:
%   P(mu | lambda) = N(mu | mu0, n0*lambda)
%   P(lambda) = Wishart(lambda | alpha0, beta0)
% where lambda = inv(Sigma) is the precision matrix of mu.
% n0 is a scale factor, beta0 is a precision matrix.
% Pass in [] to generate a default value for a hyperparameter.
% mu and Sigma will be set to their prior expected values.
% See "Bayesian Theory", Bernardo and Smith (2000), p441.

if nargin==0
  % This occurs if we are trying to load an object from a file.
  CPD = init_fields;
  CPD = class(CPD, 'root_gaussian_CPD', generic_CPD(0));
  return;
elseif isa(bnet, 'root_gaussian_CPD')
  % This might occur if we are copying an object.
  CPD = bnet;
  return;
end
CPD = init_fields;

ns = bnet.node_sizes;
d = ns(self);

if nargin < 5,
  prior = [];
  if isempty(mu), mu = randn(d, 1); end
  if isempty(Sigma), Sigma = eye(d); end
else
  if isempty(mu0), mu0 = zeros(d, 1); end
  if isempty(n0), n0 = 0.1; end
  if isempty(alpha0), alpha0 = (d-1)/2 + 1; end % Wishart requires 2 alpha > d-1
  if isempty(beta0), beta0 = eye(d); end

  prior.mu = mu0;
  prior.n = n0;
  prior.alpha = alpha0;
  prior.beta = beta0;

  % set params to their mean
  mu = prior.mu;
  Sigma = prior.beta/prior.alpha; % mean of Wishart is E[lambda] = alpha*inv(beta)
end

CPD.self = self;
CPD.mu = mu;
CPD.Sigma = Sigma;
CPD.prior = prior;

clamped = 0;
CPD = class(CPD, 'root_gaussian_CPD', generic_CPD(clamped));

%%%%%%%%%%%

function CPD = init_fields()
% This ensures we define the fields in the same order
% no matter whether we load an object from a file,
% or create it from scratch. (Matlab requires this.)

CPD.self = [];
CPD.mu = [];
CPD.Sigma = [];
CPD.prior = [];

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%MAPEX Train and execute arbitrary untrained mapping
%
%    B = A*(W*MAPEX) = A*(A*W)
%    B = A*MAPEX(W,PAR1,PAR2, ...)
%    B = A*MAPEX(UMAP,PAR1,PAR2, ...)
%
% INPUT
%    A        Dataset or datafile or double
%    W        Untrained mapping
%    UMAP     String with name of untrained mapping
%    PAR1     Parameters of untrained mapping W or UMAP
%
% OUTPUT
%    B               Resulting dataset, datafile or double array
%
% DESCRIPTION
% This routine facilitates the construction of shortcuts by training and
% executing an untrained mapping by the same data. It is typically useful
% for mappings like PROXM and SCALEM that are often trained and exectuted
% with the same data.
%
% EXAMPLE
% mink1 = proxm('m',1)*mapex
% mink1 = mapex(proxm,'m',1)    % the same
% mink1 = mapex('proxm','m',1)  % the same
% Herewith D = A*mink1 computes a dissimilarity matrix based on the
% Minkowsky_1 metric between all objects in A.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, PROXM, SCALEM

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

function out = mapex(p1,varargin)

  if nargin == 0
    out = prmapping(mfilename,'combiner');
  elseif nargin == 1 && ismapping(p1)
    out = prmapping(mfilename,'fixed',p1);
  elseif nargin >= 1 && ischar(p1)
    out = prmapping(mfilename,'fixed',{p1,varargin{:}});
  elseif nargin > 1 && ismapping(p1)
    out = prmapping(mfilename,'fixed',{getmapping_file(p1),varargin{:}});
  elseif (isdataset(p1) || isa(p1,'double')) && ismapping(varargin{1})
    out = p1*(p1*varargin{1});
  elseif (isdataset(p1) || isa(p1,'double')) && nargin == 2
    out = p1*feval(varargin{1},p1);
  elseif (isdataset(p1) || isa(p1,'double')) && nargin > 1
    out = p1*feval(varargin{1},p1,varargin{2:end});
  else
    error('Illegal input')
  end

author: marianux

repo: ecg-kit

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


text:
function counts = numNeighbors(grains)
% returns the number of neighboring grains
%
% Input
%  grains - @grain2d
%
% Output
%  counts - number of neighbors per grain
%

% get list of neighbouring grains
pairs = grains.neighbors('full');

% get the number of neighbours per grain
ng = max(grains.id);
countIds = full(sparse(pairs(pairs<=ng),1,1,ng,1));

% rearange with respect to grain order
counts = countIds(grains.id);

author: mtex-toolbox

repo: mtex

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


text:
function [ fn ] = getBatchWrapper_ucf101_flow(opts, numThreads, trainopts)

fn = @(imdb,batch,valopts,frameList) getBatch(imdb,batch,opts,numThreads, trainopts, valopts) ;
end

% -------------------------------------------------------------------------
function [im,labels] = getBatch(imdb, batch, opts, numThreads, trainopts, moreopts)
% -------------------------------------------------------------------------
opts.nFrames = imdb.images.nFrames(batch);

images = cell(1,numel(opts.nFrames));
labels = imdb.images.label(batch) ;
for k = 1:numel(opts.nFrames)
  images{k}  = strcat([imdb.images.name{batch(k)}(1:end-4) filesep], images{k} ) ;
end
if ~isempty(moreopts) || nargin < 6
  for f = fieldnames(moreopts)'
    f = char(f) ;
    trainopts.(f) = moreopts.(f);
  end
end

opts.nFramesPerVid = trainopts.nFramesPerVid;
opts.numAugments = trainopts.numAugments;
opts.frameSample = trainopts.frameSample;
opts.augmentation =  trainopts.augmentation;

if isfield(trainopts, 'temporalStride')
  opts.temporalStride =  trainopts.temporalStride(randi(numel(trainopts.temporalStride))) ; % shuffle;
end
if isfield(trainopts, 'nFrameStack')
  opts.nFrameStack = trainopts.nFrameStack;
else
  opts.nFrameStack = trainopts.nFramesPerVid;
end
if isfield(trainopts, 'keepFramesDim')
  keepFramesDim = trainopts.keepFramesDim;
else
  keepFramesDim = false;
end
if isfield(trainopts, 'cheapResize')
  opts.cheapResize =  trainopts.cheapResize;
end

if isfield(trainopts, 'prefetch')
    opts.prefetch = trainopts.prefetch;
end
if isfield(trainopts, 'frameList')
    opts.frameList = trainopts.frameList;
end

im = cnn_ucf101_get_flow_batch(images, opts, ...
                            'numThreads', numThreads, ...
                            'augmentation', opts.augmentation , ...
                          'imageDir', imdb.flowDir) ;
if iscell(im), return; end

if nargout == 1 % to work with dag code
  if ndims(im) > 4 && ~keepFramesDim
    sz = size(im);
    nFrames = sz(5:end);
    im = permute(im, [1 2 3 5 4]);
    im = reshape(im, sz(1), sz(2), sz(3), []);
  end
  im = {'input', im, 'label', labels } ;
end

end

author: feichtenhofer

repo: st-resnet


text:
% VL_NNROIPOOL  CNN region of interest pooling.
%   Y = VL_NNROIPOOL(X, ROIS) pools each feature channel in X in
%   the specified regions of interest ROIS. ROIS is a 5 x K array
%   containing K regions. Each region has five coordinates `[t, u0,
%   v0, u1, v1]` where `u0`, `v0` is the upper-left corner of a ROI,
%   `u1`, `v1` is the bottom-right corner, and `t` is the index of the
%   image that contains the region. Spatial coordiantes start at (1,1),
%   with `u` indexing the horizontal axis and `v` the vertical one.
%   The image indeces ranges from 1 to the number of images stored
%   in the tensor X.
%
%   If X has C feature channels, then the output Y is a 1 x 1 x C x K
%   array, with one image instance per region. Arguments can be SINGLE
%   or DOUBLE and CPU or GPU arrays; however, they must all be of the
%   same type (unless empty).
%
%   DZDX = VL_NNROIPOOL(X, ROIS, DZDY) computes the derivative of
%   the layer projected on DZDY with respect to X.
%
%   VL_NNROIPOOL(___, 'opt', value, ...) accepts the following
%   options:
%
%   `Method`:: `'max'`
%     Choose between `'max'` and `'avg'` (average) pooling.
%
%   `Subdivisions`:: `[1 1]`
%     Specifies the number [SH,SW] of vertical and horizontal tiles of
%     a region. This makes the output a SH x SW x C x K array.
%
%   `Transform`:: `1`
%     Specifies a spatial transformation to apply to region vertices before
%     they are applied to the input tensor. If T is a scalar, then
%     the transformation is a scaling centered at the origin:
%
%        u' = T (u - 1) + 1,
%        v' = T (v - 1) + 1.
%
%     If T is a 2D vector, then different scaling factors for the
%     `u` and `v` can be specified. Finally, if T is a 2 x 2 matrix, then:
%
%        u' = T(1,1) u + T(1,2) v + T(1,3),
%        v' = T(2,1) u + T(2,2) v + T(2,3).
%
%     Note that only the upper-left and bottom-right corners of each
%     rectangular region are transformed. Thus this is mostly useful
%     for axis-aligned transformations; the generality of the expression
%     allows, however, to swap `u` and `v`, which may be needed
%     to match different conventions for the box coordiantes.
%
%   See also: VL_NNPOOL().

% Copyright (C) 2016 Hakan Bilen, Abishek Dutta, and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function timing_MX
% TIMING_MX  Speed of MX as performed by MULTIPROD and by a nested loop.
%    TIMING_MX compares the speed of matrix expansion as performed by
%    MULTIPROD and an equivalent nested loop. The results are shown in the
%    manual (fig. 2).
%    Notice that MULTIPROD enables array expansion which generalizes matrix
%    expansion to arrays of any size, while the loop tested in this
%    function works only for this specific case, and would be much slower
%    if it were generalized to N-D arrays.

% Checking whether needed software exists
message = sysrequirements_for_testing('timeit');
if message
    disp ' ', error('testing_memory_usage:Missing_subfuncs', message)
end

% Matrix expansion example (fig. 2)
disp ' '
disp 'Timing matrix expansion (see MULTIPROD manual, figure 2)'
disp ' '

a = rand(2, 5);
b = rand(5, 3, 1000, 10);

fprintf ('Size of A:  %0.0fx%0.0f\n', size(a))
fprintf ('Size of B: (%0.0fx%0.0f)x%0.0fx%0.0f\n', size(b))
disp ' ', disp 'Please wait...'
disp ' '

f1 = @() loop(a,b);
f2 = @() multiprod(a,b);

t1 = timeit(f1)*1000;
fprintf('LOOP(A, B):      %10.4f milliseconds\n', t1)
t2 = timeit(f2)*1000;
fprintf('MULTIPROD(A, B): %10.4f milliseconds\n', t2)

disp ' '
fprintf('MULTIPROD performed matrix expansion %6.0f times faster than a plain loop\n', t1/t2)
disp ' '

function C = loop(A,B)
for i = 1:1000
    for j = 1:10
        C(:,:,i,j) = A * B(:,:,i,j);
    end
end


% params
totalNum = 1502861;
bValue = 2000;
lamda = 0.1;
qFile = '/home/yinhuan/Data/mapModel/park/weightVector.txt';
visFilesDir = '/home/yinhuan/Data/mapModel/park/visMatrix/';
maxQ = 110;
minQ = 1;

splitLength_0 = 100;
saveResultsDir_0 = '/home/yinhuan/Data/mapModel/park/iter_b_2000/0/';
saveReIndexFile_0 = '/hozcme/yinhuan/Data/mapModel/park/iter_b_2000/index/0.txt';
saveNewQFile_0 = '/home/yinhuan/Data/mapModel/park/iter_b_2000/weight/0.txt';

splitLength_1 = 200;
saveResultsDir_1 = '/home/yinhuan/Data/mapModel/park/iter_b_2000/1/';
saveReIndexFile_1 = '/home/yinhuan/Data/mapModel/park/iter_b_2000/index/1.txt';
saveNewQFile_1 = '/home/yinhuan/Data/mapModel/park/iter_b_2000/weight/1.txt';

splitLength_2 = 400;
saveResultsDir_2 = '/home/yinhuan/Data/mapModel/park/iter_b_2000/2/';
saveReIndexFile_2 = '/home/yinhuan/Data/mapModel/park/iter_b_2000/index/2.txt';
saveNewQFile_2 = '/home/yinhuan/Data/mapModel/park/iter_b_2000/weight/2.txt';

splitLength_3 = 800;
saveResultsDir_3 = '/home/yinhuan/Data/mapModel/park/iter_b_2000/3/';
saveReIndexFile_3 = '/home/yinhuan/Data/mapModel/park/iter_b_2000/index/3.txt';
saveNewQFile_3 = '/home/yinhuan/Data/mapModel/park/iter_b_2000/weight/3.txt';

splitLength_4 = 800;
saveResultsDir_4 = '/home/yinhuan/Data/mapModel/park/iter_b_2000/4/';
saveReIndexFile_4 = '/home/yinhuan/Data/mapModel/park/iter_b_2000/index/4.txt';
saveNewQFile_4 = '/home/yinhuan/Data/mapModel/park/iter_b_2000/weight/4.txt';

%% iteration of optimization

visCells_origin = fromVisDirToCells(visFilesDir);

[ epsilon_soft_0, time_sum_0  ] = Uniform_loopCompress_Cells( lamda, qFile, visCells_origin, splitLength_0, totalNum, bValue, saveResultsDir_0 );
min_cost_0 = get_min_cost(minQ, maxQ, saveResultsDir_0, qFile, lamda, epsilon_soft_0);
[totalNum_0, visCells_0] = deleteZeroPoints_Cells( qFile, visCells_origin, saveResultsDir_0, saveReIndexFile_0, saveNewQFile_0 );

[ epsilon_soft_1, time_sum_1  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_0, visCells_0, splitLength_1, totalNum_0, bValue, saveResultsDir_1 );
min_cost_1 = get_min_cost(minQ, maxQ, saveResultsDir_1, saveNewQFile_0, lamda, epsilon_soft_1);
[totalNum_1, visCells_1] = deleteZeroPoints_Cells( saveNewQFile_0, visCells_0, saveResultsDir_1, saveReIndexFile_1, saveNewQFile_1 );

[ epsilon_soft_2, time_sum_2  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_1, visCells_1, splitLength_2, totalNum_1, bValue, saveResultsDir_2 );
min_cost_2 = get_min_cost(minQ, maxQ, saveResultsDir_2, saveNewQFile_1, lamda, epsilon_soft_2);
[totalNum_2, visCells_2] = deleteZeroPoints_Cells( saveNewQFile_1, visCells_1, saveResultsDir_2, saveReIndexFile_2, saveNewQFile_2 );

[ epsilon_soft_3, time_sum_3  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_2, visCells_2, splitLength_3, totalNum_2, bValue, saveResultsDir_3 );
min_cost_3 = get_min_cost(minQ, maxQ, saveResultsDir_3, saveNewQFile_2, lamda, epsilon_soft_3);
[totalNum_3, visCells_3] = deleteZeroPoints_Cells( saveNewQFile_2, visCells_2, saveResultsDir_3, saveReIndexFile_3, saveNewQFile_3 );

[ epsilon_soft_4, time_sum_4  ] = Uniform_loopCompress_Cells( lamda, saveNewQFile_3, visCells_3, splitLength_4, totalNum_3, bValue, saveResultsDir_4 );
min_cost_4 = get_min_cost(minQ, maxQ, saveResultsDir_4, saveNewQFile_3, lamda, epsilon_soft_4);
[totalNum_4, visCells_4] = deleteZeroPoints_Cells( saveNewQFile_3, visCells_3, saveResultsDir_4, saveReIndexFile_4, saveNewQFile_4 );

%% save the results

compressIndex_3 = anal_reindex_last(saveResultsDir_4, saveReIndexFile_3);
compressIndex_2 = anal_reindex_middle(compressIndex_3, saveReIndexFile_2);
compressIndex_1 = anal_reindex_middle(compressIndex_2, saveReIndexFile_1);
compressIndex_0 = anal_reindex_middle(compressIndex_1, saveReIndexFile_0);

dlmwrite('/home/yinhuan/Data/mapModel/park/iter_b_2000/compressIndex_final.txt', compressIndex_0, 'precision', '%d');


text:
function [obj] = regObjective(obj, funcs)
    % This method registers the information of a NLP function as a
    % cost function
    %
    % Parameters:
    %  funcs: a list of objective functions @type NlpFunction
    %
    % Parameters:
    %  funcs: new NLP variables @type NlpFunction
    %
    % @note We always assume ''funcs'' to be a 2-D array of objects or
    % table consists of array objects. It also covers 1-D array or scalar
    % variables.
    %
    % @see NlpFunction

    if iscell(funcs)
        % validate the class of each objects
        assert(all(cellfun(@(x)isa(x,'NlpFunction'), funcs)), ...
            'Each variable must be an object of ''NlpFunction'' class or inherited subclasses');

        % remove empty (Dimension==0) objects from the the array
        funcs = funcs(~cellfun(@(x)x.Dimension==0,funcs));
        % contecate to the variable array
        obj.CostArray = [obj.CostArray; vertcat(funcs{:})];

    elseif isa(funcs, 'NlpFunction')
        % remove empty (Dimension==0) objects from the the array
        funcs = funcs(~arrayfun(@(x)x.Dimension==0,funcs));

        % contecate to the variable array
        obj.CostArray = [obj.CostArray; funcs(:)];
    elseif istable(funcs)

        % convert table content to an array
        funcs = transpose(table2array(funcs));
        % remove empty (Dimension==0) objects from the the array
        funcs = funcs(~arrayfun(@(x)x.Dimension==0,funcs));

        % contecate to the variable array
        obj.CostArray = [obj.CostArray; funcs(:)];

    else
        error('Unsupported variable type found: %s\n', class(funcs));
    end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function SDP = chordal_relax_mesh_registration(problem)
%% Apply a chordal sparse second-order relaxation to mesh registration
%% Generate SDP relaxations with multiple smaller blocks
%% This relaxation is likely to be looser than the single block relaxation
%% Depending on multivariate polynomial package SPOT
%% Heng Yang, July 28, 2021

fprintf('\n===================================================================')
fprintf('\nApplying Chordal SDP relaxation to mesh registration')
fprintf('\n===================================================================\n')
t0              = tic;

%% define POP variables
N                   = problem.N;
normalM             = problem.normalM;
normalP             = problem.normalP;
pointM              = problem.pointM;
pointP              = problem.pointP;
pointNoiseBoundSq   = problem.pointNoiseBoundSq;
normalNoiseBoundSq  = problem.normalNoiseBoundSq;
tBound              = problem.translationBound;
tBoundSq            = tBound^2; % t'*t <= tBoundSq
barc2               = 1.0;

nrPrimalVars    = 9 + 3 + N;
p               = msspoly('p',nrPrimalVars);
r               = p(1:9);
R               = reshape(r,3,3);
col1            = r(1:3);
col2            = r(4:6);
col3            = r(7:9);
t               = p(10:12);
theta           = p(13:nrPrimalVars);
x               = [r;t];

%% compute the cost function
residuals = [];
for i = 1:N
    pointM_i            = pointM(:,i);
    pointP_i            = pointP(:,i);
    normalM_i           = normalM(:,i);
    normalP_i           = normalP(:,i);
    residual_point_i    = ( normalM_i' * (R*pointP_i - pointM_i + t) )^2;
    residual_normal_i   = normalP_i'*normalP_i + normalM_i'*normalM_i - 2*normalP_i'*R'*normalM_i;
    residuals           = [residuals;...
                           (residual_point_i/pointNoiseBoundSq + residual_normal_i/normalNoiseBoundSq)/2];
end
f_cost = [];
for i = 1:N
    f_cost              = [f_cost;(1+theta(i))/2 * residuals(i) + (1-theta(i))/2 * barc2];
end

%% Define the equality and inequality constraints
h_x = [1.0-col1'*col1;...
        1.0-col2'*col2;...
        1.0-col3'*col3;... % columns unit length
        col1'*col2;...
        col2'*col3;...
        col3'*col1;... % colums orthogonal
        cross(col1,col2) - col3;...
        cross(col2,col3) - col1;...
        cross(col3,col1) - col2]; % columns righthandedness

h_theta = [];
for i = 1:N
    h_theta =[h_theta; 1-theta(i)^2];
end

g_x = tBoundSq - t'*t; % Translation bounded

%% Formulate the chordal sparse second-order relaxation
%% the 0-th block [1;x] * [1;x]'
basis0          = [1;x];
n0              = length(basis0);
basis_x0        = 1;

pop0            = [mykron(basis_x0,h_x);...
                    mykron(basis0,basis0)];
[~,degmat,coef_all] = decomp(pop0);
coef_all            = coef_all';
dim_loc0        = length(basis_x0) * length(h_x);
n0delta         = triangle_number(n0);
nterms          = size(degmat,1);
m_mom0          = n0delta - nterms;

assert(m_mom0==0,'The zero-th blk should have 0 moment constraints.')

coef_mom    = coef_all(:,dim_loc0+1:end);
coef_mom    = coef_mom';
B           = {};
B_normalize = {};

for i = 1:nterms
    [row,~,~]   = find(coef_mom(:,i));
    SDP_coli    = floor((row-1)./n0) + 1;
    SDP_rowi    = mod(row-1,n0) + 1;
    nnz         = length(SDP_rowi);

    Bi          = sparse(SDP_rowi,SDP_coli,ones(nnz,1),n0,n0);
    B{end+1}    = Bi;
    B_normalize{end+1} = Bi/nnz;
end

coef_loc0       = coef_all(:,1:dim_loc0);
A0_local        = {};

for i = 1:dim_loc0
    [rowi,~,vi] = find(coef_loc0(:,i));
    Ai      = sparse(n0,n0);
    for j   = 1:length(rowi)
        Ai  = Ai + vi(j) * B_normalize{rowi(j)};
    end
    A0_local = [A0_local;{Ai}];
end
A0_0    = sparse([1],[1],[1],n0,n0);
% The first block satisfies A0(X0) = b0;
A0      = [{A0_0};A0_local];
b0      = sparse(1,1,1,length(A0),1);

%% the 1-N blocks [1;x;theta(i);theta(i)*x]' * [1;x;theta(i);theta(i)*x]
%% Since there is an inequality constraint too, it will generate 2*N blocks
Aall            = {};
A1all           = {};
ball            = [];
Call            = {};
A0append        = {};
for blkidx = 1:N
    basis       = [1;x;theta(blkidx);theta(blkidx)*x];
    n           = length(basis);
    basis_x     = monomials(theta(blkidx),1:2);
    basis_theta = monomials(x,0:2);
    basis_g     = [1;theta(blkidx)];

    out         = gen_chordal_subblk_pcr(...
                    basis,basis_x,h_x,basis_theta,h_theta(blkidx),g_x,basis_g,f_cost(blkidx));
    Acell       = out.A;
    A1          = out.A1;
    b           = out.b;
    C           = out.C;

    % add constraint that the top-left [1;x]*[1;x]' block is the same as
    % the 0-th block
    A0blk = {};
    for i = 1:n0
        for j = i:n0
            if i == j
                A0i  = sparse(i,j,-1,n0,n0);
                Ai   = sparse(i,j,1,n,n);
            else
                A0i  = sparse([i,j],[j,i],[-0.5,-0.5],n0,n0);
                Ai   = sparse([i,j],[j,i],[0.5,0.5],n,n);
            end
            A0blk    = [A0blk;{A0i}];
            Acell    = [Acell;{Ai}];
        end
    end

    ball             = [ball;b;sparse(n0delta,1)];
    A0append{end+1}  = A0blk;
    Aall{end+1}      = Acell;
    Call             = [Call;C];
    A1all{end+1}     = A1;
end

%% Convert to standard SDPT3 format
b           = [b0;ball];
blk         = cell(2*N+1,2);
blk{1,1}    = 's'; blk{1,2} = n0;
n1          = out.blk{2,2};
n1delta     = triangle_number(n1);
for i = 1:N
    blk{2*i,1} = 's';
    blk{2*i,2} = n;
    blk{2*i+1,1} = 's';
    blk{2*i+1,2} = n1;
end

A0t     = sparsesvec(blk(1,:),A0);
for i = 1:N
    A0t = [A0t,...
            sparse(n0delta,out.m),...
            sparsesvec(blk(1,:),A0append{i})];
end
ndelta  = triangle_number(n);
At    = {A0t};
for i = 1:N
    Ait = [sparse(ndelta,length(b0)),...
            sparse(ndelta,(i-1)*length(Aall{i})),...
            sparsesvec(blk(2*i,:),Aall{i}),...
            sparse(ndelta,(N-i)*length(Aall{i}))];

    A1it = [sparse(n1delta,length(b0)),...
            sparse(n1delta,(i-1)*length(Aall{i})),...
            sparse(n1delta,out.m_mom+out.m_loc),...
            sparsesvec(blk(2*i+1,:),A1all{i}),...
            sparse(n1delta,n0delta),...
            sparse(n1delta,(N-i)*length(Aall{i}))];
    At  = [At;{Ait};{A1it}];
end

SDP.blk = blk;
SDP.At  = At;
SDP.n   = n;
SDP.m   = length(b);
SDP.C   = [{sparse(n0,n0)};Call];
SDP.b   = b;

tf    = toc(t0);
fprintf('\nDone in %g seconds.\n',tf);
fprintf('===================================================================\n')

%% Convert to sedumi
fprintf('Convert to sedumi')
t0 = tic;

sK.s  = [n0];
for i = 1:N
    sK.s        = [sK.s,n,n1];
end

A0t     = sparsevec(blk(1,:),A0);
n0sq    = n0^2;
for i = 1:N
    A0t = [A0t,...
            sparse(n0sq,out.m),...
            sparsevec(blk(1,:),A0append{i})];
end

nsq     = n^2;
n1sq    = n1^2;
At      = {A0t};
for i = 1:N
    Ait = [sparse(nsq,length(b0)),...
            sparse(nsq,(i-1)*length(Aall{i})),...
            sparsevec(blk(2*i,:),Aall{i}),...
            sparse(nsq,(N-i)*length(Aall{i}))];

    A1it = [sparse(n1sq,length(b0)),...
            sparse(n1sq,(i-1)*length(Aall{i})),...
            sparse(n1sq,out.m_mom+out.m_loc),...
            sparsevec(blk(2*i+1,:),A1all{i}),...
            sparse(n1sq,n0delta),...
            sparse(n1sq,(N-i)*length(Aall{i}))];
    At  = [At;{Ait};{A1it}];
end

sdata.K     = sK;
sdata.At    = cat(1,At{:});
sdata.b     = b;

sc          = [];
for i = 1:length(SDP.C)
    sc      = [sc;sparsevec(blk(i,:),SDP.C(i))];
end
sdata.c     = sc;

SDP.sedumi   = sdata;

tf    = toc(t0);
fprintf('\nDone in %g seconds.\n',tf);
fprintf('===================================================================\n')

end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function z_hat = estimate_states(uu, time, dt, p_drone, p_physics)
% ESTIMATE_STATES: Estimate the Drone states using gyros, accels, pressure
% sensors, and GPS.
% Inputs:
%   uu: input measurements
%   p_drone: drone parameters
%
% Outputs:
%   pnhat    - estimated North position,
%   pehat    - estimated East position,
%   hhat     - estimated altitude,
%   Vahat    - estimated airspeed,
%   alphahat - estimated angle of attack
%   betahat  - estimated sideslip angle
%   phihat   - estimated roll angle,
%   thetahat - estimated pitch angel,
%   chihat   - estimated course,
%   phat     - estimated roll rate,
%   qhat     - estimated pitch rate,
%   rhat     - estimated yaw rate,
%   Vghat    - estimated ground speed,
%   wnhat    - estimate of North wind,
%   wehat    - estimate of East wind
%   psihat   - estimate of heading angle
%

% Rename inputs
y_gyro_x      = uu(1);
y_gyro_y      = uu(2);
y_gyro_z      = uu(3);
y_accel_x     = uu(4);
y_accel_y     = uu(5);
y_accel_z     = uu(6);
y_static_pres = uu(7);
y_diff_pres   = uu(8);

y_gps_n       = uu(9);
y_gps_e       = uu(10);
y_gps_h       = uu(11);
y_gps_Vg      = uu(12);
y_gps_course  = uu(13);

% Persistent variables
persistent  kf1;
%persistent  kf2;
persistent x;
persistent kf_P;

% Initilize estimated states at t=0
if time == 0
    kf1.x_prev = [p_drone.phi_T, p_drone.theta_T]';
    kf1.P_prev = zeros(2,2);
    Vg_T = p_drone.Va_T;
    chi_T = p_drone.psi_T;
    %kf2.x_prev = [P.pn_T; P.pe_T; Vg_T; chi_T; P.wn; P.we; P.psi_T];
    %kf2.P_prev = zeros(7,7);
    phi_hat     = p_drone.phi_T;
    theta_hat   = p_drone.theta_T;
    x = [0; 0];
    kf_P = [0, 0;
        0  , 0];
end

% Not estimating these states
alphahat = 0;
betahat  = 0;
bxhat    = 0;
byhat    = 0;
bzhat    = 0;

pnhat=0;
pehat=0;
hhat=0;
Vahat=0;
alphahat=0;
betahat=0;
phihat=0;
thetahat=0;
chihat=0;
phat=0;
qhat=0;
rhat=0;
Vghat=0;
wnhat=0;
wehat=0;
psihat=0;
bxhat=0;
byhat=0;
bzhat=0;

%% Estimate states by low pass filtering

phat     = lp_filter( y_gyro_x, 0, p_drone.a_gyro, dt);
qhat     = lp_filter( y_gyro_y, 0, p_drone.a_gyro, dt);
rhat     = lp_filter( y_gyro_z, 0, p_drone.a_gyro, dt);

hhat     = (p_drone.pres0 - lp_filter( y_static_pres, 0, p_drone.a_static_pres, dt))/...
    (p_physics.rho*p_physics.gravity);
Vahat    = sqrt(2/p_physics.rho*lp_filter(y_diff_pres, 0, p_drone.a_diff_pres, dt));

%thetahat = asin(lp_filter(y_accel_x, 0, P.a_accel, P.dt)/p_physics.gravity);
%phihat   = atan(lp_filter(y_accel_y, 0, P.a_accel, P.dt)/lp_filter(y_accel_z, 0, P.a_accel, P.dt));

pnhat    = lp_filter(y_gps_n, 0, p_drone.a_gps_pos, dt);
pehat    = lp_filter(y_gps_e, 0, p_drone.a_gps_pos, dt);
chihat   = lp_filter(y_gps_course, 0, p_drone.a_gps_speed, dt);
Vghat    = lp_filter(y_gps_Vg, 0, p_drone.a_gps_speed, dt);

%% Estimated states by kalman filter

% Reconstruct theta and phi measured
% theta_m  = asin((y_accel_x)/p_physics.gravity);
% phi_m    = atan(y_accel_y/y_accel_z);
% last_measure1 = [ t       phi_m;           % not really t, but time of sensor measurement
%                   t       theta_m];        % not really t, ...
%
% input_z_hat = [phihat; thetahat; phat; qhat; rhat; Vahat];
% [F1, G1, H1, u1, Q1, R1] = ekf_phi_theta( input_z_hat, P);
% [kf1] = ext_kalman_filter(F1, G1, H1, u1, Q1, R1, t, last_measure1, kf1, P);

% Extract x to output variables
% phihat = kf1.x(1);
% thetahat = kf1.x(2);
x_prev = x;
phi_hat = x(1);
theta_hat = x(2);
kf_P_prev = kf_P;

z = [y_accel_x ; y_accel_y; y_accel_z];

% Compute f, F, h and H

cp = cos(phi_hat);
sp = sin(phi_hat);
ct = cos(theta_hat);
st = sin(theta_hat);
tt = tan(theta_hat);

p = y_gyro_x;
q = y_gyro_y;
r = y_gyro_z;
Va = Vahat;
g = p_physics.gravity;

f = [  p + q*sp*tt + r*cp*tt; ...
    q*cp - r*sp];

F = [  q*cp*tt - r*sp*tt,   (q*sp-r*cp)/(ct^2); ...
    - q*sp - r*cp,          0                ];

h = [ q*Va*st + g*st; ...
    r*Va*ct - p*Va*st - g*ct*sp; ...
    -q*Va*ct - g*ct*cp];

H = [  0,                 q*Va*ct + p_physics.gravity*ct;...
    - p_physics.gravity*cp*ct,   - r*Va*st - p*Va*ct + p_physics.gravity*sp*st;...
    p_physics.gravity*sp*ct,   (q*Va+p_physics.gravity*cp)*st               ];

%% Compute prediction (a priori)

x = x_prev+dt*f;    % a priori state estimate
kf_P = kf_P_prev + (F*kf_P_prev + kf_P_prev*F' + p_drone.Q);    % a priori estimate covariance

%% Compute update (a posteriori)

if(norm(z-p_physics.gravity) < 0.5)
    z_tilde = z - H*x;              % innovation (a priori residual)
    S = p_drone.R + H*kf_P*(H');    % innovation covariance
    K = kf_P*(H')/S;                % optimal kalman gain
    x = x + K*z_tilde;              % a posteriori state estimate
    kf_P = (eye(2) - K*H)*kf_P;     % a posteriori estimate covariance
end

z_res = z - H*x;                    % a posteriori residual

phihat      = x(1);
thetahat    = x(2);

%% Output data

z_hat = [
    pnhat;...
    pehat;...
    hhat;...
    0
    0
    0
    Vahat;...
    alphahat;...
    betahat;...
    phihat;...
    thetahat;...
    chihat;...
    phat;...
    qhat;...
    rhat;...
    Vghat;...
    wnhat;...
    wehat;...
    psihat;...
    bxhat;...
    byhat;...
    bzhat
    ];

end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function corr_coef = RegreesionTID2013(objectiveValues,mos)
%this script is used to calculate the pearson linear correlation
%coefficient and root mean sqaured error after regression

%get the objective scores computed by the IQA metric and the subjective
%scores provided by the dataset
% matData = load('VSIOnTID2013.mat');
% VSIOnTID2013 = matData.VSIOnTID2013;
% objectiveValues = VSIOnTID2013(:,1);
% mos = VSIOnTID2013(:,2);

%plot objective-subjective score pairs
% p = plot(objectiveValues,mos,'+');
% set(p,'Color','blue','LineWidth',1);

%initialize the parameters used by the nonlinear fitting function
beta(1) = 10;
beta(2) = 0;
beta(3) = mean(objectiveValues);
beta(4) = 0.1;
beta(5) = 0.1;
%fitting a curve using the data
[bayta ehat,J] = nlinfit(objectiveValues,mos,@logistic,beta);
%given an objective value, predict the correspoing mos (ypre) using the fitted curve
[ypre junk] = nlpredci(@logistic,objectiveValues,bayta,ehat,J);

RMSE = sqrt(sum((ypre - mos).^2) / length(mos));%root meas squared error
corr_coef = corr(mos, ypre, 'type','Pearson'); %pearson linear coefficient

%draw the fitted curve
% t = min(objectiveValues):0.01:max(objectiveValues);
% [ypre junk] = nlpredci(@logistic,t,bayta,ehat,J);
% hold on;
% p = plot(t,ypre);
% set(p,'Color','black','LineWidth',2);
% legend('Images in TID2013','Curve fitted with logistic function', 'Location','NorthWest');
% xlabel('Objective score by VSI');
% ylabel('MOS');

author: HuiZeng


text:
% EX_COLLOCATION_LAPLACE_RING: solve the Laplace problem with a NURBS discretization by isogeometric collocation

% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_ring.txt';

% Type of boundary conditions for each side of the domain
problem_data.nmnn_sides   = [1 2];
problem_data.drchlt_sides = [3 4];

% Physical parameters
problem_data.c_diff  = @(x, y) ones(size(x)); % Only used in Galerkin, not in collocation

% Source and boundary terms
problem_data.f = @(x, y) exp(x).*((x.^2 + y.^2 - 1).*sin(x.*y) - 2*y.*cos(x.*y));
problem_data.g = @test_ring_mixed_bc_g_nmnn;
problem_data.h = @(x, y, ind) exp(x).*sin(x.*y);

% Exact solution (optional)
problem_data.uex     = @(x, y) exp(x) .* sin (x.*y);
problem_data.graduex = @(x, y) cat (1, ...
               reshape ( exp(x).*(sin(x.*y) + y.*cos(x.*y)), [1, size(x)]), ...
               reshape (exp(x).*x.*cos(x.*y), [1, size(x)]));

% Discretization parameters (p and h)
method_data.degree     = [4 4]; % Degree of the splines in each direction
method_data.regularity = [3 3]; % Regularity of the splines, should be at least C^1.
method_data.nsub       = [9 9]; % Divide each subinterval of the original knot span in nsub subintervals
method_data.pts_case   = 1;     % Collocation points. 1: Greville abscissae, 2: clustered superconvergent points

% Solve with collocation method
[geometry, msh_coll, space_coll, u_coll] = solve_laplace_collocation (problem_data, method_data);

% Solve with Galerkin, for comparison
method_data.nquad = method_data.degree + 1;
[~, msh_gal, space_gal, u_gal] = solve_laplace_iso (problem_data, method_data);

% Plot of solution
vtk_pts = {linspace(0, 1, 20), linspace(0, 1, 20)};
figure;
[eu, F] = sp_eval (u_coll, space_coll, geometry, vtk_pts);
[X, Y]  = deal (squeeze(F(1,:,:)), squeeze(F(2,:,:)));
subplot (1,3,1)
surf (X, Y, eu)
title ('Collocation solution'), axis tight
[eu, F] = sp_eval (u_gal, space_gal, geometry, vtk_pts);
[X, Y]  = deal (squeeze(F(1,:,:)), squeeze(F(2,:,:)));
subplot (1,3,2)
surf (X, Y, eu)
title ('Galerkin solution'), axis tight
subplot (1,3,3)
surf (X, Y, problem_data.uex (X,Y))
title ('Exact solution'), axis tight

% Compute errors of collocation and Galerkin. Since it involves quadrature,
%  it is computed using the mesh and space objects from the Galerkin method.
disp ('Error in H1 and L2 norms, for isogeometric collocation')
[error_h1_coll, error_l2_coll] = sp_h1_error (space_gal, msh_gal, u_coll, problem_data.uex, problem_data.graduex);
disp([error_l2_coll error_h1_coll])

disp ('Error in H1 and L2 norms, for isogeometric Galerkin')
[error_h1_gal, error_l2_gal] = sp_h1_error (space_gal, msh_gal, u_gal, problem_data.uex, problem_data.graduex);
disp ([error_l2_gal error_h1_gal])

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function   s1 = sigma_8(k, k1, psm1, h)
% sigma_8.m calculates  the integrant for the expected rms mass overdensity in a sphere of radius 8 Mpc/h.
% input : k, wavenumber in units h/Mpc and the vectors k1 resp psm1, the mass power spectrum as a
% function of k1
% output : the integrant to obtain sigma_8, the rms mass overdensity sampled with
% a sphere of radius 8/h Mpc.

% D Vangheluwe 15 july 2005

% radius of the sphere in units Mpc/h
r0 = 8/h;
x = k * r0;
% calculate the point in the mass power spectrum at k :
psm = spline(k1, psm1, k);

% calculate the integrant:
s1 = (3*(sin(x) - x .* cos(x)) ./ (x .^3)) .^2 .* psm .* k .^2;

return

author: Sable

repo: mcbench-benchmarks


text:
%kkthreshbelow 'Threshold Data Values that are Below Specified Cutoff'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kthreshbelow.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input file'
% OutputFile: o ' Output ', required: 'output file'
%
% Example: o = kkthreshbelow(i, {'i','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kthreshbelow - Threshold Data Values that are Below Specified Cutoff
%
%  DESCRIPTION
% The "Threshold Below" operator generates a two-level output by
% thresholding the "Input" data object (i).  Data values less than
% the "Cutoff Value" (lc) are set to the "LESS THAN" output value
% (fval), and values greater than or equal to the "Cutoff Value" are set
% to the "GREATER THAN or EQUAL" Output value (tval).
%
% The following histogram defines the results of the thresholding operation.
%
%
% \f(CW
%                                 |<-F->|<----T---->|
%                                 |     | :     .   |
%                                 |    .:.:.  .:::  |
%                                 |:. .::::::.::::: |
%                                 |::::::::::::::::.|
%                                 ------|------------
%                                       Cutoff
%
%         F:  if value <= Cutoff,  Output = fval
%         T:  if value >  Cutoff,  Output = tval
% "
%
%
%  "Data Type" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "Map Data" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_1input
%
%  "Validity Mask" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "Explicit Location and Time Data" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "Failure Modes"
% .cI $DATAMANIP/repos/shared/man/sections/fail_1input
%
% Executing "Threshold Above" runs the program \fIkclip\fP with the -thresh
% option set to TRUE (-thresh 1).
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::kthreshabove, DATAMANIP::kthreshin, DATAMANIP::kthreshout,
% DATAMANIP::kmegathresh, DATAMANIP::kclip
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkthreshbelow(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkthreshbelow(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output'};
maxval={0,0};
minval={0,0};
istoggle=[0,0];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kclip"  -thresh'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
addpath('../../matlab/');

for n=[1000 5000 10000 20000]
C = int32( round(rand(n)*1e6) );
tic; [rho,varrho,u,v] = hungarianLSAP(C); toc
end

disp(['min cost = ',num2str(sum(u)+sum(v))]);

author: optimaltransport


text:
function W = weightingscheme(cases,wfct,varargin)

%
% INTERNAL FUNCTION
%

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

switch wfct

    case {'WHuber','whuber'}
        b = varargin{1};
        cases = abs(cases);
        temp = 1;
        W = temp.*(cases<b) + temp.*(cases>=b).*(b./cases);

    case {'WLogistic','wlogistic'}
        W = tanh(cases)./cases;

    case {'WHampel','whampel'}
       cases = abs(cases);
       % defaults for c1 and c2
       c1=2.5; c2=3;
       %[c1,c2] = adaptweight(cases);
       dc = c2-c1;
       temp = 1;
       W = temp.*(cases<=c1) + ...
           temp.*(cases<=c2 & cases>c1).*((c2-cases)./dc) + ...
           temp.*(cases>c2).*10e-8;

    case {'WMyriad','wmyriad'}
        %K = 0.5*iqr(cases);
        K = varargin{1};
        W = K^2./(K^2+cases.^2);
end

author: peterhcharlton

repo: RRest


text:
% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function v = double(q)
%Quaternion.double Convert a quaternion object to a 4-element vector
%
% V = Q.double() is a 4-vector comprising the quaternion
% elements [s vx vy vz].

    v = [q.s q.v];
endfunction

author: petercorke

repo: robotics-toolbox-matlab


text:
function vl_test_print(varargin)

addpath(fullfile(vl_rootnn(), 'examples', 'mnist'));

net = cnn_mnist_init('networkType', 'dagnn');
net.print(varargin{:});

end


text:
classdef LaunchVehicleViewProfileRefFrameData < matlab.mixin.SetGet
    %LaunchVehicleViewProfileUserRefFrameData Summary of this class goes here
    %   Detailed explanation goes here

    properties
        %these are for vehicle trajectory data and not the point itself!
        timesArr(1,:) cell = {};
        xInterps(1,:) cell = {};
        yInterps(1,:) cell = {};
        zInterps(1,:) cell = {};
        vxInterps(1,:) cell = {};
        vyInterps(1,:) cell = {};
        vzInterps(1,:) cell = {};

        markerPlot = matlab.graphics.GraphicsPlaceholder.empty(1,0)
        refFrame AbstractGeometricRefFrame
        viewFrame AbstractReferenceFrame
    end

    methods
        function obj = LaunchVehicleViewProfileRefFrameData(refFrame, viewFrame)
            obj.refFrame = refFrame;
            obj.viewFrame = viewFrame;
        end

        function addData(obj, times, rVects, vVects)
            obj.timesArr(end+1) = {times};

            if(length(times) >= 3)
                method = 'pchip';
            else
                method = 'linear';
            end

            obj.xInterps{end+1} = griddedInterpolant(times, rVects(:,1), method, 'linear');
            obj.yInterps{end+1} = griddedInterpolant(times, rVects(:,2), method, 'linear');
            obj.zInterps{end+1} = griddedInterpolant(times, rVects(:,3), method, 'linear');
            obj.vxInterps{end+1} = griddedInterpolant(times, vVects(:,1), method, 'linear');
            obj.vyInterps{end+1} = griddedInterpolant(times, vVects(:,2), method, 'linear');
            obj.vzInterps{end+1} = griddedInterpolant(times, vVects(:,3), method, 'linear');
        end

        function plotRefFrameAtTime(obj, time, hAx)
            for(i=1:length(obj.timesArr))
                times = obj.timesArr{i};

                if(time >= min(floor(times)) && time <= max(ceil(times)))
                    xInterp = obj.xInterps{i};
                    x = xInterp(time);

                    yInterp = obj.yInterps{i};
                    y = yInterp(time);

                    zInterp = obj.zInterps{i};
                    z = zInterp(time);

                    vxInterp = obj.vxInterps{i};
                    vx = vxInterp(time);

                    vyInterp = obj.vyInterps{i};
                    vy = vyInterp(time);

                    vzInterp = obj.vzInterps{i};
                    vz = vzInterp(time);

                    vehElemSet = CartesianElementSet(time, [x;y;z], [vx;vy;vz], obj.viewFrame);
                    [posOffsetOrigin, ~, ~, rotMatToInertial] = obj.refFrame.getRefFrameAtTime(time, vehElemSet, obj.viewFrame);

                    scaleFactor = obj.refFrame.scaleFactor;
                    xAxis = [posOffsetOrigin(:), posOffsetOrigin(:) + rotMatToInertial(:,1)*scaleFactor];
                    yAxis = [posOffsetOrigin(:), posOffsetOrigin(:) + rotMatToInertial(:,2)*scaleFactor];
                    zAxis = [posOffsetOrigin(:), posOffsetOrigin(:) + rotMatToInertial(:,3)*scaleFactor];

                    if(isempty(obj.markerPlot))
                        hold(hAx,'on');
                        obj.markerPlot(end+1) = plot3(hAx, xAxis(1,:), xAxis(2,:), xAxis(3,:), 'Color',obj.refFrame.xAxisColor.color, 'LineWidth',obj.refFrame.xAxisLineWidth, 'LineStyle',obj.refFrame.xAxisLineSpec.linespec);
                        obj.markerPlot(end+1) = plot3(hAx, yAxis(1,:), yAxis(2,:), yAxis(3,:), 'Color',obj.refFrame.yAxisColor.color, 'LineWidth',obj.refFrame.yAxisLineWidth, 'LineStyle',obj.refFrame.yAxisLineSpec.linespec);
                        obj.markerPlot(end+1) = plot3(hAx, zAxis(1,:), zAxis(2,:), zAxis(3,:), 'Color',obj.refFrame.zAxisColor.color, 'LineWidth',obj.refFrame.zAxisLineWidth, 'LineStyle',obj.refFrame.zAxisLineSpec.linespec);
                        hold(hAx,'off');
                    else
                        obj.markerPlot(1).XData = xAxis(1,:);
                        obj.markerPlot(1).YData = xAxis(2,:);
                        obj.markerPlot(1).ZData = xAxis(3,:);

                        obj.markerPlot(2).XData = yAxis(1,:);
                        obj.markerPlot(2).YData = yAxis(2,:);
                        obj.markerPlot(2).ZData = yAxis(3,:);

                        obj.markerPlot(3).XData = zAxis(1,:);
                        obj.markerPlot(3).YData = zAxis(2,:);
                        obj.markerPlot(3).ZData = zAxis(3,:);
                    end
                end
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function y = norms( x, p, dim )

%NORMS   Internal cvx version.

%
% Check second argument
%

error( nargchk( 1, 3, nargin ) );
if nargin < 2 || isempty( p ),
    p = 2;
elseif ~isnumeric( p ) || numel( p ) ~= 1 || ~isreal( p ),
    error( 'Second argument must be a real number.' );
elseif p < 1 || isnan( p ),
    error( 'Second argument must be between 1 and +Inf, inclusive.' );
end

%
% Check third argument
%

sx = size( x );
if nargin < 3 || isempty( dim ),
    dim = cvx_default_dimension( sx );
elseif ~cvx_check_dimension( dim, false ),
    error( 'Third argument must be a valid dimension.' );
end
if isempty( x ) || length( sx ) < dim || sx( dim ) == 1,
    p = 1;
end

%
% Type check
%

persistent remap1 remap2
if isempty( remap2 ),
    remap1 = cvx_remap( 'constant', 'log-convex' );
    remap2 = cvx_remap( 'affine', 'log-convex' );
end
xc = reshape( cvx_classify( x ), sx );
if ~all( remap2( xc( : ) ) ),
    error( 'Disciplined convex programming error:\n   Invalid computation: norms( {%s}, ... )', cvx_class( x, true, true ) );
end

%
% Compute norms
%

switch p,
    case 1,
        y = sum( abs(x), dim );
    case Inf,
        y = max( abs(x), [], dim );
    otherwise,
        nd = length( sx );
        nx = sx( dim );
        sy = sx;
        sy( dim ) = 1;
        tt = all( remap1( xc ), dim );
        if all( tt( : ) ),
            y = sum( abs(x) .^ p, dim ) .^ (1/p);
        elseif any( tt( : ) ),
            if dim ~= 1,
                perm = [ dim, 1:dim-1, dim+1:length(sx) ];
                x  = permute( x,  perm );
                tt = permute( tt, perm );
                sy = sy( perm );
            else
                perm = [];
            end
            nv = prod( sy );
            x  = reshape( x, nx, nv );
            y  = cvx( [ 1, nv ], [] );
            xt = cvx_subsref( x, ':', tt );
            y  = cvx_subsasgn( y, tt, sum( (abs(xt)).^p, 1 ) .^ (1/p) );
            tt = ~tt;
            xt = cvx_subsref( x, ':', tt );
            y  = cvx_subsasgn( y, tt, norms( xt, p, 1 ) );
            y  = reshape( y, sy );
            if ~isempty( perm ),
                y = ipermute( y, perm );
            end
        elseif p == 2,
            cvx_begin
                epigraph variable y( sy )
                { cvx_accept_convex(x), y } == lorentz( sx, dim, ~isreal( x ) ); %#ok
            cvx_end
        else
            sw = sx;
            sw(nd+1) = 2;
            cvx_begin
                variable z( sx )
                epigraph variable y( sy )
                if isreal(x), cmode = 'abs'; else cmode = 'cabs'; end
                { cat( nd+1, z, cvx_expand_dim( y, dim, nx ) ), cvx_accept_convex(x) } ...
                    == geo_mean_cone( sw, nd+1, [1/p,1-1/p], cmode ); %#ok
                sum( z, dim ) == y;
            cvx_end
        end
end

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

repo: REST


text:
function [u,X,T,uf,t]=forcmove(a,v,tmax,nt)

% [u,X,T,uf,t]=forcmove(a,v,tmax,nt)
% See Article 9.4
% This function computes the dynamic response
% of a taut string subjected to a force moving
% along the string with constant speed. The
% string is fixed at x=0 and and x=+infinity.
% The system is initially at rest when the
% force is imposed at the left end and moves
% to the right at constant speed. If the force
% speed is larger than the wave propagation
% speed, then no disturbance occurs ahead of
% the force. If the force moves slower than the
% wave propagation speed, then the deflection
% propagates ahead of the force at the speed
% of wave propagation.
%
% v     - speed of the moving load
% a     - speed of wave propagation in the
%         string
% tmax  - maximum time for which the
%         solution is computed
% u     - matrix of deflection values where
%         time and position vary row-wise and
%         column-wise, respectively
% T,X   - matrices of time and position values
%         corresponding to the deflection
%         matrix U
% uf    - deflection values where the force acts
% t     - vector of times (same as colunms of T)

if nargin==0, a=.8; v=1; tmax=10; nt=15; end

% Obtain solution values and plot results
[u,X,T,uf,t]=ustring(a,v,tmax,nt);
if a>v, xf=X(:,2); uf=u(:,2); xw=X(:,3);
else, xf=X(:,3); uf=u(:,3); end
close, subplot(211)
waterfall(X,T,-u), xlabel('x axis')
ylabel('time'), zlabel('deflection')
titl=['INVERTED MOTION SURFACE FOR  A = ',...
      num2str(a),'  AND  V = ',num2str(v)];
title(titl), grid on, hold on
plot3(xf,t,-uf,'.k',xf,t,-uf,'k')
colormap([0 0 0]), view([-10,30]), shg
umin=min(u(:)); umax=max(u(:)); xmax=X(1,4);
range=[0,xmax,2*umin,2*umax]; hold on
Titl='T = %4.2f'; subplot(212) , axis off

% Use a dense set of points for animation
nt=80; [uu,XX,TT,uuf,tt]=ustring(a,v,tmax,nt);
umax=max(abs(uu(:))); uu=uu/umax; uuf=uuf/umax;
XX=XX/xmax; range=[0,1,-1,0.5]; h=.4;
arx=h*[0,.02,-.02,0,0]; ary=h*[0,.25,.25,0,1];
for j=1:nt
   uj=uu(j,:); xj=XX(j,:);
   xfj=v/xmax*tt(j); ufj=uuf(j);
   plot(xj,uj,'k',xfj+arx,ufj+ary,'-k')
   axis off, time=(sprintf(Titl,tt(j)));
   text(.45,.25,time), axis(range), drawnow
   pause(.05), figure(gcf), if j<nt, cla, end
end
hold off;

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

function [u,X,T,uf,t]=ustring(a,v,tmax,nt)
% [u,X,T,uf,t]=ustring(a,v,tmax,nt)
% This function computes the deflection u(x,t)
% of a semi-infinite string subjected to a
% moving force. The equation for the normalized
% deflection
% u(x,t)=1/a/(a^2-v^2)*((v-a-v*abs(x-a*t)...
%                           +a*abs(x-v*t));
% a    - speed of wave propagation in the string
% v    - speed of the force moving to the right
% tmax - maximum time for computing the solution
% nt   - number of time increments computed
% uu   - array of displacement values normalized
%        by dividing by a factor equal to the force
%        magnitude over twice the density per unit
%        length. Position varies column-wise and
%        time varies row-wise in the array.
% X,T  - position and time arrays for the solution
% uf   - deflection vector under the force
% t    - time vector for the solution (same as the
%        columns of T)
%
t=linspace(0,tmax,nt)'; xmax=1.05*tmax*max(a,v);
u=zeros(nt,4); nx=4; X=zeros(nt,nx); X(:,nx)=xmax;
c=1/a/(a^2-v^2); xw=a*t; xf=v*t; T=repmat(t,1,4);
uw=c*xw*(v-a+abs(v-a)); uf=c*xf*(v-a-abs(v-a));
if a>v
   X(:,2)=xf; X(:,3)=xw; u(:,2)=uf;
else
   X(:,2)=xw; X(:,3)=xf; u(:,2)=uw;
end


text:
%  Figure 3.1      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%script to generate Fig. 3.1
%% fig3_01.m Example 3.4
% Frequency response
clf;
k=1;
num=1;                      % form numerator
den=[1 k];                  % form denominator
sys=tf(num,den);            % form system
% frequency range
w=logspace(-2,2);
[mag,phase]=bode(sys,w);    % compute frequency response
%plot frequency response
subplot(2,1,1); loglog(w,mag(:));xlabel('\omega (rad/sec)');
ylabel('M');title('Figure 3.1: Magnitude, phase');grid;
subplot(2,1,2); semilogx(w,phase(:));xlabel('\omega (rad/sec)');
ylabel('\phi (deg)');grid;