text:
function [success] = downsampleHdf5Movie(inputFilePath, varargin)
	% [success] = downsampleHdf5Movie(inputFilePath, varargin)
	%
	% Downsamples a movie (HDF5, SLD, etc.) in inputFilePath piece by piece and appends it to already created hdf5 file. Useful for large files that otherwise would not fit into RAM.
	%
	% Biafra Ahanonu
	% started: 2013.12.19
	%
	% Part of append code based on work by Dinesh Iyer (http://www.mathworks.com/matlabcentral/newsreader/author/130530)
	% note: checked between this implementation and imageJ's scale, they are nearly identical (subtracted histogram mean is 1 vs each image mean of ~1860, so assuming precision error).
	%
	% inputs
	%	inputFilePath - Str: path to movie to downsample into HDF5.
	% outputs
	% 	success - Binary: 1 = completed successfully, 0 = did not complete successfully.
	% options
	%	inputDatasetName = hierarchy where data is stored in HDF5 file

	% changelog
		% 2014.01.18 - improved method of obtaining the newFilename
		% 2014.06.16 - updated output notifications to user
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2022.07.10 [19:24:29] - Added Bio-Formats support.
		% 2022.07.15 [13:34:28] - Switch chunking to automatic by default, generally improved reading from data later than small chunks.
	% TODO
		% Use handles to reduce memory load when doing computations.

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% Str: name of the input hierarchy in the HDF5 files
	options.inputDatasetName = '/1';
	% Str: name of the output hierarchy in the HDF5 files
	options.outputDatasetName = '/1';
	% Vector: frames to use.
	options.frameList = 1;
	% Int: amount to downsample movie. Value of 1 = no downsample (e.g. just file conversion).
	options.downsampleFactor = 4;
	% Int: max size of a chunk in Mbytes. For HDF5 files.
	options.maxChunkSize = 20000;
	% Int: bytes to MB conversion. DO NOT CHANGE.
	options.bytesToMB = 1024^2;
	% Int: number of frames to use when loading file in chunks. For Bio-Formats files!
	options.numFramesSubset = 1000;
	% interval over which to show waitbar
	options.waitbarInterval = 1000;
	% Str: new filename to save to for the downsampled movie.
	options.newFilename = '';
	% Str: downsample to different folder. Leave blank to use same folder as input file.
	options.saveFolder = [];
	% second if want to do another downsample without loading another file
	options.saveFolderTwo = [];
	% Int: amount to downsample movie to 2nd folder. Value of 1 = no downsample (e.g. just file conversion).
	options.downsampleFactorTwo = 2;
	% Str: new filename to save to for the downsampled movie in 2nd folder.
	options.newFilenameTwo = '';
	% Int: Defines gzip compression level (0-9). 0 = no compression, 9 = most compression.
	options.deflateLevel = 1;
	% Int: chunk size in [x y z] of the dataset, leave empty for auto chunking
	options.dataDimsChunkCopy = []; % [128 128 1]
	% Str: movie type: hdf5, bioformats. Leave blank for automatic determination.
	options.movieType = '';
	% Int: Bio-Formats series number to load.
	options.bfSeriesNo = 1;
	% Int: Bio-Formats channel number to load.
	options.bfChannelNo = 1;
	% Int: Bio-Formats z dimension to load.
	options.bfZdimNo = 1;
	% Str: Bio-Formats series name. Function will search via regexp for this series within a Bio-Formats file.
	options.bfSeriesName = '';
	% Binary: 1 = whether to display info on command line.
	options.displayInfo = 1;
	% Binary: 1 = waitbar/progress bar is shown, 0 = no progress shown.
	options.waitbarOn = 1;
	% get options
	options = getOptions(options,varargin);
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	%     eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	[pathstr,name,ext] = fileparts(inputFilePath);
	inputFileEmpty = 0;
	if isempty(options.newFilename)
		inputFileEmpty = 1;
		options.newFilename = [pathstr filesep 'concat_' name '.h5'];
	end
	if isempty(options.newFilenameTwo)
		options.newFilenameTwo = [pathstr filesep 'concat_' name '.h5'];
	end

	% Check that options.dataDimsChunkCopy

	% Check type of movie being read in.
	if isempty(options.movieType)
		[options.movieType] = ciapkg.io.getMovieFileType(inputFilePath);
	end

	startTime = tic;

	fileIdOpen = [];

	switch options.movieType
		case 'hdf5'
			% movie dimensions and subsets to analyze
			[subsets, dataDim] = getSubsetOfDataToAnalyze(inputFilePath, options, varargin);
		case 'bioformats'
			disp(['Opening connection to Bio-Formats file:' inputFilePath])
			startReadTime = tic;
			fileIdOpen = bfGetReader(inputFilePath);
			toc(startReadTime)

			% If user inputs a series name, find within the file and overwrite default series name.
			if ~isempty(options.bfSeriesName)
				[bfSeriesNo] = ciapkg.bf.getSeriesNoFromName(inputFilePath,options.bfSeriesName,'fileIdOpen',fileIdOpen);
				if isnan(bfSeriesNo)
					disp(['Series "' options.bfSeriesName '" not found in ' inputFilePath])
					disp('Returning...')
					return;
				else
					options.bfSeriesNo = bfSeriesNo;
				end
			end

			% Convert for 0-base indexing.
			bfSeriesNoTrue = options.bfSeriesNo-1;

			fileIdOpen.setSeries(bfSeriesNoTrue);
			omeMeta = fileIdOpen.getMetadataStore();

			dim = struct;
			dataDim.x = omeMeta.getPixelsSizeX(bfSeriesNoTrue).getValue(); % image width, pixels
			dataDim.y = omeMeta.getPixelsSizeY(bfSeriesNoTrue).getValue(); % image height, pixels
			dataDim.z = fileIdOpen.getSizeT;

			thisStr = char(omeMeta.getImageName(bfSeriesNoTrue));

			disp(['Loading series: ' thisStr])
			display(dataDim)

			% If filename is empty
			if inputFileEmpty==1
				options.newFilename = [pathstr filesep 'concat_' thisStr '.h5'];
			end

			disp('===')
			disp(['Series in file: ' inputFilePath])
			nSeries = fileIdOpen.getSeriesCount();
			for seriesNo = 1:nSeries
				thisStr = char(omeMeta.getImageName(seriesNo-1));
				disp([num2str(seriesNo) ': "' thisStr '"'])
			end
			disp('===')

			% Get the subsets of the movie to analyze
			subsetSize = options.numFramesSubset;
			movieLength = dataDim.z;
			numSubsets = ceil(movieLength/subsetSize)+1;
			subsets = round(linspace(1,movieLength,numSubsets));
		otherwise
	end

	subsetsDiff = diff(subsets);
	% to compensate for flooring of linspace
	subsetsDiff(end) = subsetsDiff(end)+1;
	display(['subsets: ' num2str(subsets)]);
	display(['subsets diffs: ' num2str(subsetsDiff)]);

	% Change output filename if user request a different folder than where original file is.
	if ~isempty(options.saveFolder)
		[pathstr,name,ext] = fileparts(options.newFilename);
		options.newFilename = [options.saveFolder filesep name '.h5'];
	end
	if ~isempty(options.saveFolderTwo)
		[pathstr,name,ext] = fileparts(options.newFilename);
		options.newFilenameTwo = [options.saveFolderTwo filesep name '.h5'];
	end

	% Pre-establish connection to file to save time.
	switch options.movieType
		case 'hdf5'
			% hdf5FileWorkerConstant.Value = H5F.open(inputMovie);
		case 'bioformats'
			if isempty(fileIdOpen)
				disp(['Opening connection to Bio-Formats file:' inputFilePath])
				startReadTime = tic;
				fileIdOpen = bfGetReader(inputFilePath);
				toc(startReadTime)
			end
		otherwise
	end
	display(repmat('+',1,21))
	disp(['Saving to: ' options.newFilename])
	% display(['saving to: ' options.newFilename])

	try
		nSubsets = (length(subsets)-1);
		for currentSubset = 1:nSubsets
			loopStartTime = tic;
			% get current subset location and size
			currentSubsetLocation = subsets(currentSubset);
			lengthSubset = subsetsDiff(currentSubset);
			% convert offset to C-style offset for low-level HDF5 functions
			offset = [0 0 currentSubsetLocation-1];
			block = [dataDim.x dataDim.y lengthSubset];
			display('---')
			% display(sprintf(['current location: ' num2str(round(currentSubsetLocation/dataDim.z*100)) '% | ' num2str(currentSubsetLocation) '/' num2str(dataDim.z) '\noffset: ' num2str(offset) '\nblock: ' num2str(block)]));
			fprintf('current location: %d%% | %d/%d \noffset: %s \nblock: %s\n',round(currentSubsetLocation/dataDim.z*100),currentSubsetLocation,dataDim.z,mat2str(offset),mat2str(block));

			switch options.movieType
				case 'hdf5'
					% load subset of HDF5 file into memory
					inputMovie = readHDF5Subset(inputFilePath,offset,block,'datasetName',options.inputDatasetName);
				case 'bioformats'
					thisFrameList = currentSubsetLocation:(currentSubsetLocation+lengthSubset-1);
					disp(['Frames ' num2str(thisFrameList(1)) ' to ' num2str(thisFrameList(end))])
					[inputMovie] = local_readBioformats(inputFilePath,thisFrameList,fileIdOpen,options);
				otherwise
			end

			% split into second movie if need be
			if ~isempty(options.saveFolderTwo)
				inputMovieTwo = inputMovie;
			end

			% Downsample section of the movie, keep in memory
			if options.downsampleFactor==1
			else
				downsampleMovieNested('downsampleDimension', 'space','downsampleFactor',options.downsampleFactor,'waitbarInterval',options.waitbarInterval);
			end
			% thisMovie = uint16(thisMovie);
			display(['subset class: ' class(inputMovie)])

			% For snapshots with a single frame
			if size(inputMovie,3)==1
				inputMovie(:,:,2) = inputMovie;
			end

			% save the movie
			if currentSubset==1
				createHdf5File(options.newFilename, options.outputDatasetName, inputMovie,'deflateLevel',options.deflateLevel,'dataDimsChunkCopy',options.dataDimsChunkCopy);
			else
				appendDataToHdf5(options.newFilename, options.outputDatasetName, inputMovie);
			end
			toc(loopStartTime);
			display(sprintf(['downsample dims: ' num2str(size(inputMovie)) '\n-------']));

			if ~isempty(options.saveFolderTwo)
				display('secondary downsample in progress...')
				inputMovie = inputMovieTwo;
				% downsample section of the movie, keep in memory
				downsampleMovieNested('downsampleDimension', 'space','downsampleFactor',options.downsampleFactorTwo,'waitbarInterval',options.waitbarInterval);
				% save the movie
				if currentSubset==1
					createHdf5File(options.newFilenameTwo, options.outputDatasetName, inputMovie,'deflateLevel',options.deflateLevel,'dataDimsChunkCopy',options.dataDimsChunkCopy);
				else
					appendDataToHdf5(options.newFilenameTwo, options.outputDatasetName, inputMovie);
				end
			end
			% clear inputMovie;
		end
		success = 1;
		options.inputDatasetName = options.outputDatasetName;
		[subsets dataDim] = getSubsetOfDataToAnalyze(options.newFilename, options, varargin);
		display(repmat('+',1,7))
		display(['final HDF5 dims: ' num2str(cell2mat(struct2cell(dataDim))')]);
		display(repmat('+',1,7))
	catch err
		display(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		display(repmat('@',1,7))
		success = 0;
	end
	toc(startTime);

	function downsampleMovieNested(varargin)
		% downsamples a movie in either space or time, uses floor to calculate downsampled dimensions.
		% biafra ahanonu
		% started 2013.11.09 [09:31:32]
		%
		% inputs
			% inputMovie: a NxMxP matrix
		% options
			% downsampleType
			% downsampleFactor - amount to downsample in time
		% changelog
			% 2013.12.19 added the spatial downsampling to the function.
		% TODO

		%========================
		% default options
		nestedoptions.downsampleDimension = 'time';
		nestedoptions.downsampleType = 'bilinear';
		nestedoptions.downsampleFactor = 4;
		% exact dimensions to downsample in Z
		nestedoptions.downsampleZ = [];
		nestedoptions.waitbarOn = 1;
		% number of frames in each movie to load, [] = all, 1:500 would be 1st to 500th frame.
		nestedoptions.frameList = [];
		% whether to convert movie to double on load, not recommended
		nestedoptions.convertToDouble = 0;
		% name of HDF5 dataset name to load
		nestedoptions.inputDatasetName = '/1';
		% interval over which to show waitbar
		nestedoptions.waitbarInterval = 1000;
		% get user options, else keeps the defaults
		nestedoptions = getOptions(nestedoptions,varargin);
		% unpack options into current workspace
		% fn=fieldnames(options);
		% for i=1:length(fn)
		%     eval([fn{i} '=options.' fn{i} ';']);
		% end
		%========================
		% load the movie within downsample function
		if strcmp(class(inputMovie),'char')|strcmp(class(inputMovie),'cell')
			inputMovie = loadMovieList(inputMovie,'convertToDouble',nestedoptions.convertToDouble,'frameList',nestedoptions.frameList,'inputDatasetName',nestedoptions.inputDatasetName);
		end

		switch nestedoptions.downsampleDimension
			case 'time'
				switch nestedoptions.downsampleType
					case 'bilinear'
						% we do a bit of trickery here: we can downsample the movie in time by downsampling the X*Z 'image' in the Z-plane then stacking these downsampled images in the Y-plane. Would work the same of did Y*Z and stacked in X-plane.
						downX = size(inputMovie,1);
						downY = size(inputMovie,2);
						if isempty(nestedoptions.downsampleZ)
							downZ = floor(size(inputMovie,3)/nestedoptions.downsampleFactor);
						else
							downZ = nestedoptions.downsampleZ;
						end
						% pre-allocate movie
						% inputMovieDownsampled = zeros([downX downY downZ]);
						% this is a normal for loop at the moment, if convert inputMovie to cell array, can force it to be parallel
						reverseStr = '';
						for frame=1:downY
						   downsampledFrame = imresize(squeeze(inputMovie(:,frame,:)),[downX downZ],'bilinear');
						   % to reduce memory footprint, place new frame in old movie and cut off the unneeded frames after
						   inputMovie(1:downX,frame,1:downZ) = downsampledFrame;
						   % inputMovie(:,frame,:) = downsampledFrame;
							if frame==1||mod(frame,nestedoptions.waitbarInterval)==0&nestedoptions.waitbarOn==1|frame==downY
								reverseStr = cmdWaitbar(frame,downY,reverseStr,'inputStr','temporally downsampling matrix');
							end
						end
						inputMovie = inputMovie(:,:,1:downZ);
						drawnow;
					otherwise
						return;
				end
			case 'space'
				switch nestedoptions.downsampleType
					case 'bilinear'
						% we do a bit of trickery here: we can downsample the movie in time by downsampling the X*Z 'image' in the Z-plane then stacking these downsampled images in the Y-plane. Would work the same of did Y*Z and stacked in X-plane.
						downX = floor(size(inputMovie,1)/nestedoptions.downsampleFactor);
						downY = floor(size(inputMovie,2)/nestedoptions.downsampleFactor);
						downZ = size(inputMovie,3);
						% pre-allocate movie
						% inputMovieDownsampled = zeros([downX downY downZ]);
						% this is a normal for loop at the moment, if convert inputMovie to cell array, can force it to be parallel
						reverseStr = '';
						for frame=1:downZ
						   downsampledFrame = imresize(squeeze(inputMovie(:,:,frame)),[downX downY],'bilinear');
						   % to reduce memory footprint, place new frame in old movie and cut off the unneeded frames after
						   inputMovie(1:downX,1:downY,frame) = downsampledFrame;
						   % inputMovieDownsampled(1:downX,1:downY,frame) = downsampledFrame;
							if frame==1||mod(frame,nestedoptions.waitbarInterval)==0&nestedoptions.waitbarOn==1|frame==downZ
								reverseStr = cmdWaitbar(frame,downZ,reverseStr,'inputStr','spatially downsampling matrix');
							end
						end
						inputMovie = inputMovie(1:downX,1:downY,:);
						drawnow;
					otherwise
						return;
				end
			otherwise
				display('incorrect dimension option, choose time or space');
		end
		% display(' ');
	end
end

function [outputMovie] = local_readBioformats(thisMoviePath,thisFrameList,bfreaderTmp,options)
	% Reads in Bio-Formats data in chunks, uses open ID to speed up read times.

	% Setup movie class
	numMovies = 1;

	% Account for 0-base indexing.
	bfSeriesNoTrue = options.bfSeriesNo-1;

	% bfreaderTmp = bfGetReader(thisMoviePath);
	bfreaderTmp.setSeries(bfSeriesNoTrue);

	omeMeta = bfreaderTmp.getMetadataStore();
	stackSizeX = omeMeta.getPixelsSizeX(bfSeriesNoTrue).getValue(); % image width, pixels
	stackSizeY = omeMeta.getPixelsSizeY(bfSeriesNoTrue).getValue(); % image height, pixels
	stackSizeZ = omeMeta.getPixelsSizeZ(bfSeriesNoTrue).getValue();
	nChannels = omeMeta.getChannelCount(bfSeriesNoTrue);

	% Get the number of time points (frames).
	nFramesHere = bfreaderTmp.getSizeT;

	xyDims = [stackSizeX stackSizeY];

	if isempty(thisFrameList)
		nFrames = nFramesHere;
		framesToGrab = 1:nFrames;
	else
		nFrames = length(thisFrameList);
		framesToGrab = thisFrameList;
	end
	vidHeight = xyDims(2);
	vidWidth = xyDims(1);

	iPlane = bfreaderTmp.getIndex(0, options.bfChannelNo-1, 0)+1;
	tmpFrame = bfGetPlane(bfreaderTmp, iPlane);
	imgClass = class(tmpFrame);

	% Preallocate movie structure.
	if numMovies==1
		outputMovie = zeros(vidHeight, vidWidth, nFrames, imgClass);
	else
		tmpMovie = zeros(vidHeight, vidWidth, nFrames, imgClass);
	end

	% Read one frame at a time.
	reverseStr = '';
	iframe = 1;
	nFrames = length(framesToGrab);

	zPlane = 1;
	iZ = 1;
	dispEvery = round(nFrames/20);
	% Adjust for zero-based indexing
	% framesToGrab = framesToGrab-1;

	% figure;
	reverseStr = '';
	for t = framesToGrab
		% Assume 1-based adjust to 0-based indexing
		iPlane = bfreaderTmp.getIndex(zPlane-1, options.bfChannelNo-1, t-1)+1;
		tmpFrame = bfGetPlane(bfreaderTmp, iPlane);
		% figure;
		% imagesc(tmpFrame);
		if numMovies==1
			outputMovie(:,:,iZ) = tmpFrame;
		else
			tmpMovie(:,:,iZ) = tmpFrame;
		end
		iZ = iZ+1;

		if options.displayInfo==1
			reverseStr = ciapkg.view.cmdWaitbar(iZ,nFrames,reverseStr,'inputStr','Loading bio-formats file: ','waitbarOn',options.waitbarOn,'displayEvery',dispEvery);
		end
		% pause(0.01)
	end
end

function [subsets dataDim] = getSubsetOfDataToAnalyze(inputFilePath, options, varargin)
	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	% get HDF5 info
	hinfo = hdf5info(inputFilePath);
	hinfo.GroupHierarchy.Datasets;
	% find dataset name location
	% datasetNames = {hinfo.GroupHierarchy.Datasets.Name};
	% thisDatasetName = strmatch(options.inputDatasetName,datasetNames);
	% hReadInfo = hinfo.GroupHierarchy.Datasets(thisDatasetName);
	hReadInfo = getHdf5Info(hinfo,options);
	dataDim.x = hReadInfo.Dims(1);
	dataDim.y = hReadInfo.Dims(2);
	dataDim.z = hReadInfo.Dims(3);
	% estimate size of movie in Mbytes
	testFrame = readHDF5Subset(inputFilePath,[0 0 0],[dataDim.x dataDim.y 1],'datasetName',options.inputDatasetName);
	testFrameInfo = whos('testFrame');
	estSizeMovie = (testFrameInfo.bytes/options.bytesToMB)*dataDim.z;
	numSubsets = ceil(estSizeMovie/options.maxChunkSize)+1;
	% get the subsets of the 3D matrix to analyze
	subsets = floor(linspace(1,dataDim.z,numSubsets));
end
function hReadInfo = getHdf5Info(hinfo,options)
	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	try
		datasetNames = {hinfo.GroupHierarchy.Datasets.Name};
		thisDatasetName = strmatch(options.inputDatasetName,datasetNames);
		hReadInfo = hinfo.GroupHierarchy.Datasets(thisDatasetName);
	catch
		try
			datasetNames = {hinfo.GroupHierarchy.Groups.Datasets.Name};
			thisDatasetName = strmatch(options.inputDatasetName,datasetNames);
			hReadInfo = hinfo.GroupHierarchy.Groups.Datasets(thisDatasetName);
		catch
			nGroups = length(hinfo.GroupHierarchy.Groups);
			datasetNames = {};
			for groupNo = 1:nGroups
				datasetNames{groupNo} = hinfo.GroupHierarchy.Groups(groupNo).Datasets.Name;
			end
			thisDatasetName = strmatch(options.inputDatasetName,datasetNames);
			hReadInfo = hinfo.GroupHierarchy.Groups(thisDatasetName).Datasets;
		end
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function C = inv(S)
% compliance to stiffness tensor
%
% Input
%  S - @complianceTensor
%
% Output
%  C - @stiffnessTensor
%

C = stiffnessTensor(inv@tensor(S));

end

% this can be done more explicitely by
function test

M = matrix(S,'voigt');

D = M(1,1,:) .* M(2,2,:) .* M(3,3,:) ...
  - M(1,1,:) .* M(2,3,:).^2 ...
  - M(2,2,:) .* M(1,3,:) .* M(1,3,:) ...
  - M(3,3,:) .* M(1,2,:).^2 ...
  + 2 * M(1,2,:) .* M(2,3,:) .* M(1,3,:);

C11 = (M(2,2,:) .* M(3,3,:)-M(2,3,:) .* M(2,3,:))/D;
C12 = (M(1,3,:) .* M(2,3,:)-M(1,2,:) .* M(3,3,:))/D;
C13 = (M(1,2,:) .* M(2,3,:)-M(1,3,:) .* M(2,2,:))/D;
C22 = (M(1,1,:) .* M(3,3,:)-M(1,3,:) .* M(1,3,:))/D;
C23 = (M(1,2,:) .* M(1,3,:)-M(2,3,:) .* M(1,1,:))/D;
C33 = (M(1,1,:) .* M(2,2,:)-M(1,2,:) .* M(1,2,:))/D;

% Enter tensor as 6 by 6 matrix,M line by line.
M = [[  C11   C12   C13    0     0     0];...
    [   C12   C22   C23    0     0     0];...
    [   C13   C23   C33    0     0     0];...
    [   0      0      0   1./M(4,4,:)    0     0];...
    [   0      0      0    0    1./M(5,5,:)    0];...
    [   0      0      0    0     0   1./M(6,6,:)]];
%
C = reshape(stiffnessTensor(M,cs_tensor,'density',rho),size(S));

end

author: mtex-toolbox

repo: mtex

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


text:
function vscl = vscale( f )
%VSCALE   Estimate the vertical scale of a function.
%   VSCALE(F) estimates the vertical scale (also known as the dynamical range)
%   of a function. This is required because a TRIGTECH does not store its
%   interpolation data. If F is an array-valued TRIGTECH with K columns, then
%   the result is a row vector of length K.

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

% Get the coefficients:
c = f.coeffs;

% Check isempty:
if ( isempty( c ) )
    vscl = 0;
elseif ( size(c, 1) == 1 )
    vscl = abs(c);
else
    % Compute values:
    vals = f.coeffs2vals(c);
    % Take max:
    vscl = max(abs(vals), [], 1);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [GF2] = GreenFunction2(r0,r,f0)

parameter; % Define parameters
k0    = 2*pi*f0/c; % fixed wavenumber

GF2   = zeros(1,N+1); % Initialisation of array vector
mode0 = zeros(1,N+1);

% Construct a line from r0 to r
diff_r_r0 = r - r0;    %distance between r and r0
rz0 = [r(1);r(2);0];   %Projection onto x-y plane
diff_rz0_r = rz0 - r0; %distance between r in z-axis and r0

% Calculate the angles between the lines
angle1 = acos((diff_r_r0'*diff_rz0_r)/(norm(diff_r_r0)*norm(diff_rz0_r)));
angle2 = acos((diff_rz0_r'*rz0)/(norm(rz0)*norm(diff_rz0_r)));

Lvec = norm(diff_r_r0);  %length of line
Nvec = (0:1:N);
xpos = Nvec.*Lvec/N*cos(angle1)*sin(angle2) + r0(1);
ypos = Nvec.*Lvec/N*cos(angle1)*cos(angle2) + r0(2);
zpos = Nvec.*Lvec/N*sin(angle1)             + r0(3);

for d = 1:N+1

%%% Calculate the AXIAL MODE %%%
    %lx
    	N = fix(2*2*lx*f0/c); %number of axial modes below frequency f0
        nx = 1:N;
        rx = xpos(d);
    	[modeshape,k_term] = AxialMode(nx,lx,rx,r0(1),k0);
    	GF2(d) = GF2(d) + (-1/V)*sum(sum(modeshape./k_term));
    %ly
        N = fix(2*2*ly*f0/c); %number of axial modes below frequency f0
        ny = 1:N;
        ry = ypos(d);
    	[modeshape,k_term] = AxialMode(ny,ly,ry,r0(2),k0);
    	GF2(d) = GF2(d) + (-1/V)*sum(sum(modeshape./k_term));
    %lz
        N = fix(2*2*lz*f0/c); %number of axial modes below frequency f0
        nz = 1:N;
        rz = zpos(d);
    	[modeshape,k_term] = AxialMode(nz,lz,rz,r0(3),k0);
    	GF2(d) = GF2(d) + (-1/V)*sum(sum(modeshape./k_term));

%%% Calculate the TANGENTIAL MODE %%%
    %lx,ly
    	N = fix(2*pi*lx*ly*f0^2/c^2); %number of tangential modes below frequency f0
    	nx = 1:N; ny = nx;
        rx = xpos(d); ry = ypos(d);
    	[modeshape,k_term] = TangentialMode(nx,ny,lx,ly,rx,ry,r0(1),r0(2),k0);
    	GF2(d) = GF2(d) + (-1/V)*sum(sum(modeshape./k_term));
    %ly,lz
		N = fix(2*pi*ly*lz*f0^2/c^2); %number of tangential modes below frequency f0
        ny = 1:N; nz = ny;
        ry = ypos(d); rz = zpos(d);
    	[modeshape,k_term] = TangentialMode(ny,nz,ly,lz,ry,rz,r0(2),r0(3),k0);
    	GF2(d) = GF2(d) + (-1/V)*sum(sum(modeshape./k_term));
    %lx,lz
    	N = fix(2*pi*lx*lz*f0^2/c^2); %number of tangential modes below frequency f0
        nx = 1:N; nz = nx;
        rx = xpos(d); rz = zpos(d);
    	[modeshape,k_term] = TangentialMode(nx,nz,lx,lz,rx,rz,r0(1),r0(3),k0);
    	GF2(d) = GF2(d) + (-1/V)*sum(sum(modeshape./k_term));

%%% Calculate the OBLIQUE MODE %%%
    %lx,ly,lz
    	N = fix(2*4*pi*lx*ly*lz*f0^3/(3*c^3)); %number of oblique modes below frequency f0
        ny = 1:N; nz = ny;
        rx = xpos(d); ry = ypos(d); rz = zpos(d);
        for nx = 1:N,
            [modeshape,k_term] = ObliqueMode(nx,ny,nz,lx,ly,lz,rx,ry,rz,r0(1),r0(2),r0(3),k0);
            GF2(d) = GF2(d) + (-1/V)*sum(sum(modeshape./k_term));
        end;

%%% Calculate the (0,0,0) mode %%%
    mode0 = -8/(V*k0.^2);

%%% Calculate the Green function plus the (0,0,0) mode
    GF2(d) = GF2(d) + mode0;
end;

author: Sable

repo: mcbench-benchmarks


text:
function [PB] = bit2PB(bit)
% Convert computery things from bits to petabytes.
% Chad A. Greene 2012
PB = bit*2^-53;


text:
%PREX_DATAFILE  PRTools example of the datafile usage

delfigs
echo on

prdatafiles
            % makes sure that datafiles are available
						% will download them from the PRTools website when needed

a = highway
            % load datafile
            % 100 observations of 5 images and a label image
						% R,G,B
						% and two pixel features (comparisons with previous frames)

b = gendat(a,0.06)
            % random selection of 6 observations

figure; show(b); drawnow

c = selectim(b,[1 2 3])
            % select RGB

figure; show(c); drawnow
showfigs

x = b(2,:)
           % select one observation
figure; show(x,3); drawnow
showfigs

y = data2im(x);
           % select features by retrieving images

data = squeeze(y(:,:,1:5)); % data
lab  = squeeze(y(:,:,6));   % labels
datasize = size(data);

z = im2feat(data);  % store images as feature, pixels become objects
z = setlabels(z,lab(:));
z = setobjsize(z,datasize(1:2));
            % stores pixels as objects with 5 features

trainset = gendat(z,[1000,1000]);
w = qdc(trainset) % train classifier on trainset only
d = z*w*classc;   % classify entire image
figure; show(z*w*classc);
figure; imagesc(d*classim);
figure; imagesc(lab);
showfigs

echo off

author: marianux

repo: ecg-kit

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


text:
function [iHi,iLo,iCr] = findextrema(x)
%FINDEXTREMA find indices of local extrema and zero-crossings.
%   [IMAX,IMIN,ICRS] = FINDEXTREMA(X) returns the indices of local maxima
%   in IMAX, minima in IMIN and zero-crossing in ICRS for input vector X.
%
%   Example:
%       x = [0 1 0 0 -1 -1 -2 -2 1 0];
%       [i1,i2,i3] = findextrema(x)
%       % Returns:
%       %   i1 = 2 9
%       %   i2 = 8
%       %   i3 = 1 4 8 10
%
%   See also FIND.

% Siyi Deng; 05-29-2009;
% sdeng@uci.edu; UCI HNL;

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

%% Generate Plot for demo;
% x = [1 2 1 3 3 3 3 4 4 4 4 3 3 3 2 2 2 1 3 3 2 2 6 6 6 5 5 4 4 4 1 1 2]-3;
% figure; plot(x); set(gca,'xtick',0:numel(x)+2); ylim([-4 4])
% [iP,iT,iC] = findextrema(x),
% hold on; plot(iP,x(iP),'r*'); plot(iT,x(iT),'g*'); plot(iC,x(iC),'ks');
% legend({'Data','Peak','Trough','Zero-Crossing'},'location','NorthWest');

%% Code starts here;
if ~isvector(x), error('Input must be a vector;'); end
dx = diff(x);
gz = dx > 0;
lz = dx < 0;
ez = dx == 0;
hi = gz(1:end-1) & lz(2:end);
hasCorner = any(ez);
if hasCorner
    cornerVec = double(gz(1:end-1) & ez(2:end))+...
        double(ez(1:end-1) & lz(2:end)).*2+...
        double(ez(1:end-1) & gz(2:end)).*110+...
        double(lz(1:end-1) & ez(2:end)).*100;
    cornerLoc = find(cornerVec);
    cornerType = diff(cornerVec(cornerLoc));
    n = find(cornerType == 1); % plateu;
    hi(ceil((cornerLoc(n+1)+cornerLoc(n))/2)) = true;
end
iHi = find(hi)+1;
if nargout > 1
    lo = gz(2:end) & lz(1:end-1);
    if hasCorner
        u = find(cornerType == 10); % valley;
        lo(ceil((cornerLoc(u+1)+cornerLoc(u))/2)) = true;
    end
    iLo = find(lo)+1;
end
if nargout > 2
    xc = (x(1:end-1).*x(2:end)) < 0;
    xz = false(length(x)+2,1);
    xz(2:end-1) = x == 0;
    if any(xz)
        xc(fix((find(~xz(1:end-1) & xz(2:end))+...
            find(~xz(2:end) & xz(1:end-1)))/2)) = true;
    end
    iCr = find(xc);
end

end % FINDEXTREMA;


text:
classdef GenericLoss < dagnn.ElementWise

  properties
    display = true;
  end

  properties (Transient)
    average = 0
    numAveraged = 0
  end

  methods
    function outputs = account(obj, inputs, outputs)
      if obj.display
        n = obj.numAveraged ;
        m = n + size(inputs{1},4) ;
        obj.average = (n * obj.average + gather(outputs{1})) / m ;
        obj.numAveraged = m ;
      end
    end

    function reset(obj)
      obj.average = 0 ;
      obj.numAveraged = 0 ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes, paramSizes)
      outputSizes{1} = [1 1 1 inputSizes{1}(4)] ;
    end

    function rfs = getReceptiveFields(obj)
      % the receptive field depends on the dimension of the variables
      % which is not known until the network is run
      rfs(1,1).size = [NaN NaN] ;
      rfs(1,1).stride = [NaN NaN] ;
      rfs(1,1).offset = [NaN NaN] ;
      rfs(2,1) = rfs(1,1) ;
    end

  end
end

author: a-nagrani

repo: VGGVox

save_path: github-repos/MATLAB/a-nagrani-VGGVox


text:
function y = nanmedian(x,dim)
%NANMEDIAN - Median value, ignoring NaN values.
%
%See median

if nargin==1,
  dim = min(find(size(x)~=1));
  if isempty(dim), dim = 1; end
end
if isempty(x), y = []; return, end

siz= size(x);
n= siz(dim);

%% Permute and reshape so that DIM becomes the row dimension of a 2-D array
perm= [dim:max(length(size(x)),dim) 1:dim-1];
x= reshape(permute(x,perm),n,prod(siz)/n);

%% Do it columnwise (unefficient)
y= zeros(1, size(x,2));
for ii= 1:size(x,2),
  idx= find(~isnan(x(:,ii)));
  if isempty(idx),
    y(ii)= NaN;
  else
    y(ii)= median(x(idx,ii));
  end
end

%% Permute and reshape back
siz(dim)= 1;
y= ipermute(reshape(y,siz(perm)),perm);

author: bbci


text:
% MAP_PROJECTION - Set the projection
%
% 'GLOBAL-ELLIPSE'
% 'GLOBAL-RECT'
% 'TESTBED'

function map_projection(projection)

if nargin < 1
  projection = 'global-ellipse';
end

switch lower(projection)
 case 'global-ellipse'
  m_proj('mollweide', 'clongitude', 0); % straight latitudes (equal area)
  %m_proj('hammer-aitoff', 'clon', 0); % curved latitudes (equal area)
 case 'global-rect'
  m_proj('miller');
 case 'testbed'
  m_proj('Mercator', 'lon',[22.5,26.8], 'lat',[59.74,61.0])
 otherwise
  error('Unknown projection')
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function fem1d_bvp_linear_test02 ( )

%*****************************************************************************80
%
%% FEM1D_BVP_LINEAR_TEST02 carries out test case #2.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/fem1d_bvp_linear/fem1d_bvp_linear_test02.m
%
%  Discussion:
%
%    Use A2, C2, F2, EXACT2, EXACT_UX2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dianne O'Leary,
%    Scientific Computing with Case Studies,
%    SIAM, 2008,
%    ISBN13: 978-0-898716-66-5,
%    LC: QA401.O44.
%
  n = 11;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_BVP_LINEAR_TEST02\n' );
  fprintf ( 1, '  Solve -( A(x) U''(x) )'' + C(x) U(x) = F(x)\n' );
  fprintf ( 1, '  for 0 < x < 1, with U(0) = U(1) = 0.\n' );
  fprintf ( 1, '  A2(X)  = 1.0\n' );
  fprintf ( 1, '  C2(X)  = 2.0\n' );
  fprintf ( 1, '  F2(X)  = X * ( 5 - X ) * exp ( X )\n' );
  fprintf ( 1, '  U2(X)  = X * ( 1 - X ) * exp ( X )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes = %d\n', n );
%
%  Geometry definitions.
%
  x_first = 0.0;
  x_last = 1.0;
  x = linspace ( x_first, x_last, n );
  x = x(:);

  u = fem1d_bvp_linear ( n, @a2, @c2, @f2, x );

  uexact = exact2 ( x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     I    X         U         Uexact    Error\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %4d  %8f  %8f  %8f  %8e\n', ...
      i, x(i), u(i), uexact(i), abs ( u(i) - uexact(i) ) );
  end

  e1 = l1_error ( n, x, u, @exact2 );
  e2 = l2_error_linear ( n, x, u, @exact2 );
  h1s = h1s_error_linear ( n, x, u, @exact_ux2 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  l1 norm of error  = %g\n', e1 );
  fprintf ( 1, '  L2 norm of error  = %g\n', e2 );
  fprintf ( 1, '  Seminorm of error = %g\n', h1s  );

  return
end
function value = a2 ( x )

%*****************************************************************************80
%
%% A2 evaluates A function #2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of A(X).
%
  value = 1.0;

  return
end
function value = c2 ( x )

%*****************************************************************************80
%
%% C2 evaluates C function #2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of C(X).
%
  value = 2.0;

  return
end
function value = exact2 ( x )

%*****************************************************************************80
%
%% EXACT2 evaluates exact solution #2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of U(X).
%
  value = x .* ( 1.0 - x ) .* exp ( x );

  return
end
function value = exact_ux2 ( x )

%*****************************************************************************80
%
%% EXACT_UX2 evaluates the derivative of exact solution #2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of dUdX(X).
%
  value = ( 1.0 - x - x .* x ) .* exp ( x );

  return
end
function value = f2 ( x )

%*****************************************************************************80
%
%% F2 evaluates right hand side function #2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of F(X).
%
  value = x .* ( 5.0 - x ) .* exp ( x );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = grDilate(varargin)
%GRDILATE Morphological dilation on graph.
%
%   LBL2 = grDilate(EDGES, LBL1)
%   Each label of the graph is assigned the highest label of its
%   neighbours, or it keeps the same label this one is bigger.
%
%   Example
%   grDilate
%
%   See also
%   grErode, grOpen, grClose
%

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2006-01-20
% Copyright 2006-2022 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas)

if length(varargin) == 2
    edges   = varargin{1};
    lbl     = varargin{2};
elseif length(varargin) == 3
    edges   = varargin{2};
    lbl     = varargin{3};
else
    error('Wrong number of arguments in "grDilate"');
end

lbl2 = zeros(size(lbl));

uni = unique(edges(:));
for n = 1:length(uni)
    neigh = grAdjacentNodes(edges, uni(n));
    lbl2(uni(n)) = max(lbl([uni(n); neigh]));
end

varargout{1} = lbl2;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function sphere = makeSphere(Nx, Ny, Nz, radius, plot_sphere, binary)
%MAKESPHERE     Create a binary map of a sphere within a 3D grid.
%
% DESCRIPTION:
%       makeSphere creates a binary map of a spherical shell (using an
%       extension of the midpoint circle algorithm) within a
%       three-dimensional grid. The sphere position is denoted by 1's in
%       the matrix with 0's elsewhere. If the Boolean input parameter
%       "binary" is set to false (the default), the sphere map is returned
%       as a double precision matrix. If it is set to true, the map is
%       returned as a logical matrix.
%
% USAGE:
%       sphere = makeSphere(Nx, Ny, Nz, radius)
%       sphere = makeSphere(Nx, Ny, Nz, radius, plot_sphere)
%       sphere = makeSphere(Nx, Ny, Nz, radius, plot_sphere, binary)
%       sphere = makeSphere(Nx, Ny, Nz, radius, [], binary)
%
% INPUTS:
%       Nx, Ny, Nz      - size of the 3D grid [grid points]
%       radius          - sphere radius [grid points]
%
% OPTIONAL INPUTS:
%       plot_sphere     - Boolean controlling whether the sphere is
%                         plotted using voxelPlot (default = false)
%       binary          - Boolean controlling whether the sphere map is
%                         returned as a double precision matrix (false) or
%                         a logical matrix (true) (default = false)
%
% OUTPUTS:
%       sphere          - 3D binary map of a sphere
%
% ABOUT:
%       author          - Bradley Treeby
%       date            - 16th June 2009
%       last update     - 20th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also makeBall, makeCartSphere, makeCircle, makeSphericalSection

% This file is part of k-Wave. k-Wave 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.
%
% k-Wave 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 Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% check for plot_sphere input
if nargin < 5 || isempty(plot_sphere)
    plot_sphere = false;
end

% check for binary input
if nargin < 6 || isempty(binary)
    binary = false;
end

% enforce a centered sphere
cx = floor(Nx/2)+1;
cy = floor(Ny/2)+1;
cz = floor(Nz/2)+1;

% preallocate the storage variable
if binary
    sphere = false(Nx, Ny, Nz);
else
    sphere = zeros(Nx, Ny, Nz);
end

% create a guide circle from which the individal radii can be extracted
guide_circle = makeCircle(Ny, Nx, cy, cx, radius);

% step through the guide circle points and create partially filled discs
centerpoints = (cx - radius):cx;
reflection_offset = length(centerpoints):-1:2;
for centerpoint_index = 1:length(centerpoints)

    % extract the current row from the guide circle
    row_data = guide_circle(:, centerpoints(centerpoint_index));

    % add an index to the grid points in the current row
    row_index = row_data.*(1:length(row_data)).';

    % calculate the radius
    swept_radius = (max(row_index) - min(row_index(row_index ~= 0)))/2;

    % create a circle to add to the sphere
    circle = makeCircle(Ny, Nz, cy, cz, swept_radius);

    % make an empty fill matrix
    if binary
        circle_fill = false(Ny, Nz);
    else
        circle_fill = zeros(Ny, Nz);
    end

    % fill in the circle line by line
    fill_centerpoints = (cz - swept_radius):(cz + swept_radius);
    for fill_centerpoint_index = 1:length(fill_centerpoints)

        % extract the first row
        row_data = circle(:, fill_centerpoints(fill_centerpoint_index));

        % add an index to the grid points in the current row
        row_index = row_data.*(1:length(row_data)).';

        % calculate the diameter
        start_index = min(row_index(row_index ~= 0));
        stop_index = max(row_index);

        % count how many points on the line
        num_points = sum(row_data);

        % fill in the line
        if start_index ~= stop_index && (stop_index - start_index) >= num_points
            circle_fill(start_index + num_points/2:stop_index - num_points /2, fill_centerpoints(fill_centerpoint_index)) = 1;
        end
    end

    % remove points from the filled circle that existed in the previous
    % layer
    if centerpoint_index == 1
        sphere(centerpoints(centerpoint_index), :, :) = circle + circle_fill;
        prev_circle = circle + circle_fill;
    else
        prev_circle_alt = circle + circle_fill;
        circle_fill = circle_fill - prev_circle;
        circle_fill(circle_fill < 0) = 0;
        sphere(centerpoints(centerpoint_index), :, :) = circle + circle_fill;
        prev_circle = prev_circle_alt;
    end

    % create the other half of the sphere at the same time
    if centerpoint_index ~= length(centerpoints)
        sphere(cx + reflection_offset(centerpoint_index) - 1, :, :) = sphere(centerpoints(centerpoint_index), :, :);
    end
end

% plot results
if plot_sphere
    voxelPlot(double(sphere));
end

repo: Aero-matlab


text:
function [ output_args ] = CudaConvolution( input_args )
% CUDACONVOLUTION Download description
% Please download the
%
% 64 Bit CUDA 4.2 MEX files here.
%
% and download the
%
% 32 Bit CUDA 4.2 MEX files here.
%
% and download the
%
% NVIDIA driver version 190.38 for your OS from here.
%
realKernel = double(1:5);
imagKernel = double(1:5);
realKernel(:)= 0.127;
imagKernel(:)= 0;

% Construct the filter kernel
Filter_Kernel = complex(realKernel,imagKernel);

imagSignal = double(1:2^15);

time = 0:(1/2^15):(1-1/2^15);
idealSignal = sin(2*pi*5*time);

realSignal = idealSignal + randn(size(time))/2;
imagSignal(:) = 0;

% Construct the signal
Signal_1D_Complex = complex(realSignal,imagSignal);

% Convolution of kernel and signal on the CUDA device
filteredSignal = cuFilter(Signal_1D_Complex,Filter_Kernel);

plot(realSignal,'b');hold;
plot(real(filteredSignal),'m');
plot(idealSignal,'g');


text:
function kappa = polynomialCurveCurvature(t, varargin)
%POLYNOMIALCURVECURVATURE Compute the local curvature of a polynomial curve
%
%   KAPPA = polynomialCurveCurvature(T, XCOEF, YCOEF)
%   XCOEF and YCOEF are row vectors of coefficients, in the form:
%       [a0 a1 a2 ... an]
%   KAPPA is the local curvature of the polynomial curve, computed for
%   position T. If T is a vector, KAPPA has the same length as T.
%
%   KAPPA = polynomialCurveCurvature(T, COEFS)
%   COEFS is either a 2xN matrix (one row for the coefficients of each
%   coordinate), or a cell array.
%
%   Example
%   polynomialCurveCurvature
%
%   See also
%   polynomialCurves2d, polynomialCurveLength, polynomialCurveDerivative
%
%
% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2007-02-23
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

%% Extract input parameters

% polynomial coefficients for each coordinate
var = varargin{1};
if iscell(var)
    xCoef = var{1};
    yCoef = var{2};
elseif size(var, 1)==1
    xCoef = varargin{1};
    yCoef = varargin{2};
else
    xCoef = var(1,:);
    yCoef = var(2,:);
end

%% compute derivative

% compute first derivatives of the polynomials
dx  = polynomialDerivate(xCoef);
dy  = polynomialDerivate(yCoef);

% compute second derivatives
sx  = polynomialDerivate(dx);
sy  = polynomialDerivate(dy);

% convert to polyval convention
dx  = dx(end:-1:1);
dy  = dy(end:-1:1);
sx  = sx(end:-1:1);
sy  = sy(end:-1:1);

% compute local first and second derivatives
xp  = polyval(dx, t);
yp  = polyval(dy, t);
xs  = polyval(sx, t);
ys  = polyval(sy, t);

% compute local curvature of polynomial curve
kappa  = (xp.*ys - xs.*yp) ./ power(xp.*xp + yp.*yp, 3/2);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function a = pds_random ( n, key )

%*****************************************************************************80
%
%% PDS_RANDOM returns a random positive definite symmetric matrix.
%
%  Discussion:
%
%    The matrix returned will have eigenvalues in the range [0,1].
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%    A is positive definite: 0 < x'*A*x for nonzero x.
%
%    The eigenvalues of A will be real.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, integer KEY, a positive value that selects the data.
%
%    Output, real A(N,N), the matrix.
%

%
%  Get a random set of eigenvalues.
%
  seed = key;
  [ lambda, seed ] = r8vec_uniform_01 ( n, seed );
%
%  Get a random orthogonal matrix Q.
%
  q = orth_random ( n, key );
%
%  Set A = Q * Lambda * Q'.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n
      a(i,j) = 0.0;
      for k = 1 : n
        a(i,j) = a(i,j) + q(i,k) * lambda(k) * q(j,k);
      end
    end
  end

  return
end


text:
%% Example
% real time control of the KUKA iiwa 7 R 800
% the impedence control is turned on
% Moving first joint of the robot, using a sinisoidal function

% The external torques are plotted in real-time during the test
% the feedback from the measured torques can be used to for a closed loop
% control

% An example script, it is used to show how to use the different
% functions of the KUKA Sunrise matlab toolbox

% First run the following script in Matlab
% Then start the client on the KUKA iiwa controller

% Mohammad SAFEEA, 24th of April 2018

close all,clear all;clc;

warning('off')

ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
global t_Kuka;
t_Kuka=net_establishConnection( ip );

if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
  warning('Connection could not be establised, script aborted');
  return;
else

      %% Move point to point to an initial position
        jPos={0,0,0,-pi/2,0,pi/2,0};
      setBlueOff(t_Kuka); % turn Off blue light

      relVel=0.15;
      movePTPJointSpace( t_Kuka , jPos, relVel); % move to initial configuration
     %% Pause for 3 seocnds
     pause(3);
        %% Start direct servo in joint space
        massOfTool=0.0; % the mass of the tool attached to flange in Kg
        cOMx=0; % X coordinate of the center of mass of the tool in (mm)
        cOMy=0; % Y coordinate of the center of mass of the tool in (mm)
        cOMz=40; % Z coordinate of the center of mass of the tool in (mm)
        cStiness=900; % cartizian stifness
        rStifness=80; % rotational stifness
        nStifness=50; % null space stifness

        % Start the realtime control with impedence
        realTime_startImpedanceJoints(t_Kuka,massOfTool,cOMx,cOMy,cOMz,...
        cStiness,rStifness,nStifness);

       w=0.6; % motion constants, frequency rad/sec
       A=0.2; % motion constants, amplitude of motion

       a=datevec(now);
       t0=a(6)+a(5)*60+a(4)*60*60; % calculate initial time

     dt=0;
     tstart=t0;
     counter=0;
     duration=20*60; %1 minutes
     %% Control loop
     % real time plot handles
         colors={'k','b','r','g','c','m','y'};
         figureHandle=figure('Units','inches','Position',[0 0 5 3.75]);
         numOfSamples=120; % number of samples to show in the plot
         timeVec=1:numOfSamples;
         tawVec=zeros(7,numOfSamples);
         plotHandle=[];
         for i=1:7
            plotHandle=[plotHandle,plot(timeVec,tawVec(i,:),colors{i},'LineWidth',2)];
            hold on;
         end
        % format the plot
        ylim([0,8]);
        temp=title('External torques');
        set(temp,'FontSize',16);
        temp=xlabel('Time (seconds)');
        set(temp,'FontSize',14);
        temp=ylabel('Distance (m)');
        set(temp,'FontSize',14);

       while(dt<duration)
         %% perform trajectory calculation here
          a=datevec(now);
          time=a(6)+a(5)*60+a(4)*60*60;
          dt=time-t0;

          temp=A*(1-cos(w*dt));
            for dacount=1:7
                jPosCommand{dacount}=jPos{dacount}+temp;
            end
          counter=counter+1;
          %% Send joint positions to robot
          taw=sendJointsPositionsExTorque( t_Kuka ,jPosCommand);

          tawTemp=zeros(7,1);
          for i=1:7
              tawTemp(i)=taw{i};
              torque_array(i,counter)=taw{i};
          end
          % delete first measurment
          tawVec(:,1)=[];
          % add new measurment to end of arrays
          tawVec=[tawVec,tawTemp];
          % update plots data
          for i=1:6
            set(plotHandle(i),'YData',tawVec(i,:));
          end
          set(plotHandle(7),'YData',tawVec(i,:));
          % show plot data
          set(figureHandle,'Visible','on');

       end

       %% Stop the realtime control with impedence motion
       realTime_stopImpedanceJoints( t_Kuka );

       fprintf('\n Motion stopped \n');
       pause(2);
       %% turn off light
       setBlueOff(t_Kuka);

      %% turn off the server
       net_turnOffServer( t_Kuka );

       fclose(t_Kuka);

end
warning('on')

repo: KST-Kuka-Sunrise-Toolbox


%MVDR_FF: Minimum Variance Distortionless Response Beamformer using
%Free-field Steering Vector
%--------------------
% Input parameters
%--------------------
%nChannels :   The steering vector, or relative transfer function [numChannels x numFreqs ]
%dMic : The distance between two microphones [a real positive value]
%DOA: The direction of arrival in degree [ 0~360 ]
%freqVec: The frequencies vector [ numFreqs]
%Phi_u: The PSD matrix of the interfering signal [numChannels x numChannels x numFreqs ]
%--------------------
% Output parameters
%--------------------
%H_mvdr: An mvdr beamformer [ numChannels x numFreqs ]
function H_mvdr = mvdrFF(nChannels, dMic, DOA, freqVec, Phi_u)

c = 340;
theta0 = DOA * pi / 180;
d = exp(-1j*[0:nChannels-1]'*2*pi*freqVec*cos(theta0)*dMic/c) ./ nChannels;  % steering vector

H_mvdr = zeros(nChannels, length(freqVec));
%Phi_u_inv = zeros(nChannels, nChannels);
for f = 1:length(freqVec)
   Phi_u_inv = pinv(Phi_u(:,:,f));
   df = d(:, f);
   H_mvdr(:, f) = (Phi_u_inv * df) / (df' * Phi_u_inv * df);
end

end


text:
function [yd3] = gal2yd3(gal)
% Convert volume from US liquid gallons to cubic yards.
% Chad Greene 2012
yd3 = gal*0.0049511316873;


text:
function sysmodel = SparseRegression(Y,Yp,U,dt, options_method, lambda)

Nstates = size(Y,1);
Ninputs = size(U,1);

% Sparse regression
switch options_method.sparsify
    case 'LASSO'
        % using the l1 norm to promote sparsity
        G = zeros(Nstates+Ninputs,Nstates);
        for i = 1:Nstates
            [G(:,i), FitInfo] = lasso([Y;U]',Yp(i,:)','Lambda',lambda); %,'Lambda',0.00000001
            %     lassoPlot(G{i},FitInfo);
        end
    case 'ILSTH'
        G = sparsifyDynamics([Y;U]',Yp',lambda,Nstates);
end

G = G';
A = G(1:Nstates,1:Nstates);
B = G(1:Nstates,Nstates+1:Nstates+Ninputs);
C = eye(Nstates,Nstates);
D = zeros(Nstates,1);
sysmodel = ss(A,B,C,D,dt);
end

author: eurika-kaiser

repo: SINDY-MPC


text:
function cond = diagonal_condition ( n, x )

%*****************************************************************************80
%
%% DIAGONAL_CONDITION returns the L1 condition of the DIAGONAL matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real X(N), the diagonal entries of A.
%
%    Output, real COND, the L1 condition.
%
  cond = max ( abs ( x(1:n) ) ) / min ( abs ( x(1:n) ) );

  return
end


text:
function [ n_data, n, z, fx ] = polylogarithm_values ( n_data )

%*****************************************************************************80
%
%% POLYLOGARITHM_VALUES returns some values of the polylogarithm.
%
%  Discussion:
%
%    The polylogarithm of n and z is defined as
%
%      f[n,z] = Sum ( 1 <= k < infinity ) z^k / k^n
%
%    In Mathematica, the function can be evaluated by:
%
%      PolyLog[n,z]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the exponent of the denominator.
%
%    Output, real Z, the base of the numerator.
%
%    Output, real FX, the value of the function.
%
  n_max = 12;

  fx_vec = [ ...
     0.1644934066848226E+01, ...
     0.1202056903159594E+01, ...
     0.1000994575127818E+01, ...
     0.5822405264650125E+00, ...
     0.5372131936080402E+00, ...
     0.5002463206060068E+00, ...
     0.3662132299770635E+00, ...
     0.3488278611548401E+00, ...
     0.3334424797228716E+00, ...
     0.1026177910993911E+00, ...
     0.1012886844792230E+00, ...
     0.1000097826564961E+00 ];

  n_vec = [ ...
     2, 3, 10, 2, 3, 10, 2, 3, 10, 2, 3, 10 ];

  z_vec = [ ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.5000000000000000E+00, ...
     0.5000000000000000E+00, ...
     0.5000000000000000E+00, ...
     0.3333333333333333E+00, ...
     0.3333333333333333E+00, ...
     0.3333333333333333E+00, ...
     0.1000000000000000E+00, ...
     0.1000000000000000E+00, ...
     0.1000000000000000E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    z = 0.0;
    fx = 0.0;
  else
    n = n_vec(n_data);
    z = z_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function [val, remaining] = keyval(key, varargin)

% KEYVAL returns the value that corresponds to the requested key in a
% key-value pair list of variable input arguments
%
% Use as
%   [val] = keyval(key, varargin)
%
% See also VARARGIN

% Undocumented option
%   [val] = keyval(key, varargin, default)

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

% what to return if the key is not found
emptyval = [];

if nargin==3 && iscell(varargin{1})
  emptyval = varargin{2};
  varargin = varargin{1};
end

if nargin==2 && iscell(varargin{1})
  varargin = varargin{1};
end

if mod(length(varargin),2)
  error('optional input arguments should come in key-value pairs, i.e. there should be an even number');
end

% the 1st, 3rd, etc. contain the keys, the 2nd, 4th, etc. contain the values
keys = varargin(1:2:end);
vals = varargin(2:2:end);

% the following may be faster than cellfun(@ischar, keys)
valid = false(size(keys));
for i=1:numel(keys)
  valid(i) = ischar(keys{i});
end

if ~all(valid)
  error('optional input arguments should come in key-value pairs, the optional input argument %d is invalid (should be a string)', i);
end

hit = find(strcmpi(key, keys));
if isempty(hit)
  % the requested key was not found
  val = emptyval;
elseif length(hit)==1
  % the requested key was found
  val = vals{hit};
else
  error('multiple input arguments with the same name');
end

if nargout>1
  % return the remaining input arguments with the key-value pair removed
  keys(hit) = [];
  vals(hit) = [];
  remaining = cat(1, keys(:)', vals(:)');
  remaining = remaining(:)';
end

author: PatternRecognition

repo: OpenBMI


text:
function [output, creationString] = prtUtilObjectGuiSimple(obj)
% PRTUTILOBJECTGUISIMPLE Provides a graphical method to change an object
%
% Syntax: prtUtilStructureGui(BaseStructure,DefaultStructure,DefinitionStructure)
%
% Internal
% xxx Need Help xxx

% Numberic types bounds etc do not work
newFunctionOutput = [];

% Create default creationString
if isstruct(obj)
    guiPropNames = fieldnames(obj);
    name = inputname(1);
else
    thisClass = class(obj);
    m = meta.class.fromName(thisClass);
    name = thisClass;

    isInThisClass = cellfun(@(c)strcmp(c.DefiningClass.Name,thisClass),m.Properties);
    isSetable = cellfun(@(c)strcmp(c.SetAccess,'public'),m.Properties);
    isHidden = cellfun(@(c)c.Hidden,m.Properties);
    provideInGui = isInThisClass & ~isHidden & isSetable;

    GuiProperties = m.Properties(provideInGui);

    nGuiProperties = length(GuiProperties);
    % Create default structure

    if nGuiProperties == 0
        msgbox(sprintf('%s has no properties that can be set using this GUI.',name),'No Properties','modal');
        output = obj;
        creationString = sprintf('%s()',class(obj));

        return
    end

    guiPropNames = cellfun(@(c)c.Name,GuiProperties,'uniformoutput',false);
end

guiPropValues = cell(nGuiProperties,1);
for iProp = 1:nGuiProperties
    guiPropValues{iProp} = obj.(guiPropNames{iProp});
end

% Determine types of each field in DefaultStructure
%   This is rather difficult to do robustly
%   The if else try below is just hueristic
types = cell(length(guiPropValues),1);
for iProp = 1:length(guiPropValues)
    cVal = guiPropValues{iProp};

    if isnumeric(cVal) || islogical(cVal)
        types{iProp} = 'numeric';

    elseif ischar(cVal)
        types{iProp} = 'char';

    elseif isstruct(cVal);
        types{iProp} = 'structure';

    elseif isa(cVal,'function_handle')
        types{iProp} = 'function';

    else
        % Not a default MATLAB data type
        % Use the structure version to iteratively call this function
        types{iProp} = 'structure';

    end
end

% At this point we can't handle sub-objects or sub-structures
% So we remove them from the GUI
removeFields = strcmpi(types,'structure');
guiPropNames = guiPropNames(~removeFields);
guiPropValues = guiPropValues(~removeFields);
types = types(~removeFields);

nGuiProperties = length(guiPropValues);

structInputs = cell(nGuiProperties*2,1);
structInputs(1:2:end) = guiPropNames;
structInputs(2:2:end) = guiPropValues;
DefaultStructure = struct(structInputs{:});

OutputStructure = DefaultStructure;
BaseStructure = DefaultStructure;

%backgroundColor = get(0,'defaultFigureColor');
backgroundColor = [0.95 0.95 0.95];

DefaultMainPanelOptions.style = 'panel';
DefaultMainPanelOptions.units = 'normalized';
DefaultMainPanelOptions.position = [0.01 0.1 0.98 0.89];
DefaultMainPanelOptions.title = '';
DefaultMainPanelOptions.borderType = 'none';
DefaultMainPanelOptions.backgroundColor = backgroundColor;

MainPanelOptions = DefaultMainPanelOptions;

% Default MainFigureOptions
ss = get(0,'screensize');

windowSize = [754 600];

% Center the window
sizePads = round((ss(3:4)-windowSize));
sizePads(1) = sizePads(1)/2; % We should use 2 right?
sizePads(2) = sizePads(2)/2;
pos = cat(2,sizePads,windowSize);

DefaultMainFigureOptions.style = 'figure';
DefaultMainFigureOptions.units = 'pixels';
DefaultMainFigureOptions.position = pos;
DefaultMainFigureOptions.menuBar = 'none';
DefaultMainFigureOptions.toolbar = 'none';
DefaultMainFigureOptions.numberTitle = 'off';
DefaultMainFigureOptions.name = cat(2,'PRT Action Editor - ',name);
DefaultMainFigureOptions.color = backgroundColor;
DefaultMainFigureOptions.windowStyle = 'modal';
DefaultMainFigureOptions.CloseRequestFcn = @cancelCallback;
DefaultMainFigureOptions.backingstore = 'off';
DefaultMainFigureOptions.DoubleBuffer = 'off';

MainFigureOptions = DefaultMainFigureOptions;

% Finish making the primary elements of the gui
MainFigureOptions.Children.MainPanel = MainPanelOptions;

buttonHeight = 0.06;
MainFigureOptions.Children.Apply.style = 'pushbutton';
MainFigureOptions.Children.Apply.units = 'normalized';
MainFigureOptions.Children.Apply.position = [0.2550 0.02 0.23 buttonHeight];%[0.1 0.02 0.23 0.06];
MainFigureOptions.Children.Apply.string = 'Accept';
MainFigureOptions.Children.Apply.callback = @applyCallback;

MainFigureOptions.Children.Defaults.style = 'pushbutton';
MainFigureOptions.Children.Defaults.units = 'normalized';
MainFigureOptions.Children.Defaults.position = [0.5350 0.02 0.23 buttonHeight];%[0.38 0.02 0.23 0.06];
MainFigureOptions.Children.Defaults.string = 'Restore Defaults';
MainFigureOptions.Children.Defaults.callback = @defaultCallback;

% MainFigureOptions.Children.Cancel.style = 'pushbutton';
% MainFigureOptions.Children.Cancel.units = 'normalized';
% MainFigureOptions.Children.Cancel.position = [0.66 0.02 0.23 0.06];
% MainFigureOptions.Children.Cancel.string = 'Cancel';
% MainFigureOptions.Children.Cancel.callback = @cancelCallback;

%%

GuiInfo.usableVerticalRange = [0.1 0.85];
GuiInfo.label.horizontal = [0.01 0.32];
GuiInfo.value.horizontal = [0.34 0.65];
GuiInfo.edit.horizontal = [0.67 0.99];

GuiInfo.maxHeight = 0.05;

%% Create the field boxes

fields = fieldnames(BaseStructure);
nFields = length(fields);

itemHeight = min((GuiInfo.usableVerticalRange(2)-GuiInfo.usableVerticalRange(1))./(nFields+(nFields-1)/2), GuiInfo.maxHeight);

if itemHeight==GuiInfo.maxHeight
    % Things are maxed out in size.
    % So we resize the window to make it look nice.
    newWindowHeight = (itemHeight*nFields + GuiInfo.usableVerticalRange(1) + (1-GuiInfo.usableVerticalRange(2)))*windowSize(2);

    itemHeight = itemHeight*windowSize(2)/newWindowHeight;

    buttonHeight = buttonHeight*windowSize(2)/newWindowHeight;

    MainFigureOptions.Children.Apply.position(4) = buttonHeight;
    MainFigureOptions.Children.Defaults.position(4) = buttonHeight;

    windowSize(2) = newWindowHeight;

    sizePads = round((ss(3:4)-windowSize));
    sizePads(1) = sizePads(1)/2; % We should use 2 right?
    sizePads(2) = sizePads(2)/2;
    pos = cat(2,sizePads,windowSize);

    MainFigureOptions.position = pos;
end

bottom = max(GuiInfo.usableVerticalRange(1),mean(GuiInfo.usableVerticalRange) - itemHeight*nFields/2 - (nFields-1)/2*itemHeight/2);

verticalStarts = fliplr(bottom + cumsum([0 0.5*ones(1,nFields-1)])*itemHeight + (0:(nFields-1))*itemHeight);

if length(verticalStarts) > 1
    editOffSet = [0 (verticalStarts(1)-verticalStarts(2)-itemHeight)/2 -0.01 0];
else
    editOffSet = [0 itemHeight/4 -0.01 0];
end

fontSize = 9;

for iField = 1:nFields

    %if mod(iField,2)
        %cBackgroundColor = max(MainPanelOptions.backgroundColor-0.05,[0 0 0]);
    %else
        cBackgroundColor = MainPanelOptions.backgroundColor;
    %end

    % Label string area
    MainFigureOptions.Children.MainPanel.Children.(['label' fields{iField}]).style = 'text';
    MainFigureOptions.Children.MainPanel.Children.(['label' fields{iField}]).units = 'normalized';
    MainFigureOptions.Children.MainPanel.Children.(['label' fields{iField}]).position = [GuiInfo.label.horizontal(1) verticalStarts(iField) GuiInfo.label.horizontal(2)-GuiInfo.label.horizontal(1) itemHeight];
    MainFigureOptions.Children.MainPanel.Children.(['label' fields{iField}]).fontweight = 'bold';
    MainFigureOptions.Children.MainPanel.Children.(['label' fields{iField}]).horizontalAlignment = 'right';
    MainFigureOptions.Children.MainPanel.Children.(['label' fields{iField}]).backgroundcolor = cBackgroundColor;
    MainFigureOptions.Children.MainPanel.Children.(['label' fields{iField}]).string = [guiPropNames{iField}  ':'];

    % Value string area
    MainFigureOptions.Children.MainPanel.Children.(['value' fields{iField}]).style = 'text';
    MainFigureOptions.Children.MainPanel.Children.(['value' fields{iField}]).units = 'normalized';
    MainFigureOptions.Children.MainPanel.Children.(['value' fields{iField}]).position = [GuiInfo.value.horizontal(1) verticalStarts(iField) GuiInfo.value.horizontal(2)-GuiInfo.value.horizontal(1) itemHeight];
    MainFigureOptions.Children.MainPanel.Children.(['value' fields{iField}]).horizontalAlignment = 'left';
    MainFigureOptions.Children.MainPanel.Children.(['value' fields{iField}]).backgroundcolor = cBackgroundColor;

    % Crerate the uicontrol object
    switch types{iField}
        case 'numeric'
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).style = 'edit';
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).units = 'normalized';
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).position = [GuiInfo.edit.horizontal(1) verticalStarts(iField) GuiInfo.edit.horizontal(2)-GuiInfo.edit.horizontal(1) itemHeight]+editOffSet;
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).callback = {@numericCallback fields{iField}};
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).fontSize = fontSize;

        case 'char'
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).style = 'edit';
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).units = 'normalized';
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).position = [GuiInfo.edit.horizontal(1) verticalStarts(iField) GuiInfo.edit.horizontal(2)-GuiInfo.edit.horizontal(1) itemHeight]+editOffSet;
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).callback = {@charCallback fields{iField}};
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).fontSize = fontSize;

        case 'function'
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).style = 'edit';
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).units = 'normalized';
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).position = [GuiInfo.edit.horizontal(1) verticalStarts(iField) GuiInfo.edit.horizontal(2)-GuiInfo.edit.horizontal(1) itemHeight]+editOffSet;
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).callback = {@functionCallback fields{iField}};
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).fontSize = fontSize;

        case 'boolean'
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).style = 'checkbox';
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).units = 'normalized';
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).position = [GuiInfo.edit.horizontal(1)+(GuiInfo.edit.horizontal(2)-GuiInfo.edit.horizontal(1))/2.25 verticalStarts(iField)+itemHeight/2 (GuiInfo.edit.horizontal(2)-GuiInfo.edit.horizontal(1))/4 itemHeight-itemHeight/3]; %[GuiInfo.edit.horizontal(1) verticalStarts(iField) GuiInfo.edit.horizontal(2)-GuiInfo.edit.horizontal(1) itemHeight];
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).callback = {@booleanCallback fields{iField}};
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).string = '';
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).backgroundColor = backgroundColor;
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).fontSize = fontSize;

        case 'structure'
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).style = 'pushbutton';
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).units = 'normalized';
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).position = [GuiInfo.edit.horizontal(1) verticalStarts(iField) GuiInfo.edit.horizontal(2)-GuiInfo.edit.horizontal(1) itemHeight]+editOffSet;
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).callback = {@structureCallback fields{iField}};
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).string = 'More Options';
            MainFigureOptions.Children.MainPanel.Children.(['edit' fields{iField}]).fontSize = fontSize;
    end

    MainFigureOptions.Children.MainPanel.Children.(['hline' num2str(iField)]).style = 'frame';
    MainFigureOptions.Children.MainPanel.Children.(['hline' num2str(iField)]).units = 'normalized';
    MainFigureOptions.Children.MainPanel.Children.(['hline' num2str(iField)]).position = [GuiInfo.label.horizontal(1) verticalStarts(iField) GuiInfo.edit.horizontal(2)-GuiInfo.label.horizontal(1) 1];
end

MainFigureOptions.Children.MainPanel.Children.('vline1').style = 'frame';
MainFigureOptions.Children.MainPanel.Children.('vline1').units = 'normalized';
MainFigureOptions.Children.MainPanel.Children.('vline1').position = [GuiInfo.label.horizontal(1) min(verticalStarts) 1 max(verticalStarts)-min(verticalStarts)+itemHeight*3/2];

MainFigureOptions.Children.MainPanel.Children.('vline2').style = 'frame';
MainFigureOptions.Children.MainPanel.Children.('vline2').units = 'normalized';
MainFigureOptions.Children.MainPanel.Children.('vline2').position = [mean([GuiInfo.label.horizontal(2) GuiInfo.value.horizontal(1)]) min(verticalStarts) 1 max(verticalStarts)-min(verticalStarts)+itemHeight*3/2];

MainFigureOptions.Children.MainPanel.Children.('vline3').style = 'frame';
MainFigureOptions.Children.MainPanel.Children.('vline3').units = 'normalized';
MainFigureOptions.Children.MainPanel.Children.('vline3').position = [mean([GuiInfo.value.horizontal(2) GuiInfo.edit.horizontal(1)]) min(verticalStarts) 1 max(verticalStarts)-min(verticalStarts)+itemHeight*3/2];

MainFigureOptions.Children.MainPanel.Children.('vline4').style = 'frame';
MainFigureOptions.Children.MainPanel.Children.('vline4').units = 'normalized';
MainFigureOptions.Children.MainPanel.Children.('vline4').position = [GuiInfo.edit.horizontal(2) min(verticalStarts) 1 max(verticalStarts)-min(verticalStarts)+itemHeight*3/2];

MainFigureOptions.Children.MainPanel.Children.(['hline' num2str(nFields+1)]).style = 'frame';
MainFigureOptions.Children.MainPanel.Children.(['hline' num2str(nFields+1)]).units = 'normalized';
MainFigureOptions.Children.MainPanel.Children.(['hline' num2str(nFields+1)]).position = [GuiInfo.label.horizontal(1) max(verticalStarts)+itemHeight*3/2 GuiInfo.edit.horizontal(2)-GuiInfo.label.horizontal(1) 1];

Handles = prtUtilSuicontrol(MainFigureOptions);

set(Handles.handle,'windowStyle','normal');

% Set the lines which are actually frames to have a width of 1 pixel
for iV = 1:4
    v = getpixelposition(Handles.MainPanel.(['vline' num2str(iV)]));
    setpixelposition(Handles.MainPanel.(['vline' num2str(iV)]),[v(1) v(2) 1 v(4)]);
end

for iField = 1:nFields+1
   v = getpixelposition(Handles.MainPanel.(['hline' num2str(iField)]));
   setpixelposition(Handles.MainPanel.(['hline' num2str(iField)]),[v(1) v(2) v(3) 1]);
end

modifiedProperties = false(nFields,1);
modifiedValues = cell(nFields,1);

updateValueFields;
updateEditFields;

uiwait(Handles.handle)

% Now that OutputStructure is set we need to set the values
modifiedValues = modifiedValues(modifiedProperties);
modifiedNames = guiPropNames(modifiedProperties);

output = obj;
for iProp = 1:length(modifiedNames)
    output.(modifiedNames{iProp}) = OutputStructure.(modifiedNames{iProp});
end

if isstruct(obj)
    creationString = '';
else
    creationString = sprintf('%s(',class(obj));
    for iProp = 1:length(modifiedNames)
        creationString = cat(2,creationString,'''',modifiedNames{iProp},''',',modifiedValues{iProp},',');
    end
    if strcmpi(creationString(end),',');
        creationString = creationString(1:(end-1));
    end
    creationString = cat(2,creationString,')');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function numericCallback(myHandles,eventData,fieldName) %#ok

        shouldRevert = false;
        newNumber = [];
        try
            eval(sprintf('newNumber = %s;',get(myHandles,'string')))
        catch ME %#ok<NASGU>
            errordlg('This is not valid MATLAB syntax for a numeric value.','Invalid numeric value.','modal')
            shouldRevert = true;
        end

        try
            objTest = obj;
            objTest.(fieldName) = newNumber;
        catch ME
            errordlg(ME.message,'modal')
            shouldRevert = true;
        end

        if shouldRevert
            updateEditFields;
        else

            propInd = find(strcmp(guiPropNames,fieldName),1,'first');
            modifiedProperties(propInd) = true;
            modifiedValues{propInd} = get(myHandles,'string');

            OutputStructure.(fieldName) = newNumber;
            updateValueFields;
            updateEditFields;
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function charCallback(myHandles,eventData,fieldName) %#ok

        shouldRevert = false;
        newString = [];
        try
            eval(sprintf('newString = ''%s'';',get(myHandles,'string')))
        catch ME %#ok<NASGU>
            errordlg('This is not valid MATLAB syntax for a string value.','Invalid string','modal')
            shouldRevert = true;
        end

        try
            objTest = obj;
            objTest.(fieldName) = newString;
        catch ME
            errordlg(ME.message,'modal')
            shouldRevert = true;
        end

        if shouldRevert
           updateEditFields;
        else
            propInd = find(strcmp(guiPropNames,fieldName),1,'first');
            modifiedProperties(propInd) = true;
            modifiedValues{propInd} = get(myHandles,'string');

            OutputStructure.(fieldName) = newString;
            updateValueFields;
            updateEditFields;
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function functionCallback(myHandles,eventData,fieldName) %#ok

        shouldRevert = false;

        cString = get(myHandles,'string');
        if length(cString) < 1
            errordlg('This is not valid MATLAB syntax for a function.','Invalid function','modal');
            shouldRevert = true;
        end

        if ~shouldRevert
            try
                eval(sprintf('newFunctionOutput = %s;',cString));

                try
                    objTest = obj;
                    objTest.(fieldName) = newFunctionOutput;
                catch ME
                    errordlg(ME.message,'modal')
                    shouldRevert = true;
                end

            catch ME %#ok<NASGU>
                errordlg('This is not valid MATLAB syntax for a function.','Invalid function','modal')
                shouldRevert = true;
            end
        end

        if shouldRevert
            updateEditFields;
        else
            propInd = find(strcmp(guiPropNames,fieldName),1,'first');
            modifiedProperties(propInd) = true;
            modifiedValues{propInd} = cString;

            OutputStructure.(fieldName) = eval(cString);
            updateValueFields;
            updateEditFields;
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     function structureCallback(myHandles,eventData,fieldName) %#ok
%
%         %newMainFigureOptions = rmfield(newMainFigureOptions,'CloseRequestFcn'); % This is really important. We will need the handle to the newly create cancelCallback not this one.
%         %newMainPanelOptions = MainPanelOptions;
%
%         [OutputStructure.(fieldName), subCreationString] = prtUtilObjectGuiSimple(OutputStructure.(fieldName));
%
%         updateValueFields;
%         updateEditFields;
%
%     end
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function booleanCallback(myHandles,eventData,fieldName) %#ok

        OutputStructure.(fieldName) = get(myHandles,'value');

        propInd = find(strcmp(guiPropNames,fieldName),1,'first');
        modifiedProperties(propInd) = true;
        if OutputStructure.(fieldName)
            modifiedValues{propInd} = 'true';
        else
            modifiedValues{propInd} = 'false';
        end

        updateValueFields;
        updateEditFields;
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function defaultCallback(myHandles,eventdata) %#ok

        OutputStructure = DefaultStructure;

        updateValueFields;
        updateEditFields;

    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function cancelCallback(myHandles,eventdata) %#ok

        set(Handles.handle,'visible','off');

        OutputStructure = BaseStructure;

        delete(Handles.handle);
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function applyCallback(myHandles,eventdata) %#ok
        delete(Handles.handle);
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function updateValueFields
        for jField = 1:nFields
            switch types{jField}
                case 'numeric'
                    set(Handles.MainPanel.(['value' fields{jField}]),'string',mat2str(OutputStructure.(fields{jField})));
                case 'char'
                    set(Handles.MainPanel.(['value' fields{jField}]),'string',OutputStructure.(fields{jField}));
                case 'function'
                    set(Handles.MainPanel.(['value' fields{jField}]),'string',func2str(OutputStructure.(fields{jField})));
                case 'boolean'
                    trueFalse = {'False','True'};
                    set(Handles.MainPanel.(['value' fields{jField}]),'string',trueFalse{OutputStructure.(fields{jField})+1});
                case 'structure'
                    set(Handles.MainPanel.(['value' fields{jField}]),'string','To view and edit additional parameters use the button to the right.');

            end
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function updateEditFields
        for jField = 1:nFields
            switch types{jField}
                case 'numeric'
                    set(Handles.MainPanel.(['edit' fields{jField}]),'string',mat2str(OutputStructure.(fields{jField})));
                case 'char'
                    set(Handles.MainPanel.(['edit' fields{jField}]),'string',OutputStructure.(fields{jField}));
                case 'function'
                   set(Handles.MainPanel.(['edit' fields{jField}]),'string',func2str(OutputStructure.(fields{jField})));
                case 'boolean'
                    set(Handles.MainPanel.(['edit' fields{jField}]),'value',OutputStructure.(fields{jField}));
                case 'structure'
                    % We shouldn't have to do anything.
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% mcs.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [xbest,fbest,xmin,fmi,ncall,ncloc,flag]=mcs(fcn,data,u,v,prt,
% smax,nf,stop,iinit,local,gamma,hess)
% MCS global optimization for function defined by fcn in the
% n-dimensional box [u,v]
%
% Input:
% fcn = 'fun' 	name of function fun(data,x), x an n-vector
% data		data vector (or other data structure)
% [u,v]       	box in which the optimization is carried out (u, v
%             	n-vectors)
% prt		print level
% 		prt = 0: no printing
% 		prt = 1: # sweep, minimal nonempty level, # f-calls,
% 		best point and function value (default)
% 		prt > 1: only meaningful for test functions with known
% 		global minimizers
% 		in addition levels and function values of boxes
% 		containing the global minimizers of a test function
% smax        	number of levels (default: 5*n+10)
% nf         	maximum number of function evaluations (default: 50*n^2)
% stop         	stop(1) in ]0,1[:  relative error with which the known
%		 global minimum of a test function should be found
%		 stop(2) = fglob known global minimum of a test function
%		 stop(3) = safeguard parameter for absolutely small
%		 fglob
%		stop(1) >= 1: the program stops if the best function
%		 value has not been improved for stop(1) sweeps
%		stop(1) = 0: the user can specify a function value that
%		 should be reached
%                stop(2) = function value that is to be achieved
%             	(default: stop = 3*n)
% iinit       	parameter defining the initialization list
%             	= 0        corners and midpoint (default for finite u,v)
%             	= 1        safeguarded version *default otherwise)
% 		= 2        5u/6 + v/6, u/6 + 5v/6 and midpoint
%		= 3        initialization list with line searches
%             	otherwise  self-defined init. list (to be stored in
%			   init0.m)
%		for a self-defined initialization list, the user should
% 		provide an m-script file init0.m containing a matrix x0
%		with n rows and at least 3 columns and two n-vectors l
%		and L
%		the ith column of x0 contains the initialization list
%		values for the ith coordinate, their number is L(i), and
%		x0(i,l(i)) is the ith coordinate of the initial point
% local		local = 0: no local search
%		otherwise: maximal number of steps in local search
%		(default: 50)
% gamma		stopping criterion for local search (default: eps)
%           	the local search is stopped if abs(g)'*max(abs(x),
%		abs(xold)) < gamma*(f0-f)
% hess		sparsity pattern of the Hessian for local search
%		(default: hessian = ones(n,n))
%
% Output:
% xbest(1:n)  	current best point
% fbest    	function value at xbest
% xmin        	matrix with n rows; the columns are the points in the
%             	'shopping basket' (i.e. good points resp. local
%		minimizers)
% fmi         	function values corresponding to the 'shopping basket';
%		fmi(i) is the function value at xmin(:,i)
% ncall       	number of function evaluations
% ncloc		number of function evaluations used for local search
% flag        	specifies which stopping criterion has been used
%             	= 0  a (known) global minimum fglob of a test function
%                    has been found with the required relative error
%		     relerr
%             	= 1  the division procedure has been completed
%             	= 2  the maximum number nf of function calls has been
%                    reached without finding a known minimum with the
%                    required relative error or completing the division
%                    procedure
%		= 3  stop(1) sweeps without progress (for stop(1) >= 1)
%
% Uses the following m-files (directly or indirectly):
% addloc.m
% basket.m
% basket1.m
% bounds.m
% chkloc.m
% chrelerr.m
% chvtr.m
% csearch.m 	called by lsearch.m
% exgain.m
% fbestloc.m
% genbox.m
% hessian.m	called by triple.m
% init.m
% initbox.m
% initlist.m
% gls.m and its subprograms   called by lsearch.m
% lsearch.m
% minq.m and its subprograms  called by lsearch.m
% neighbor.m
% polint.m	called by exgain.m and initbox.m
% polint1.m	called by triple.m
% quadmin.m
% quadpol.m
% range.m	called by lsearch.m
% splinit.m
% split.m
% split1.m
% split2.m   	called by splrnk.m
% splrnk.m
% strtsw.m
% subint.m
% triple.m	called by lsearch.m
% updtf.m    	called by vertex.m
% updtoptl.m
% updtrec.m   	called by splinit.m and split.m
% vert1.m    	called by vertex.m
% vert2.m    	called by vertex.m
% vert3.m    	called by vertex.m
% vertex.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xbest,fbest,xmin,fmi,ncall,ncloc,flag]=mcs(fcn,data,u,v,prt,smax,nf,stop,iinit,local,gamma,hess)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% global variables
global foptbox nbasket nboxes ncall nglob nsweep nsweepbest optlevel  record xglob xloc
% foptbox(1:nglob)  function value(s) of the box(es) containing the (a)
%             	global minimizer of a test function
% nbasket   	counter for boxes in the 'shopping basket'
% nboxes      	counter for boxes not in the 'shopping basket'
% nglob       	number of global minimizers of a test function
% nloc		(for local ~= 0) counter of points that have been used
% 		as starting points for a local search
% nsweep      	sweep counter
% nsweepbest    number of sweep in which fbest was updated for the last
%		time
% optlevel    	level(s) of the box(es) containing the (a) global
%             	minimum of a test function
% record(1:smax-1) record(i) points to the best non-split box at level i
%             	(record list)
% xglob(1:n,1:nglob)  xglob(:,i), i=1:nglob, are the global minimizers
% of a test function in [u,v]
% xloc(1:n,:)	(for local ~= 0) columns are the points that have been
%		used as starting points for local search
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = length(u);

% check box bounds
if ~isempty(find(v<u))
  error('incompatible box bounds')
elseif ~isempty(find(u==v))
  error('degenerate box bound')
end

% default values for the input parameters
if nargin < 5, prt = 1; end
if nargin < 6, smax = 5*n+10; end
if nargin < 7, nf = 50*n^2; end
if nargin < 8, stop = 3*n; end
if nargin < 9,
  if isempty(find(isinf(u))) & isempty(find(isinf(v)))
    iinit = 0;
  else
    iinit = 1;
  end
end
if nargin < 10, local = 50; end
if nargin < 11, gamma = eps; end
if nargin < 12, hess = ones(n,n); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xmin=[];fmi=[];		% avoid warnings in Matlab5

% initial values for the numbers of function calls (total number/local
% search)
ncall = 0;
ncloc = 0;

% some parameters needed for initializing large arrays
step1 = 10000;
step = 1000;
dim = step1;

% initialization of some large arrays
isplit = zeros(1,step1);
level = zeros(1,step1);
ipar = zeros(1,step1);
ichild = zeros(1,step1);
f = zeros(2,step1);
z = zeros(2,step1);
nogain = zeros(1,step1);

% definition of the initialization list
if iinit == 0
  x0(:,1) = u;
  x0(:,2) = (u+v)/2;
  x0(:,3) = v;
  l = 2*ones(n,1);
  L = 3*ones(n,1);
elseif iinit == 1
  for i = 1:n
    if u(i) >= 0
      x0(i,1) = u(i); [x0(i,2),x0(i,3)] = subint(u(i),v(i));x0(i,2) = 0.5*(x0(i,1)+x0(i,3));
    elseif v(i) <= 0
      x0(i,3) = v(i); [x0(i,2),x0(i,1)] = subint(v(i),u(i));x0(i,2) = 0.5*(x0(i,1)+x0(i,3));
    else
      x0(i,2) = 0; [xi,x0(i,1)] = subint(0,u(i)); [xi,x0(i,3)] = subint(0,v(i));
    end
  end
  l = 2*ones(n,1);
  L = 3*ones(n,1);
elseif iinit == 2
  x0(:,1) = (5*u + v)/6;
  x0(:,2) = 0.5*(u + v);
  x0(:,3) = (u + 5*v)/6;
  l = 2*ones(n,1);
  L = 3*ones(n,1);
elseif iinit == 3
  [x0,f0,l,L,istar,ncall1] = initlist(fcn,data,u,v);
  ncall = ncall + ncall1;
else
  init0 	%self-defined initialization list
  for i=1:size(x0,2)
    if ~isempty(find(x0(:,i)<u)) | ~isempty(find(x0(:,i)>v))
      error('incorrect initialization list')
    end
  end
end

% check whether there are infinities in the initialization list
if ~isempty(find(isinf(x0))), error('infinities in ititialization list'), end

% computation of the function values f0 appertaining to the init. list
% and the pointer istar to the best point in the initialization list
if iinit ~= 3
  [f0,istar,ncall1] = init(fcn,data,x0,l,L,n);
  ncall = ncall + ncall1;
end

% definition of the base vertex of the original box
for i = 1:n
  x(i) = x0(i,l(i));
end

% definition of the opposite vertex v1 of the original box
for i = 1:n
  if abs(x(i)-u(i)) > abs(x(i)-v(i))
    v1(i) = u(i);
  else
    v1(i) = v(i);
  end
end

% initialization of the record list, the counters nboxes, nbasket, m
% and nloc, xloc and the output flag
record = zeros(smax-1,1);
nboxes = 1;
nbasket = 0;
nbasket0 = 0;
nsweep = 0;
m = n;
record(1) = 1;
nloc = 0;
xloc = [];
flag = 1;

[ipar,level,ichild,f,isplit,p,xbest,fbest] = initbox(x0,f0,l,L,istar,u,v,prt);
% generates the boxes in the initialization procedure
f0min = fbest;
if stop(1) > 0 & stop(1) < 1
  flag = chrelerr(fbest,stop);
elseif stop(1) == 0
  flag = chvtr(fbest,stop(2));
end
if ~flag,return,end
% if the (known) minimum function value fglob has been found with the
% required tolerance, flag is set to 0 and the program is terminated

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s = strtsw(smax,level,f(1,:));
% the vector record is updated, and the minimal level s containing
% non-split boxes is computed
nsweep = nsweep + 1;	% sweep counter

while s < smax & ncall + 1 <= nf
  par = record(s);   % the best box at level s is the current box
  [n0,x,y,x1,x2,f1,f2] = vertex(par,n,u,v,v1,x0,f0,ipar,isplit,ichild,z,f,l,L);
  % compute the base vertex x, the opposite vertex y, the 'neighboring'
  % vertices and their function values needed for quadratic
  % interpolation and the vector n0 indicating that the ith coordinate
  % has been split n0(i) times in the history of the box
  if s > 2*n*(min(n0)+1)
  % s 'large'
    [isplit(par),z(2,par)] = splrnk(n,n0,p,x,y);
    % splitting index and splitting value z(2,par) for splitting by
    % rank are computed
    % z(2,par) is set to Inf if we split according to the init. list
    splt = 1;  % indicates that the box is to be split
  else
    if nogain(par) % box has already been marked as not eligible for splitting
                   % by expected gain
      splt = 0;
    else
      [e,isplit(par),z(2,par)] = exgain(n,n0,l,L,x,y,x1,x2,f(1,par),f0,f1,f2);
      % splitting by expected gain
      % compute the expected gain vector e and the potential splitting
      % index and splitting value
      fexp = f(1,par) + min(e);
      if fexp < fbest
        splt = 1;
      else
        splt = 0;  % the box is not split since we expect no improvement
        nogain(par) = 1; % the box is marked as not eligible for splitting by expected gain
      end
    end
  end
  if splt == 1  % prepare for splitting
    i = isplit(par);
    level(par) = 0;
    if z(2,par) == Inf % prepare for splitting by initialization list
      m = m + 1;
      z(2,par) = m;
      [xbest,fbest,f0(:,m),xmin,fmi,ipar,level,ichild,f,flag,ncall1] = splinit(fcn,data,i,s,smax,par,x0,n0,u,v,x,y,x1,x2,L,l,xmin,fmi,ipar,level,ichild,f,xbest,fbest,stop,prt);
      ncall = ncall + ncall1;
    else  % prepare for default splitting
      z(1,par) = x(i);
      [xbest,fbest,xmin,fmi,ipar,level,ichild,f,flag,ncall1] = split(fcn,data,i,s,smax,par,n0,u,v,x,y,x1,x2,z(:,par),xmin,fmi,ipar,level,ichild,f,xbest,fbest,stop,prt);
      ncall = ncall + ncall1;
    end
    if nboxes > dim
% if the pre-assigned size of the `large' arrays has already been exceeded, these arrays are made larger
      isplit(nboxes+1:nboxes+step) = zeros(1,step);
      level(nboxes+1:nboxes+step) = zeros(1,step);
      ipar(nboxes+1:nboxes+step) = zeros(1,step);
      ichild(nboxes+1:nboxes+step) = zeros(1,step);
      z(:,nboxes+1:nboxes+step) = zeros(2,step);
      nogain(nboxes+1:nboxes+step) = zeros(1,step);
      f(:,nboxes+1:nboxes+step) = zeros(2,step);
      dim = nboxes + step;
    end
    if ~flag,break,end
  else  % splt=0: no splitting, increase the level by 1
    if s + 1 < smax
      level(par) = s + 1;
      updtrec(par,s+1,f(1,:));
    else
      level(par) = 0;
      nbasket = nbasket + 1;
      xmin(:,nbasket) = x;
      fmi(nbasket) = f(1,par);
    end
    if prt > 1
      [w1,w2] = bounds(n,n0,x,y,u,v);
      % compute lower and upper bounds of the box in order to be able
      % to check whether it contains a global minimizer
      iopt = [];
      % the vector iopt contains the indices of the global minimizers
      % contained in the box
      for iglob = 1:nglob
        if w1 <= xglob(:,iglob) & xglob(:,iglob) <= w2
          iopt = [iopt, iglob];
        end
        for iglob = 1:length(iopt)
          optlevel(iopt(iglob)) = s + 1;
        end
      end
    end
  end % of prepare for splitting
  s = s + 1;
  while s < smax
    if record(s) == 0
      s = s + 1;
    else
      break
    end
  end
  if s == smax  % if smax is reached, a new sweep is started
    if local,
      [fmi(nbasket0+1:nbasket),j] = sort(fmi(nbasket0+1:nbasket));
      xmin(:,nbasket0+1:nbasket) = xmin(:,nbasket0+j);
      xmin0 = [];
      fmi0 = [];
      for j = nbasket0+1:nbasket
        x = xmin(:,j);
        f1 = fmi(j);
        chkloc;
        if loc,
          addloc;
          [xbest,fbest,xmin,fmi,x,f1,loc,flag,ncall1] = basket(fcn,data,x,f1,xmin,fmi,xbest,fbest,stop,nbasket0);
          ncall = ncall + ncall1;
          if ~flag,break,end
          if loc,
            [xmin1,fmi1,nc,flag] = lsearch(fcn,data,x,f1,f0min,u,v,nf-ncall,stop,local,gamma,hess);
            ncall = ncall + nc;
            ncloc = ncloc + nc;
            if fmi1 < fbest
              xbest = xmin1;
              fbest = fmi1;
              nsweepbest = nsweep;
              if ~flag
                nbasket0 = nbasket0 + 1;
                nbasket = nbasket0;
                xmin(:,nbasket) = xmin1;
                fmi(nbasket) = fmi1;
                break
              end
              if stop(1) > 0 & stop(1) < 1
                flag = chrelerr(fbest,stop);
              elseif stop(1) == 0
                flag = chvtr(fbest,stop(2));
              end
              if ~flag,return,end
            end
            [xbest,fbest,xmin,fmi,loc,flag,ncall1] = basket1(fcn,data,xmin1,fmi1,xmin,fmi,xbest,fbest,stop,nbasket0);
            ncall = ncall + ncall1;
            if ~flag,break,end
            if loc,
              nbasket0 = nbasket0 + 1;
              xmin(:,nbasket0) = xmin1;
              fmi(nbasket0) = fmi1;
              fbestloc;
              if ~flag,
                nbasket = nbasket0; break
              end
            end
          end
        end
      end
      nbasket = nbasket0;
      if ~flag,break,end
    end
    s = strtsw(smax,level,f(1,:));
    if prt,
      if nsweep == 1
        fprintf('nsw  minl  ');
        if prt > 1
          fprintf('optl    fopt       ')
        end
        fprintf('nf     fbest        xbest\n')
      end
      minlevel=s;
      fprintf('%3i  %3i',nsweep,minlevel);
      if prt > 1
        fprintf('  %3i',optlevel);fprintf('  %10.3e',foptbox);
      end
      fprintf('  %5i  %10.3e',ncall,fbest);
      fprintf('  %10.4f',xbest);
      fprintf(1,'\n');
    end
    if stop(1) > 1
      if nsweep - nsweepbest >= stop(1),flag = 3; return,end
    end
    nsweep = nsweep + 1;
  end
end
if ncall >= nf
  flag = 2;
end
if local,
  if length(fmi) > nbasket
    xmin(:,nbasket+1:length(fmi)) = [];
    fmi(nbasket+1:length(fmi)) = [];
  end
end

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function [ a, b ] = p44_lim ( )

%*****************************************************************************80
%
%% P44_LIM returns the integration limits for problem 44.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end