text:
function AAM = build_model_2d_from_files(directory, varargin)
% AAM = build_model_2d_from_files(directory)
% AAM = build_model_2d_from_files(directory, 'triangulation', triangulation)
%
% FIXME: Mostly a cut/paste from build_model_2d, but this
% method was needed to avoid out-of-memory issues when dealing with large training
% sets.
%
% Builds an Active Appearance Model using shape and appearance data retrieved from
% the given directory.
% Based on the 'Active Appearance Models Revisited' paper by Iain Matthews and Simon Baker.
%
%                         PARAMETERS
%
% directory Path to a directory containing shape and appearance data in the appropriate format.
%           See the examples to see what this format is, or use one of the provided datasets.
%
%                     OPTIONAL PROPERTIES
%
% 'triangulation' A matrix property providing a triangular mesh over the
%                 shape (in the same format as the one returned by delaunay).
%                 If not provided, a triangulation will be determined using
%                 delaunay on the mean shape.
%
%
%                        RETURN VALUE
%
% AAM is a structure with the following fields:
%
% 	s0 The mean shape as a column vector [i1 i2 i3 ... j1 j2 j3 ...]'.
%
%   s  The matrix of base shapes with each shape stored on a column vector
%      with the same structure as s0
%
%   s_star The basis for the global normalizing transform. They are four column
%          vectors whose linear combination allow for translation, scale and
%          rotation of the base shape
%
%   shape_eiv The eigenvalues associated to each shape principal
%             component
%
%   warp_map A modelh-by-modelw lookup table for determining to which
%            triangle warp_map(i,j) the pixel (i,j) belongs. If 0, the
%            pixel is outside the mean shape
%
%   shape_mesh The triangles indexed by warp_map. It is a Tx3 matrix
%              containing indices to vertices in the shape. Look at
%              the delaunay function for further information
%
%   alpha_coords An image containing the alpha barycentric coordinate
%                for pixel (i,j) wrt of the triangle given by
%                shape_mesh(warp_map(i,j),:). Value is undefined if
%                warp_map(i,j) is 0
%
%   beta_coords An image containing the beta barycentric coordinate
%               for pixel (i,j) wrt of the triangle given by
%               shape_mesh(warp_map(i,j),:). Value is undefined if
%               warp_map(i,j) is 0
%
%   adj_list Array of N cells containing indices to the triangles sharing
%            a given vertex. Such that for the k-th cell, the list will
%            contain the indices to the triangles in shape_mesh sharing
%            the k-th vertex.
%
%   A0 Mean appearance as a column vector, sorted firstly by color
%      component, then j coordinate and finally i coordinate
%      Example:
%       X1 = [ 1 2; 3 4];
%       X2 = [ 5 6; 7 8];
%       X3 = [ 9 10; 11 12];
%       X(:,:,1)=X1;
%       X(:,:,2)=X2;
%       X(:,:,3)=X3;
%       reshape(X,[],1) =>
%          _________
%        1  j==1
%        3 _____   C1
%        2  j==2
%        4 _________
%        5  j==1
%        7 _____   C2
%        6  j==2
%        8 _________
%        9  j==1
%       11 _____   C3
%       10  j==2
%       12 _________
%
%   A  The matrix of base appearances, with each appearance stored on
%      the colomuns in the same manner as A0
%
%   app_eiv The eigenvalues associated to each appearance principal
%           component
%
%   dA0 The gradient of the mean appearance stored such that for pixel (i,j),
%       dA0(i,j,c,1) gives the i direction of the gradient for color c
%       and dA0(i,j,c,2) gives the j direction of the gradient for color c
%
%   dW_dp Jacobian of the warp with respect to the shape parameters p
%
%   dN_dq Jacobian of the global normalizing warp with respect of the
%         global transform parameters q
%
%   SD Steepest descent images. They are stored as row vectors (of the same length
%      as A0). The first 4 images are associated with the q
%      paramenters, the other are associated with the p parameters
%
%   H The hessian, it's a square NxN matrix where N is the number of steepest
%     descent images (or, equivalently, the sum of the number of p and
%     q parameters)
%
%   invH Inverse hessian, it's a square NxN matrix where N is the number of steepest
%     descent images (or, equivalently, the sum of the number of p and
%     q parameters)
%
%   R Simply the product between the inverse hessian and the steepest descent images,
%     invH * SD. Encodes the linear relationship between the error image and the incremental
%     parameters.
%
% Author: Luca Vezzaro (elvezzaro@gmail.com)

	shape_triangles = zeros(0, 3, 'uint32');

	i = 1;
	while i < numel(varargin)
		if strcmp(varargin{i}, 'triangulation')
			% Be sure the triangle indices are unsigned integers.
			shape_triangles = uint32(varargin{i+1});
		else
			error(sprintf('Unsupported property: %s.', varargin{i}));
		end

		i = i + 2;
	end

	tic;

	landmark_files = dir(sprintf('%s/*.mat', directory));
	disp 'Number of shapes:'
	% size of training set
	ns = numel(landmark_files)

	for j=1:ns
		% Load associated landmarks
		load(sprintf('%s/%s', directory, landmark_files(j).name));
		info = imfinfo(sprintf('%s/%s', directory, landmark_files(j).name(1:end-4)));
		shape_data(:,:,j) = xy2ij(annotations, info.Height);
	end

	% number of points
	np = size(shape_data, 1);

	% use first shape as initial mean
	mean_shape = shape_data(:,:,1);

	reference_shape = mean_shape;

	% Matrix containing the aligned shapes
	aligned_data = shape_data;

	% Less iterations should be enough, but better safe than sorry, efficiency
	% is not an issue here
	for it=1:100
		for i=1:ns
			% Align each shape to the mean
			[ D Y ] = procrustes(mean_shape, aligned_data(:,:,i));
			aligned_data(:,:,i) = Y;
		end

		% New mean shape, store in a different variable so we can compare it to
		% the old one
		new_mean_shape = mean(aligned_data, 3);
		[ D Y ] = procrustes(reference_shape, new_mean_shape);
		new_mean_shape = Y;

		mean_shape = new_mean_shape;
	end

	mean_shape = mean(aligned_data, 3);

	% Determine the region of interest
	mini = min(mean_shape(:,1));
	minj = min(mean_shape(:,2));
	maxi = max(mean_shape(:,1));
	maxj = max(mean_shape(:,2));

	% Place the origin in the upper left corner of the rectangle
	% representing the bounding box of the mean shape. Add a 1
	% pixel offset to avoid artifacts in gradient computations.
	mean_shape = mean_shape - repmat([mini - 2, minj - 2], [np, 1]);

	disp 'Width of model:'
	% Determine width and height of the model, add 1 pixel offset
	% to avoid gradient artifacts
	modelw = ceil(maxj - minj + 3)
	disp 'Height of model:'
	modelh = ceil(maxi - mini + 3)

	disp 'Time required for shape data alignment:'

	toc

	%%
	%% DEBUG STUFF
	%%
%	for i=1:ns
%		plot(aligned_data(:,2,i), aligned_data(:,1,i), 'o');
%		hold on;
%		pause;
%	end
%	plot(mean_shape(:,2), mean_shape(:,1), 'rx');
%	hold off;
%	return;

	%%
	%% END OF DEBUG STUFF
	%%

	tic

	% Definitive mean shape for AAM, stored on a column vector [i1 i2 i3 ... j1 j2 j3 ...]'
	AAM.s0 = reshape(mean_shape, 2*np, 1);

	% Prepare shape data for PCA
	shape_matrix = reshape(aligned_data, [2*np ns]) - repmat(AAM.s0, [1 ns]);
	clear aligned_data;

	% Get modes covering 98% of variation
	[pc eiv] = var_pca(shape_matrix, 0.98);
	clear shape_matrix;

	disp 'Number of Shape Modes:'
	size(pc, 2)

	AAM.shape_eiv = eiv;

	% Build the basis for the global shape transform, we do it here because
	% they are used to orthonormalize the shape principal vectors
	% It is done differently to the paper as we're using a different coordinate
	% frame. Here u -> i, v -> j
	s1_star = AAM.s0;
	s2_star(1:np) = -AAM.s0(np+1:end);
	s2_star(np+1:2*np) = AAM.s0(1:np);
	s3_star(1:np) = ones(np,1);
	s3_star(np+1:2*np) = zeros(np,1);
	s4_star(1:np) = zeros(np,1);
	s4_star(np+1:2*np) = ones(np,1);

	% Stack the basis we found before with the shape basis so
	% we can orthonormalize
	s_star_pc(:,1) = s1_star;
	s_star_pc(:,2) = s2_star;
	s_star_pc(:,3) = s3_star;
	s_star_pc(:,4) = s4_star;

	s_star_pc(:,5:size(pc,2)+4) = pc;

	% Orthogonalize the basis (should already be close to orthogonal)
	s_star_pc = gs_orthonorm(s_star_pc);

	% Basis for the global shape transform
	AAM.s_star = s_star_pc(:,1:4);
	% Basis for the shape model
	AAM.s = s_star_pc(:,5:end);

	disp 'Time required for building shape model:'
	toc

	%%
	%% DEBUG STUFF
	%%

%	disp 'Plotting principal vectors as deformation directions of the mean...'
%
%	plot_s0 = reshape(AAM.s0, [np 2]);
%
%	for i=1:size(AAM.s, 2)
%		hold off;
%		plot(plot_s0(:,2), plot_s0(:,1), 'o');
%		hold on;
%
%		fprintf('Principal vector %d\n.', i);
%
%		delta_s = reshape(AAM.s(:,i), [np 2]);
%		plot_s = plot_s0 + delta_s * 16;
%
%		for j=1:5
%			plot(plot_s(:,2), plot_s(:,1), 'r.');
%			plot_s = plot_s + delta_s * 16;
%		end
%		hold off;
%		disp 'Waiting for keypress...'
%		pause;
%	end
%	return;

	%%
	%% END OF DEBUG STUFF
	%%

	tic

	% Check if a triangulation was already provided
	if size(shape_triangles, 1) == 0
		% Triangle mesh used for warping. We determine it on the base shape.
		% Returned result is double, convert to unsigned integer.
		shape_triangles = delaunay(AAM.s0(np+1:2*np)', AAM.s0(1:np)');
		%triplot(shape_triangles, AAM.s0(1:np)', AAM.s0(np+1:2*np)');
	end

	% Bitmap mapping pixels to the triangle they belong to (the triangle used
	% for warping)
	warp_map = zeros(modelh, modelw, 'uint32');

	% Bitmap containing the first barycentric coordinate of each pixel wrt to the
	% triangle specified by warp_map
	alpha_coords = zeros(modelh, modelw);

	% Bitmap containing the second barycentric coordinate of each pixel wrt to the
	% triangle specified by warp_map
	beta_coords = zeros(modelh, modelw);

	% For each pixel, find the triangle it belongs to
	for j=1:modelw
		for i=1:modelh
			% For each triangle
			for k=1:size(shape_triangles, 1)
				t = shape_triangles(k,:);

				% Vertices of the triangle in the mean shape
				i1 = AAM.s0(t(1));
				j1 = AAM.s0(np + t(1));
				i2 = AAM.s0(t(2));
				j2 = AAM.s0(np + t(2));
				i3 = AAM.s0(t(3));
				j3 = AAM.s0(np + t(3));

				% Compute the two barjcentric coordinates
				den = (i2 - i1) * (j3 - j1) - (j2 - j1) * (i3 - i1);
				alpha = ((i - i1) * (j3 - j1) - (j - j1) * (i3 - i1)) / den;
				beta = ((j - j1) * (i2 - i1) - (i - i1) * (j2 - j1)) / den;

				if alpha >= 0 && beta >= 0 && (alpha + beta) <= 1
					% Found the triangle, save data to the bitmaps and break
					warp_map(i,j) = k;
					alpha_coords(i,j) = alpha;
					beta_coords(i,j) = beta;
					break;
				end
			end
		end
	end

	% Build adjacency list
	AAM.adj_list = cell(np, 1);

	for i=1:size(shape_triangles, 1)
		for j=1:3
			v = shape_triangles(i,j);
			AAM.adj_list{v} = [AAM.adj_list{v} i];
		end
	end

	% IMPORTANT: If the datatype is changed from 'uint32', the mex-file will
	% need to be changed!
	AAM.warp_map = uint32(warp_map);
	AAM.shape_mesh = uint32(shape_triangles);
	AAM.alpha_coords = alpha_coords;
	AAM.beta_coords = beta_coords;

	disp 'Time required for warp pre-computations:'
	toc

	%%
	%% DEBUG STUFF
	%%
%	color_map = rand(max(max(warp_map)), 3);
%	figure(1)
%	imshow(warp_map, color_map);
%
%	figure(2)
%	imshow(alpha_coords);
%
%	figure(3)
%	imshow(beta_coords);
%
%	return;

	%%
	%% END OF DEBUG STUFF
	%%

	tic

	% FIXME: This is somewhat of an hack, but the code is so slow when matrices grow inside loops...
	app = double(imread(sprintf('%s/%s', directory, landmark_files(1).name(1:end-4)))) / 255;
	app_matrix = zeros(modelw*modelh*size(app,3), ns);

	disp 'Color components:'
	nc = size(app, 3)

	for i=1:ns
		% Prepare appearance data
		app = double(imread(sprintf('%s/%s', directory, landmark_files(i).name(1:end-4)))) / 255;
		%app = convn(convn(app, exp(-(-5:5).^2)./sum(exp(-(-5:5).^2)), 'same'), exp(-(-5:5).^2)./sum(exp(-(-5:5).^2))', 'same');

		app_matrix(:,i) = reshape(pa_warp(AAM, shape_data(:,:,i), app), [(modelw*modelh*nc) 1]);
	end

	disp 'Time required for warping all input appearances:'
	toc

	clear app;
	clear landmark_files;

	%%
	%% DEBUG STUFF
	%%

%disp 'Displaying warped images...'
%
%for i=1:ns
%	fprintf('Image %d\n.', i);
%
%	imshow(reshape(app_matrix(:,i), [modelh modelw nc]));
%
%	hold on;
%
%	triplot(shape_triangles, AAM.s0(np+1:2*np)', AAM.s0(1:np)');
%	hold off;
%
%	disp 'Waiting for keypress...'
%	pause;
%end
%hold off;
%return;

	%%
	%% END OF DEBUG STUFF
	%%

	tic

	% FIXME: Should we normalize pixel values like done by Stegmann to
	% improve the generalization of the model?

	% Mean appearance as a column vector, sorted firstly by color
	% component, then j coordinate and finally i coordinate
	%Example:
	%       X1 = [ 1 2; 3 4];
	%       X2 = [ 5 6; 7 8];
	%       X3 = [ 9 10; 11 12];
	%       X(:,:,1)=X1;
	%       X(:,:,2)=X2;
	%       X(:,:,3)=X3;
	%       reshape(X,[],1) =>
	%          _________
	%        1  j==1
	%        3 _____   C1
	%        2  j==2
	%        4 _________
	%        5  j==1
	%        7 _____   C2
	%        6  j==2
	%        8 _________
	%        9  j==1
	%       11 _____   C3
	%       10  j==2
	%       12 _________

	AAM.A0 = mean(app_matrix, 2);
	mean_app = reshape(AAM.A0, [modelh modelw nc]);

	%app_size = modelw * modelh * nc;

	% Prepare apperances for PCA
	app_matrix = app_matrix - repmat(AAM.A0, [1 ns]);

	% FIXME: Goes out of memory without 'econ' flag, look at "An Introduction to Active
	% Shape Models" by Cootes for an optimized PCA
	%[pc eiv] = var_pca(app_matrix, 0.95);
  [pc eiv] = app_pca(app_matrix, 0.95);

	clear app_matrix;

	disp 'Number of Appearance Modes:'
	size(pc, 2)

	AAM.A = pc;
	AAM.app_eiv = eiv;

	disp 'Time required for building apperance model:'
	toc

	%%
	%% DEBUG STUFF
	%%

%  disp 'Displaying results of PCA on appearances...'
%
%	for i=1:size(AAM.A,2)
%   cur_app = reshape(AAM.A(:,i), [modelh modelw nc]);
%
%		fprintf('Showing mean appearance and sums of mean and principal component %d with increasing weight\n.', i);
%
%		subplot(1,3,1), imshow(mean_app)
%		subplot(1,3,2), imshow(mean_app + cur_app*sqrt(eiv(i))*3)
%		subplot(1,3,3), imshow(mean_app - cur_app*sqrt(eiv(i))*3)
%		%subplot(1,4,4), imshow(mean_app + cur_app*sqrt(eiv(i))*3)
%
%		disp 'Waiting for keypress...'
%		pause;
%	end
%	hold off;
%	return;
	%%
	%% END OF DEBUG STUFF
	%%

	tic

	% Now we get into the interesting (and tricky) stuff

	for i=1:nc
		% For each color, estimate the gradient of the mean appearance.
		% dcol is the estimated derivative wrt to columns (which for us is the
		% j coordinate), drow is wrt to rows (which for us is the i coordinate)
		[di dj] = gradient_2d(mean_app(:,:,i), AAM.warp_map);
		mean_app_gradient(:,:,i,1) = di;
		mean_app_gradient(:,:,i,2) = dj;
	end

	AAM.dA0 = mean_app_gradient;

	%%
	%% DEBUG STUFF
	%%
%	figure(1)
%imshow(mean_app);
%	mn = min(min(min(min(AAM.dA0))));
%	mx = max(max(max(max(AAM.dA0))));
%
%	plot_dA0 = (AAM.dA0 - repmat(mn, [modelh modelw nc 2])) / (mx - mn);
%	figure(2)
%	imshow(plot_dA0(:,:,:,1));
%	figure(3);
%	imshow(plot_dA0(:,:,:,2));
%	return;
	%%
	%% END OF DEBUG STUFF
	%%

	disp 'Time required for calculating image gradient:'
	toc

	tic

	dW_dp = zeros(modelh, modelw, 2, size(AAM.s,2));
	dN_dq = zeros(modelh, modelw, 2, 4);

	for j=1:modelw
		for i=1:modelh
			if warp_map(i,j) ~= 0
				% Only the vertices of the triangle containing the pixel are of relevance
				% in determining the Jacobian.
				t = shape_triangles(warp_map(i,j),:);

				% FIXME: this is not how it should be done, there is a way to juggle
				% with the barycentric coordinates we already computed to do the same
				for k=1:3
					% Derivative of vertex 'k' wrt to the shape parameters
					dik_dp = AAM.s(t(k),:);
					djk_dp = AAM.s(t(k)+np,:);

					% Derivative of vertex 'k' wrt to the global transformation parameters
					dik_dq = AAM.s_star(t(k),:);
					djk_dq = AAM.s_star(t(k)+np,:);

					% Now we need the barycentric coordinates of (i,j) computed using
					% point 'k' as the origin. So we rearrange the order of the vertices
					% (it doesn't matter if the result has clockwise or counterclockwise
					% winding)
					t2 = t;
					t2(1) = t(k);
					t2(k) = t(1);

						% Vertices of the triangle in the mean shape
					i1 = AAM.s0(t2(1));
					j1 = AAM.s0(np + t2(1));
					i2 = AAM.s0(t2(2));
					j2 = AAM.s0(np + t2(2));
					i3 = AAM.s0(t2(3));
					j3 = AAM.s0(np + t2(3));

					% Compute the two barycentric coordinates
					den = (i2 - i1) * (j3 - j1) - (j2 - j1) * (i3 - i1);
					alpha = ((i - i1) * (j3 - j1) - (j - j1) * (i3 - i1)) / den;
					beta = ((j - j1) * (i2 - i1) - (i - i1) * (j2 - j1)) / den;

					% Nonzero portion of the Jacobian of the warp wrt to the vertex 'k'.
					% TODO: Isn't this the third barycentric coordinate? It should be
					% easier and faster to compute it directly.
					dW_dij = 1 - alpha - beta;

					% This is how it's formulated on the paper (dW_dq is similar):
					% dW_dp = dW_dp + [dW_dij; 0] * dik_dp + [0; dW_dij] * djk_dp;
					% MATLAB does weird things with dimensions here, so we need squeeze()
					% to make sure the submatrix is actually a 2xn matrix and not a
					% 1x1x2xn
					dW_dp(i,j,:,:) = squeeze(dW_dp(i,j,:,:)) + dW_dij * [dik_dp; djk_dp];
					dN_dq(i,j,:,:) = squeeze(dN_dq(i,j,:,:)) + dW_dij * [dik_dq; djk_dq];
				end
			end
		end
	end

	disp 'Time required for calculating warp Jacobians:'
	toc

	%%
	%% DEBUG STUFF
	%%
%  disp 'Displaying the warp Jacobian...'
%
%	for i=1:size(dW_dp,4)
%	  jacobian_i = dW_dp(:,:,1,i);
%	  jacobian_j = dW_dp(:,:,2,i);
%
%
%		fprintf('Showing warp jacobian i and j components for parameter %d\n.', i);
%
%		subplot(1,2,1), imshow(jacobian_i)
%		subplot(1,2,2), imshow(jacobian_j)
%
%		disp 'Waiting for keypress...'
%		pause;
%	end
%	return;
	%%
	%% END OF DEBUG STUFF
	%%

	tic

	app_modes = reshape(AAM.A, [modelh modelw nc size(AAM.A, 2)]);

	SD = zeros(modelh, modelw, nc, 4 + size(dW_dp, 4));

	% Compute steepest descent images for the 4 global transformation parameters
	% TODO: the loops over c are probably optimizable, but the number of dimensions involved
	% are starting to hurt my head...
	% TODO: test me
	for i=1:4
		prj_diff = zeros(nc, size(AAM.A, 2));
		for j=1:size(AAM.A, 2)
			for c=1:nc
				prj_diff(c,j) = sum(sum(app_modes(:,:,c,j) .* (AAM.dA0(:,:,c,1) .* dN_dq(:,:,1,i) + AAM.dA0(:,:,c,2) .* dN_dq(:,:,2,i))));
			end
		end

		for c=1:nc
			SD(:,:,c,i) = AAM.dA0(:,:,c,1) .* dN_dq(:,:,1,i) + AAM.dA0(:,:,c,2) .* dN_dq(:,:,2,i);
		end

		for j=1:size(AAM.A, 2)
			for c=1:nc
				SD(:,:,c,i) = SD(:,:,c,i) - prj_diff(c,j) * app_modes(:,:,c,j);
			end
		end
	end

	% Compute steepest descent images for the shape parameters
	% TODO: the loops over c are probably optimizable, but the number of dimensions involved
	% are starting to hurt my head...
	% TODO: test me
	for i=1:size(dW_dp, 4)
		prj_diff = zeros(nc, size(AAM.A, 2));
		for j=1:size(AAM.A, 2)
			for c=1:nc
				prj_diff(c,j) = sum(sum(app_modes(:,:,c,j) .* (AAM.dA0(:,:,c,1) .* dW_dp(:,:,1,i) + AAM.dA0(:,:,c,2) .* dW_dp(:,:,2,i))));
			end
		end

		for c=1:nc
			SD(:,:,c,i+4) = AAM.dA0(:,:,c,1) .* dW_dp(:,:,1,i) + AAM.dA0(:,:,c,2) .* dW_dp(:,:,2,i);
		end

		for j=1:size(AAM.A, 2)
			for c=1:nc
				SD(:,:,c,i+4) = SD(:,:,c,i+4) - prj_diff(c,j) * app_modes(:,:,c,j);
			end
		end
	end

	AAM.SD = zeros(size(SD, 4), size(AAM.A,1));

	% FIXME: A posteriori optimization fix, integrate it more nicely with the code.
	for i=1:size(SD, 4)
		AAM.SD(i,:) = reshape(SD(:,:,:,i), 1, []);
	end

	%%
	%% DEBUG STUFF
	%%
%	mn = min(min(min(min(SD))));
%	mx = max(max(max(max(SD))));
%	SD = (SD - repmat(mn, [modelh modelw nc size(SD, 4)])) / (mx - mn);
%	for i=1:size(SD,4)
%		h=figure(i+1);
%		imshow(SD(:,:,:,i));
%	end
	%%
	%% END OF DEBUG STUFF
	%%

	clear SD;

	disp 'Time required for calculating steepest descent images:'
	toc

	tic

	% Compute the Hessian, still not sure if this is the right
	% way to treat the color information (using greyscale images
	% give pretty much the same matching results so it should be ok)
	AAM.H = AAM.SD * AAM.SD';

	AAM.invH = inv(AAM.H);

	AAM.R = AAM.invH * AAM.SD;

	disp 'Time required for determining and inverting Jacobian and calculating invH*SD:'
	toc

author: Sable


text:
% IM = mkZonePlate(SIZE, AMPL, PHASE)
%
% Make a "zone plate" image:
%     AMPL * cos( r^2 + PHASE)
% SIZE specifies the matrix size, as for zeros().
% AMPL (default = 1) and PHASE (default = 0) are optional.

% Eero Simoncelli, 6/96.

function [res] = mkZonePlate(sz, ampl, ph)

sz = sz(:);
if (size(sz,1) == 1)
  sz = [sz,sz];
end

mxsz = max(sz(1),sz(2));

%------------------------------------------------------------
%% OPTIONAL ARGS:

if (exist('ampl') ~= 1)
  ampl = 1;
end

if (exist('ph') ~= 1)
  ph = 0;
end

%------------------------------------------------------------

res = ampl * cos( (pi/mxsz) * mkR(sz,2) + ph );

repo: SelfExSR


text:
function scanList = chooseScans(view)
% scanList = schooseScans(view);
%
%   Query the user for a sub-selection of scans.

%   This function uses numScans(view) to determine the number of scans to choose from.
%   An alternate function, selectScans, actually checks for the existence of the tSeries files.
%   Use selectScans if you will be analyzing the tSeries.
%   Use chooseScans if your code does not depend on the presence/absence of the tSeries files.
%
% Output:
%  scanList: list of selected scans.
%
% djh, 3/2001, modified from selectScans

nScans = numScans(view);
scanList = 1:nScans;

for i=scanList
    scanNames{i} = ['Scan',num2str(i)];
end

%Check for zero:
if nScans == 0
  myErrorDlg('No scans!');
  return
end

% Which scans to analyze?
iSel = buttondlg('Choose scans', scanNames);
scanList = scanList(find(iSel));

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
add_paths;
data_paths;

%class = 'chair';
%class = 'sofa';
%class = 'bus';
class = 'aeroplane';
%class = 'bicycle';
%class = 'car';
%class = 'motorbike';
%class = 'diningtable';

dir_name = fullfile(pascal_result_dir, class);
dir_pascal = fullfile(PASCAL3D_dir, 'Annotations', [class,'_pascal']);
dir_save_img = fullfile(pascal_images_dir, class);
if ~exist(dir_save_img)
    mkdir(dir_save_img);
end
filenames = dir(dir_name);

maskFac = 1;
images_in = zeros(numel(filenames)-2, 256*192);
fig = figure('Visible', 'off', 'Position', [10000,10000,400,300]); % offscreen figure
set(fig,'color','white')

numCores = 8;
parpool(numCores);
for i = 3:numel(filenames)
    fprintf("%d / %d \n", i-2, numel(filenames));
    file_result = filenames(i).name;
    result = load(fullfile(dir_name, file_result), 'result');
    result = result.result;
    I = result.I.image;

    imshow(I);
    hold on;
    Hmask = imshow(ones(size(I, 1), size(I, 2)));
    Hmask.AlphaData = result.I.segmentation*maskFac;
    I_frame = frame2im(getframe(fig)); % convert plot to image (true color RGB matrix)

    I_resized = imresize(I_frame, [192, 256], 'bicubic'); % resize image
    imgfile = fullfile(dir_save_img, [file_result(1:11), '.png']);
    imwrite(I_resized, imgfile); % save image to file

    I_gray = rgb2gray(I_resized);
    imshow(I_gray)
    images_in(i-2,:) = I_gray(:);
end

h5create('images_pascal.h5','/IN',[numel(filenames)-2 256*192])
h5write('images_pascal.h5', '/IN', images_in)

fprintf("Done! \n")

author: jhonykaesemodel


text:
function [fx,J,dfdp,d2fdxdp] = f_fullDCM4fmri(Xt,Theta,ut,inF)
% DCM for fMRI evolution function (attn: embedding of HRF!)
% function [fx,dF_dX,dF_dTheta] = f_fullDCM4fmri(Xt,Theta,ut,inF)
% This function evaluates the evolution funciton of the neuronal level in
% DCM for fMRI. For the sake of HRF convolution purposes, it also
% internally calls g_fullDCM4fmri.m so that the hemodynamic states are
% updated properly.

%- Neuronal states evolution
[fx,J,dfdp] = f_dcm4fmri(Xt,Theta,ut,inF);
d2fdxdp = [];

%- HRF convolution
g_fullDCM4fmri(Xt,[],[],inF.inG);
% This call is used to update the hemodynamic states that are convolved
% outputs of the neuronal states (Xt), as well as the necessary gradient,
% i.e. Jacobian and derivative w.r.t. hemodynamic parameters.

author: MBB-team

repo: VBA-toolbox


text:
function [ know, x ] = p31_sol ( n )

%*****************************************************************************80
%
%% P31_SOL returns the solution for problem 31.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the problem.  This value
%    is only needed for those problems with variable N.
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Output, real X(N), the solution, if known.
%
  know = 1;

  x = 4.0 * ones ( n, 1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function mrk= mrk_sortChronologically(mrk, varargin)
%MRK_SORTCHRONOLOGICALLY - Sort Markers Chronologically
%
%Synopsis:
%  MRK_OUT= mrk_sortChronologically(MRK_IN, <OPT>)
%
%Arguments:
%  MRK_IN: Marker structure as received by file_loadBV
%  OPT: struct or property/value list of optional properties:
%    'RemoveVoidClasses': Void classes are removed from the list of classes.
%    'Classwise': Each class is sorted chronologically, default 0.
%
%Returns:
% MRK_OUT: Marker structure with events sorted chronologically

% Benjamin Blankertz

props= {'Sort',              0  'BOOL';
        'RemoveVoidClasses'  1  'BOOL';
        'Classwise'          0  '!BOOL';
        };
props_selectEvents= mrk_selectEvents;

if nargin==0,
  mrk= opt_catProps(props, props_selectEvents);
  return
end

opt= opt_proplistToStruct(varargin{:});
opt= opt_setDefaults(opt, props);
opt_checkProplist(opt, props, props_selectEvents);

misc_checkType(mrk, 'STRUCT(time)');

if opt.Classwise,
  nClasses= size(mrk.y,1);
  si= zeros(nClasses, length(mrk.time));
  for ci= 1:nClasses,
    idx= find(mrk.y(ci,:));
    [so,sidx]= sort(mrk.time(idx));
    si(ci, 1:length(idx))= idx(sidx);
  end
  si= si(find(si));
else
  [so,si]= sort(mrk.time);
end

opt_selectEvents= opt_substruct(opt, props_selectEvents(:,1));
mrk= mrk_selectEvents(mrk, si, opt_selectEvents, 'Sort',0);

author: bbci


text:
classdef DASCMOP9 < PROBLEM
% <multi> <real> <large/none> <constrained>
% Difficulty-adjustable and scalable constrained benchmark MOP

%------------------------------- Reference --------------------------------
% Z. Fan, W. Li, X. Cai, H. Li, C. Wei, Q. Zhang, K. Deb, and E. Goodman,
% Difficulty adjustable and scalable constrained multi-objective test
% problem toolkit, Evolutionary Computation, 2020, 28(3): 339-378.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Wenji Li

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 3;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            x_all = X(:,3:1:end);
            temp  = size(x_all,2);
            sum1  = sum((x_all - cos(0.25 * pi * temp * (X(:,1) + X(:,2)) / obj.D)) .^2,2);
            PopObj(:,1) = cos(0.5 * pi * X(:,1)) .* cos(0.5 * pi * X(:,2)) + sum1;
            PopObj(:,2) = cos(0.5 * pi * X(:,1)) .* sin(0.5 * pi * X(:,2)) + sum1;
            PopObj(:,3) = sin(0.5 * pi * X(:,1)) + sum1;
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,X)
            x_all  = X(:,3:1:end);
            temp   = size(x_all,2);
            sum1   = sum((x_all - cos(0.25 * pi * temp * (X(:,1) + X(:,2)) / obj.D)) .^2,2);
            PopCon = Constraint(X(:,1),X(:,2),sum1);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,3);
            R = R./repmat(sqrt(sum(R.^2,2)),1,3);
            X(:,2) = atan(R(:,2)./R(:,1))/0.5/pi;
            X(:,1) = acos(R(:,1)./cos(0.5*pi*X(:,2)))/0.5/pi;
            C(:,1) = -sin(20*pi*X(:,1));
            C(:,2) = -cos(20*pi*X(:,2));
            R(any(C>1e-2,2),:) = [];
            R = R + 0.5;
        end
        %% Generate the feasible region
        function R = GetPF(obj)
            [X1,X2] = meshgrid(linspace(0,1,100));
            x = cos(0.5*pi*X1).*cos(0.5*pi*X2);
            y = cos(0.5*pi*X1).*sin(0.5*pi*X2);
            z = sin(0.5*pi*X1);
            fes = all(Constraint(X1(:),X2(:),0.5)<=0,2);
            z(reshape(~fes,size(z))) = nan;
            R = {x+0.5,y+0.5,z+0.5};
        end
    end
end

function PopCon = Constraint(X1,X2,sum1)
    % set the parameters of constraints
    DifficultyFactors = [0.5,0.5,0.5];
    % Type-I parameters
    a = 20;
    b = 2 * DifficultyFactors(1) - 1;
    % Type-II parameters
    d = 0.5;
    if DifficultyFactors(2) == 0.0
        d = 0.0;
    end
    e = d - log(DifficultyFactors(2));
    if isfinite(e) == 0
        e = 1e+30;
    end
    % Type-III parameters
    r = 0.5 * DifficultyFactors(3);
    % Calculate objective values
    PopObj(:,1) = cos(0.5 * pi * X1) .* cos(0.5 * pi * X2) + sum1;
    PopObj(:,2) = cos(0.5 * pi * X1) .* sin(0.5 * pi * X2) + sum1;
    PopObj(:,3) = sin(0.5 * pi * X1) + sum1;
    % Type-I constraints
    PopCon(:,1) = b - sin(a * pi * X1);
    PopCon(:,2) = b - cos(a * pi * X2);
    % Type-II constraints
    PopCon(:,3) = -(e - sum1) .* (sum1 - d);
    if DifficultyFactors(2) == 1.0
        PopCon(:,3) = 1e-4 - abs(sum1 - e);
    end
    % Type-III constraints
    x_k = [1.0, 0.0, 0.0, 1.0 / sqrt(3.0)];
    y_k = [0.0, 1.0, 0.0, 1.0 / sqrt(3.0)];
    z_k = [0.0, 0.0, 1.0, 1.0 / sqrt(3.0)];
    for k=1:length(x_k)
        PopCon(:,3+k) = r * r - ((PopObj(:,1) - x_k(k))).^2  -...
            ((PopObj(:,2) - y_k(k))).^2 - ((PopObj(:,3) - z_k(k))).^2;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function plotMultipleAmps_SingleCondition(view)
%
% plotMultipleAmps_SingleCondition(view)
%
% Bar plot of the amplitudes for each scan, averaging across
% all pixels (in all slices) in a selection of ROIs. All y-axes are made the same. The bar heights
% and a coarse SEM can be obtained from get(gca,'UserData').
%
% gmb  5/25/98
% bw   2/19/99  Added seY field to the UserData field.
%	    seY is an estimate of the variability in the
%      amplitudes.  It is the SEM of the in the complex
%      (amp*exp(-i*ph)) representation.  The values are
%      computed in vectorMean.m
% fwc   11/07/02 plots data relative to current view
%       added plotting of multiple ROIs
%       ROI selection copied from plotMultipleTSeries.m
% Plots the mean amplitudes for each ROI in the current scan
% Based on plotMultipleAmps ARW 072803

mrGlobals;

% Select ROIs
nROIs=size(view.ROIs,2);
roiList=cell(1,nROIs);
for r=1:nROIs
    roiList{r}=view.ROIs(r).name;
end
selectedROIs = find(buttondlg('ROIs to Plot',roiList));
nROIs=length(selectedROIs);
if (nROIs==0)
    error('No ROIs selected');
end

% Plot it
selectGraphWin
clf
fontSize = 8;
headerStr = ['Mean Amplitudes'];
set(gcf,'Name',headerStr);

minylim=0;
maxylim=0;
nrows=0;
ncols=0;

nscans = numScans(view);
ROIamps=zeros(nscans,nROIs);
ROIseZ=zeros(nscans,nROIs);
ROImeanPhs=zeros(nscans,nROIs);

for r=1:nROIs

    n=selectedROIs(r);
    view = selectROI(view,n); % is there another way?
    [meanAmps,meanPhs,seZ] = vectorMeans(view);
    ROIamps(:,r)=meanAmps(:);
    ROIseZ(:,r)=seZ(:);
    ROImeanPhs(:,r)=meanPhs(:);
    %xstr{r}=[view.ROIs(selectedROIs(r)).name];

    roiName{r}=view.ROIs(selectedROIs(r)).name;
    fprintf(['\nROI #%d :',roiName{r}],r);
    xstr{r}=roiName{r};
end

% Now do the plotting

% Only plotting the current scan
   r=getCurScan(view);

   subplot(1,1,1);

    %plot the bar graph
    h=mybar(ROIamps(r,:)',ROIseZ(r,:)',xstr);
    xlabel('ROI','FontSize',fontSize);
    ylabel('Mean Amplitude','FontSize',fontSize);
    set(gca,'FontSize',ceil(fontSize*1.2));
    conditionName{r}=dataTYPES(view.curDataType).scanParams(r).annotation;
    fprintf(['\nCondition #%d :',conditionName{r}],r);
    title(conditionName{r});

%Save the data in gca('UserData')
data.y =ROIamps(r,:);

data.seY = ROIseZ(r,:); % this should probably be adapted

set(gca,'UserData',data);

return;


text:
function lam = prod_lambda_msgs(n, cs, msg, msg_type, except)

if nargin < 5, except = -1; end

lam = msg{n}.lambda_from_self;
switch msg_type
  case 'd',
   for i=1:length(cs)
     c = cs(i);
     if c ~= except
       lam = lam .* msg{n}.lambda_from_child{i};
     end
   end
 case 'g',
  if isinf(lam.precision) % isfield(lam, 'observed_val')
    return; % pass on the observed msg
  end
   for i=1:length(cs)
     c = cs(i);
     if c ~= except
       m = msg{n}.lambda_from_child{i};
       lam.precision = lam.precision + m.precision;
       lam.info_state = lam.info_state + m.info_state;
     end
   end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function y=PSNR_RGB(X,Y)

% Y= PSNR_RGB(X,Y)
% Computes the Peak Signal to Noise Ratio for two RGB images
% Class input : double [0,1] ,
% july ,25, 2012
% KHMOU Youssef

if size(X)~=size(Y)
    error('The images must have the same size');
end

%if ~isa(X,'double')
%   X=double(X)./255.00;
%end
%if  ~isa(Y,'double')
%    Y=double(Y)./255.00;
%end

% begin
d1=max(X(:));
d2=max(Y(:));
d=max(d2);
sigma=mean2((X-Y).^2);

y=10*log10((d.^2)./sigma);

author: andrewssobral


text:
function [inference] = tapas_linear_prepare_inference(data, model, inference)
%%
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

nc = numel(model.graph{2}.htheta.T);
np = size(data.u, 1);
nb = size(data.u{1}, 2);

inferece.kernel = cell(4, 1);

inference.kernel{2} = struct('k', [eye(nb)], 's', [0.5]);
inference.kernel{3} = struct('k', [eye(nb)], 's', [0.5]);

end

author: translationalneuromodeling

repo: tapas


text:
% Frequency response 3rd Order Bessel HPF
%  and G array.
% File:  c:\M_files\short_updates\UsingGarray2.m
% 02/10/07
% See Word file UsingGarray2.doc for schemtic and equations.
clear;clc;
% unit suffixes
u=1e-6;K=1e3;m=1e-3;u=1e-6;n=1e-9;p=1e-12;
% component values
R1=52*K;R2=75*K;R3=287*K;
C1=1.55*n;C2=C1;C3=C1;
Ein=1; % Unity input for (normalized) transfer function
N=3; % Number of capacitors = order of circuit
%
% Form W, Q, S, and P arrays:
%
W=[R2 -R2 0;0 R1 R3-R1;0 0 R3];
Q=-[1 0 0;1 1 0;1 1 1];S=Ein*[1;1;1];P=diag([C1 C2 C3]);
%
% Get A, B, D, & E arrays:
%
C=inv(W*P);A=C*Q;B=C*S;
% D & E
D=[0 0 0];E=0;
%
F=[0 0 R3]; % R3 is the coefficient of both iC1 and iC2 in Vo2
% which are derivative terms (iC=dVc/dt, etc.)
G=F*P;
%
% * * * * * * * * * * * * Frequency response * * * * * * * * * * * *
%
I=eye(N); % identity matrix
%
% Log frequency sweep from BF to BF+ND
%
BF=2;ND=2;PD=50;NP=ND*PD+1;Fr=logspace(BF,BF+ND,NP);
for i=1:NP
   s=2*pi*Fr(i)*j; % j = sqrt(-1)
   stm=(s*I-A)\B;
   v2=abs((D+G*A)*stm+E+G*B);Vo(i)=20*log10(v2);
   asym(i)=60*log10(Fr(i)/1000); % +60 dB/decade slope
end
%
h=plot(log10(Fr),Vo,'r',log10(Fr),asym,'k--');
set(h,'LineWidth',2);
grid on;
axis([2 4 -60 20]);
ylabel('dBV');title('Bessel HPF');
xlabel('Log Freq(Hz)');
legend('Vo','Asymptote');
figure(1);
%

repo: mcbench-benchmarks


text:
%% test modulate for all oasis functions.
col = {[0 114 178],[0 158 115], [213 94 0],[230 159 0],...
    [86 180 233], [204 121 167], [64 224 208], [240 228 66]}; % colors
plot_cvx = false;

%% example 1:  foopsi, AR1 model. This model is used when the sampling rate is low
g = 0.95;         % AR coefficient
noise = .3;
T = 3000;
framerate = 30;
firerate = 0.5;
b = 0;              % baseline
N = 1;              % number of trials
seed = 13;          % seed for genrating random variables
[y, true_c, true_s] = gen_data(g, noise, T, framerate, firerate, b, N, seed);

% case 1: all parameters are known
lambda = 2.4;
params.p = 1;
params.B = 1;
[c_oasis, s_oasis] = deconvolveCa(y, 'mcmc', params);  %#ok<*ASGLU>
[c_cvx, s_cvx] = foopsi(y, g, lambda);

figure('name', 'FOOPSI, AR1, known: g, lambda', 'papersize', [15, 4]);
plot_cvx = true;
show_results;
plot_cvx = false;

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function name = RemoveExt(filename)
%
%
%        name = RemoveExt(filename)
%
%
%        Description: remove file extension
%
%        Input:
%           -filename: the name of the file
%
%        Output:
%           -name: the filename without extension
%
%     Copyright (C) 2013  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

k = strfind(filename, '.');

if(size(k) == 0) %no extension
    name = filename;
else
    k = k(end); %get the real extension
    name = filename(1:(k - 1));
end

end

author: banterle


text:
function [slcstack] = SLC_filt(mlistack,slcstack,SHP,Coh_cal,BroNumthre,Cohthre_slc_filt)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   This file is part of TomoSAR.
%
%   TomoSAR is distributed in the hope that it will be useful,
%   but without warranty of any kind; without even the implied
%   warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%   See the Apache License for more details.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Author : Dinh Ho Tong Minh (INRAE) and Yen Nhi Ngo, Jan. 2022
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[~,~,n_slc] = size(slcstack);

mask_coh = Coh_cal > Cohthre_slc_filt;
mask_PS = SHP.BroNum>BroNumthre;
mask = and(mask_PS,mask_coh);
mask = repmat(mask,[1,1,n_slc]);

slcstack(mask) = abs(mlistack(mask)).*exp(1i*angle(slcstack(mask)));

author: DinhHoTongMinh

repo: TomoSAR


text:
function outPref = techPref(inPref)
%TECHPREF   Preference settings for TRIGTECH.
%   P = TRIGTECH.TECHPREF() returns a structure P with fields which contain the
%   default TRIGTECH preferences as field/value pairs.  This structure may be
%   passed to the TRIGTECH constructor.
%
%   P = TRIGTECH.TECHPREF(Q) does the same but replaces the default TRIGTECH
%   preferences with the values specified by the field/value pairs in the input
%   structure Q.
%
%   Note that no checking of either the input preference names or values takes
%   place and that preference names are case-sensitive.
%
%   ABSTRACT PREFERENCES REQUIRED OF ALL TECHS
%
%     chebfuneps   - Relative tolerance used in construction and subsequent
%      [2^-52]       operations.  See TRIGTECH.HAPPINESSCHECK for more details.
%
%     maxLength    - Maximum number of points used by the constructor.
%      [2^16+1]
%
%     fixedLength  - Fixed number of points used by constructor.  NaN allows
%      [NaN]         adaptive construction.
%
%     extrapolate
%       true       - Function values at endpoints may be extrapolated from
%                    interior values rather than sampled.
%      [false]     - Do not extrapolate values at endpoints.
%
%     sampleTest
%      [true]      - Tests the function at one more arbitrary point to
%                    minimize the risk of missing signals between grid
%                    points.
%       false      - Do not test.
%
%   TRIGTECH-SPECIFIC PREFERENCES
%
%     gridType     - Type of equi-spaced grid used to sample the function.
%         1        - Equally spaced grid starting at -1+h/2
%        [2]       - Equally spaced grid starting at -1
%
%     minSamples    - Minimum number of points used by the constructor.  Should
%      [17]          be of the form 2^n + 1.  If not, it is rounded as such.
%
%     refinementFunction - Define function for refining sample values.
%       'nested'         - Nest the evaluations to reduce the number of function evaluations.
%      ['resampling']    - Every function value is computed afresh as the
%                          constructor tries grids of size 9, 17, 33, etc.
%       function_handle  - A user-defined refinement function.  See REFINE.m
%
%     happinessCheck     - Define function for testing happiness.
%      ['standard']      - Standard check routine
%       'classic'        - Use the default process from Chebfun v4.
%       'strict'         - Strict tolerance for coefficients.
%       'loose'          - A looser tolerance for coefficients.
%       function_handle  - A user defined happiness. See HAPPINESSCHECK.m
%
% See also TRIGTECH, CHEBTECH, CHEBTECH1, CHEBTECH2

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

outPref.chebfuneps         = 2^-52;
outPref.gridType           = 2;
outPref.minSamples         = 17;
outPref.maxLength          = 2^16;
outPref.fixedLength        = NaN;
outPref.extrapolate        = false;
outPref.sampleTest         = true;
outPref.refinementFunction = 'nested';
outPref.happinessCheck     = 'standard';

if ( nargin == 1 )
    validPrefs = fieldnames(outPref);
    for ( givenPref = fieldnames(inPref).');
        givenPref = givenPref{1};
        if ( ~any(strcmp(givenPref, validPrefs)) )
            warning('CHEBFUN:TRIGTECH:techPref:unknownPref', ...
                ['Unrecognized input preference ''' givenPref '''.']);
        end
    end

    outPref = chebfunpref.mergeTechPrefs(outPref, inPref);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function minimize( varargin )

%MINIMIZE Specifiies a convex (or affine) objective to be maximized.

cvx_problem = evalin( 'caller', 'cvx_problem', '[]' );
if isempty( cvx_problem ) || ~isa( cvx_problem, 'cvxprob' ),
    error( 'A cvx problem does not exist in this scope.' );
end
if ~isempty( cvx_problem.direction ),
    if isequal( cvx_problem.direction, 'find' ),
        error( 'Objective functions cannot be added to sets.' );
    else
        error( 'An objective function has already been supplied.' );
    end
end

if nargin < 1,
    error( 'Objective expression missing.' );
elseif iscellstr( varargin ),
    arg = evalin( 'caller', sprintf( '%s ', varargin{:} ) );
elseif nargin > 1,
    error( 'Too many input arguments.' );
else
    arg = varargin{1};
end

if ~isa( arg, 'cvx' ) && ~isa( arg, 'double' ) && ~isa( arg, 'sparse' ),
    error( 'Cannot accept an objective of type ''%s''.', class( arg ) );
end
persistent remap
if isempty( remap ),
    remap = cvx_remap( 'convex', 'log-convex' );
end
vx = remap( cvx_classify( arg ) );
if ~all( vx ),
    error( 'Disciplined convex programming error:\n   Cannot minimize a(n) %s expression.', cvx_class(arg(vx==0),false,true) );
end

newobj( cvx_problem, 'minimize', arg );

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

repo: REST


text:
%
% The original code was downloaed from
% Ruslan Salakhutdinov\s hompeage:
%   http://http://www.mit.edu/~rsalakhu/
%
function ls = logdiff(xx,dim)
% ls = logsum(x,dim)
%
% returns the log of diff of logs
% similar to logsum.m function

if(length(xx(:))==1) ls=xx; return; end

xdims=size(xx);
if(nargin<2)
  dim=find(xdims>1);
end

alpha = max(xx,[],dim)-log(realmax)/2;
repdims=ones(size(xdims)); repdims(dim)=xdims(dim);
ls = alpha+log(diff(exp(xx-repmat(alpha,repdims)),dim));

author: kyunghyuncho

repo: deepmat


text:
% loadeep() - Load EEProbe continuous data (*.cnt).
%
% Usage:
%   >> [eep] = loadeep(file, varargin);
%
% Inputs:
%            filename    - name of ANT EEP file, including extension (*.cnt)
%
%   Optional inputs:
%           'time1'      - start time in seconds.
%           'sample1'    - start at sample1, default 0.
%                          (Overrules time1)
%           'time2'      - end time of data to be  read in seconds, default = whole file.
%           'sample2'    - end sample of data to be read, default = whole file.
%                          (Overrules time2)
% Outputs:
%   [eep]                - data structure holding continous EEG information and other
%                          relevant information.
%
% Author: Maarten-Jan Hoeve, ANT Software, Enschede, The Netherlands, 8 October 2003
%
% See also: eeglab(), pop_loadeep()
%

%123456789012345678901234567890123456789012345678901234567890123456789012

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

% $Log: loadeep.m,v $
% Revision 1.3  2005/07/21 10:12:51  mvelde
% updated information
%
% Revision 1.2  2005/06/08 08:16:37  mvelde
% converted files to unix format
%
% Revision 1.1  2004/11/26 13:22:22  jwiskerke
% Added eeglab to cvs.
%
% Revision 1.4  2003/10/24 13:34:40  Maarten-Jan Hoeve
% Added GNU Licence and updated revision history
%
% Revision 1.1.2  2003/10/17 09:55:20  mvelde
% updated: consistent copyrights, arguments/data labels, fixed some typos
%
% Advanced Neuro Technology (ANT) BV, The Netherlands, www.ant-neuro.com / info@ant-neuro.com
%

function r=loadeep(file, varargin)

if nargin < 1
	help loadeep;
	return;
end;

% defaults
if ~any(file=='.'), file=[file '.cnt']; end
[datdir,name,ext]=fileparts(file);
disp(['Loading file ' file ' ...'])

% read short piece of data to get sampling rate, channels etc
eeg = read_eep_cnt(file,1,2);

if ~isempty(varargin)
	 r=struct(varargin{:});
else r = [];
end;

% add defaults
try, r.time1;       catch, r.time1=0; end
try, r.sample1;     catch, r.sample1=[]; end
try, r.time2;       catch, r.time2=[]; end
try, r.sample2;     catch, r.sample2=[]; end

r.filename=file;
r.totsamples=eeg.nsample;
r.nchannels=eeg.nchan;
r.rate=eeg.rate;

% Create struct for holding channel labels
for i=1:r.nchannels
    chanlocs(i).labels=char(eeg.label(i));
    chanlocs(i).theta=0;
    chanlocs(i).radius=0;
    chanlocs(i).X=0;
    chanlocs(i).Y=0;
    chanlocs(i).Z=0;
    chanlocs(i).sph_theta=0;
    chanlocs(i).sph_phi=0;
    chanlocs(i).sph_radius=0;
end
r.chanlocs=chanlocs;

if r.sample1
    s1=r.sample1;
else
    s1=((r.time1 * r.rate)+1);
    r.sample1=s1;
end

if r.sample2
    s2=r.sample2;
else
    if r.time2
        s2=((r.time2 * r.rate)+1);
        if (s2 > r.totsamples)|(s2 < s1); s2=r.totsamples; end
    else
        s2=r.totsamples;
    end
end

eeg=read_eep_cnt(file,s1,s2);
r.nsmpl=eeg.npnt;
r.time=eeg.time;
r.dat=eeg.data;

author: PatternRecognition

repo: OpenBMI


text:
% eql_sps_os_test.m
% compare aspire and matlab versions of E-QPL-SPS-OS
% Copyright Apr 2000, Jeff Fessler, The University of Michigan

%
% generate data
%
if ~isvar('yi'), printm 'data'
	if has_aspire
		f.dir	= test_dir;
		f.wtf	= [f.dir 't,g.wtf'];
		f.wtr	= strrep(f.wtf, 'wtf', 'wtr');
		f.yi	= [f.dir 'yi.fld'];
		f.ci	= [f.dir 'ci.fld'];
		f.ri	= [f.dir 'ri.fld'];
		f.mask	= [f.dir 'mask.fld'];
	end
	em_test_setup
prompt
end

if ~isvar('Gb'), printm 'Gb'
	f.nblock = 5;
	Gb = Gblock(G, f.nblock);
prompt
end

if ~isvar('R'), printm 'R'
	f.nbrs = 8;
	f.l2b = 2;
	R = Robject(ig.mask, 'edge_type', 'leak', 'type_denom', 'aspire', ...
		'beta', 2^f.l2b);
prompt
end

%
% matlab iterations
%
if ~isvar('xmat'), printm 'matlab E-QPL-SPS-OS'
	f.niter = 9;
	f.pixmax = 6;
	curvs = {'pc', 'oc'};
	xinit = max(xfbp,0);
	for ic = 1:length(curvs)
		tmp = eql_sps_os(xinit(ig.mask), Gb, yi, ci, ri, R, ...
			f.niter, f.pixmax, curvs{ic});
		xmat{ic} = ig.embed(tmp);
	end
	im clf, im(xmat{1}, 'Matlab E-QPL-SPS-OS iterations')
prompt
end

if ~has_aspire, return, end

%
% aspire iterations
%
if ~isvar('xasp'), printm 'aspire E-QPL-SPS-OS'

	f.init	= [f.dir 'init.fld'];
	f.out	= [f.dir 'out.fld'];
	fld_write(f.init, xinit, 'check', 0)

	f.saver	= 'stack,1';
	f.fitype = ['2z@' f.wtr '@-'];

	for ic = 1:length(curvs)
		if exist(f.out, 'file'), delete(f.out), end
		f.alg = sprintf('ospsc,%s,%d,%d,1,0', ...
			curvs{ic}, f.nblock, sg.na);
		f.penal	= sprintf('%g,quad,%d,-', f.l2b, f.nbrs/4);
		f.method = sprintf('@%d@%s@%s', f.niter-1, f.alg, f.penal);
		f.com = sprintf(['i -chat 0 empl3 %s %s  %s %s 1 %s 1 %s -' ...
				' %s %s 0 1 %g 0 -'], ...
			f.out, f.init, f.yi, f.ci, f.ri, f.fitype, ...
			f.method, f.saver, f.pixmax);

		os_run(f.com)
		xasp{ic} = double(fld_read(f.out));
	end
end

if 1
	for ic = 1:length(xasp)
		t = vcorrcoef(xasp{ic}, xmat{ic});
		printf('corr. %g,%g', t, t-1)
		err = (xasp{ic} - xmat{ic}) / max(col(xmat{ic}));
		printf('Normalized error range %g %g', min(err(:)), max(err(:))),
	end

	im clf, im(221, xmat{1}, 'xhat matlab'), cbar
	im(222, xasp{1}, 'xhat aspire'), cbar
	im(223, (xasp{1}-xmat{1})/max(col(xmat{1})), 'aspire-matlab'), cbar

	t1 = eql_obj(xmat{1}, G, yi(:), ci(:), ri(:), R, ig.mask);
	t2 = eql_obj(xasp{1}, G, yi(:), ci(:), ri(:), R, ig.mask);

	if im
		subplot(224)
		plot(0:f.niter-1, t1-t1(1), '-o', 0:f.niter-1, t2-t1(1), '-x')
		xlabel iteration, ylabel '\Phi change', legend('mat', 'asp', 4)
		title(sprintf('E-QPL-SPS-OS, Nsubset=%d', f.nblock))
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [C, sigma] = dataset3Params(X, y, Xval, yval)
%EX6PARAMS returns your choice of C and sigma for Part 3 of the exercise
%where you select the optimal (C, sigma) learning parameters to use for SVM
%with RBF kernel
%   [C, sigma] = EX6PARAMS(X, y, Xval, yval) returns your choice of C and
%   sigma. You should complete this function to return the optimal C and
%   sigma based on a cross-validation set.
%

% You need to return the following variables correctly.
C = 1;
sigma = 0.3;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the optimal C and sigma
%               learning parameters found using the cross validation set.
%               You can use svmPredict to predict the labels on the cross
%               validation set. For example,
%                   predictions = svmPredict(model, Xval);
%               will return the predictions on the cross validation set.
%
%  Note: You can compute the prediction error using
%        mean(double(predictions ~= yval))
%

results = eye(64,3);
errorRow = 0;

for C_test = [0.01 0.03 0.1 0.3 1, 3, 10 30]
    for sigma_test = [0.01 0.03 0.1 0.3 1, 3, 10 30]
        errorRow = errorRow + 1;
        model = svmTrain(X, y, C_test, @(x1, x2) gaussianKernel(x1, x2, sigma_test));
        predictions = svmPredict(model, Xval);
        prediction_error = mean(double(predictions ~= yval));

        results(errorRow,:) = [C_test, sigma_test, prediction_error];
    end
end

sorted_results = sortrows(results, 3); % sort matrix by column #3, the error, ascending

C = sorted_results(1,1);
sigma = sorted_results(1,2);

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

end

author: Borye


text:
function Stall_margin = retreating_stall_envelope(mu_act,Cts_act)
% Datum Rotor Coefficients (mu verus Cw/s)      High Speed Rotor Coefficients (mu verus Cw/s)
% a:	-7.6568	Squared term                    a:	-7.4029	Squared term
% b:	0.2992	Linear term                     b:	1.2392	Linear term
% c:	0.5357	Constant                        c:	0.4827	Constant

highspeed = 0; %0: datum; 1: high speed

% Datum Rotor Coefficients (mu verus Cw/s)
Cdat=[-7.6568   ;
0.2992;
0.5357];

% High Speed Rotor Coefficients (mu verus Cw/s)
Chs=[-7.4029	;
1.2392	;
0.4827	];

C = highspeed*Chs+(1-highspeed)*Cdat;
a = C(1); b = C(2); c = C(3);

Cts_act=Cts_act*2; %because of the way to find Ct in the place where this data comes from (with a 1/2 in Ct calc)

mu_all = a*Cts_act.^2+Cts_act*b+c;

Cts_all = (-b-sqrt(b^2-4*a*(c-mu_act)))/(2*a);

Stall_margin = (mu_all-mu_act)/abs(mu_all);
if isreal(Cts_all)
    Ct_margin = (Cts_all-Cts_act)/abs(Cts_all);
    Stall_margin = (Stall_margin+Ct_margin)/2;
else
    Stall_margin = -abs(Stall_margin);
end

author: Sable


text:
function [X, iter] = fista(grad, Xinit, L, lambda, opts, calc_F)
% function [X, iter] = fista(grad, Xinit, L, lambda, opts, calc_F)
% * A Fast Iterative Shrinkage-Thresholding Algorithm for
% Linear Inverse Problems.
% * Solve the problem: `X = arg min_X F(X) = f(X) + lambda||X||_1` where:
%   - `X`: variable, can be a matrix.
%   - `f(X)` is a smooth convex function with continuously differentiable
%       with Lipschitz continuous gradient `L(f)` (Lipschitz constant of
%       the gradient of `f`).
% * Syntax: `[X, iter] = FISTA(calc_F, grad, Xinit, L, lambda,  opts)` where:
%   - INPUT:
%     + `grad`: a _function_ calculating gradient of `f(X)` given `X`.
%     + `Xinit`: initial guess.
%     + `L`: the Lipschitz constant of the gradient of `f(X)`.
%     + `lambda`: a regularization parameter, can be either a scalar or
%       a weighted matrix.
%     + `opts`: a _structure_ variable describing the algorithm.
%       * `opts.max_iter`: maximum iterations of the algorithm.
%           Default `300`.
%       * `opts.tol`: a tolerance, the algorithm will stop if difference
%           between two successive `X` is smaller than this value.
%           Default `1e-8`.
%       * `opts.show_progress`: showing `F(X)` after each iteration or not.
%           Default `false`.
%     + `calc_F`: optional, a _function_ calculating value of `F` at `X`
%       via `feval(calc_F, X)`.
% -------------------------------------
% Author: Tiep Vu, thv102, 4/6/2016
% (http://www.personal.psu.edu/thv102/)
% -------------------------------------
    opts = initOpts(opts);
    Linv = 1/L;
    lambdaLiv = lambda*Linv;
    x_old = Xinit;
    y_old = Xinit;
    t_old = 1;
    iter = 0;
    cost_old = 1e10;
    %% MAIN LOOP
    while  iter < opts.max_iter
        iter = iter + 1;
        x_new = shrinkage(y_old - Linv*feval(grad, y_old), lambdaLiv);
        t_new = 0.5*(1 + sqrt(1 + 4*t_old^2));
        y_new = x_new + (t_old - 1)/t_new * (x_new - x_old);
        %% check stop criteria
        e = norm1(x_new - x_old)/numel(x_new);
        if e < opts.tol
            break;
        end
        %% update
        x_old = x_new;
        t_old = t_new;
        y_old = y_new;
        %% show progress
        if opts.verbose
            if nargin ~= 0
                cost_new = feval(calc_F, x_new);
%                 if cost_new <= cost_old
%                     stt = 'YES.';
%                 else
%                     stt = 'NO, check your code.';
%                 end
%                 fprintf('iter = %3d, cost = %f, cost decreases? %s\n', ...
%                     iter, cost_new, stt);
                fprintf('iter = %3d, cost = %f\n', iter, cost_new);
                cost_old = cost_new;
            else
                if mod(iter, 5) == 0
                    fprintf('.');
                end
                if mod(iter, 10) == 0
                   fprintf('%d', iter);
                end
            end
        end
    end
    X = x_new;
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
% EX_ARTICLEV3_EXAMPLE51: simple example to generate multipatch objects of a two-patch geometry.
%
% The geometry is given by the two squares, [0,1] x [0,1] and [1,2] x [0,1],
%  defined with biquadratic splines.
%
% The code generates the multipatch msh and space, as in Example 5.1 of the paper:
%
% R. Vazquez, A new design for the implementation of isogeometric analysis
%   in Octave and Matlab: GeoPDEs 3.0. Tech, Report, IMATI-CNR, 2016
%
% Copyright (C) 2016 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

nrb(1) = nrbdegelev (nrb4surf ([0 0], [1 0], [0 1], [1 1]), [1 1]);
nrb(2) = nrbtform (nrb(1), vectrans([1 0 0]));

[geometry, boundaries, interfaces, ~, boundary_interfaces] = mp_geo_load (nrb);
npatch = numel (geometry);
for iptc = 1:npatch
  knots = geometry(iptc).nurbs.knots;
  [qn, qw] = msh_set_quad_nodes (knots, msh_gauss_nodes (geometry(iptc).nurbs.order));

  local_meshes{iptc} = msh_cartesian (knots, qn, qw, geometry(iptc), 'der2', false);
  local_spaces{iptc}  = sp_nurbs (geometry(iptc).nurbs, local_meshes{iptc});
end

msh = msh_multipatch (local_meshes, boundaries);
space = sp_multipatch (local_spaces, msh, interfaces, boundary_interfaces);

disp ('Local to global numbering in space.gnum')
disp(['PATCH 1: ', num2str(space.gnum{1}.')])
disp(['PATCH 2: ', num2str(space.gnum{2}.')])
fprintf('\n')

disp ('Boundary information in msh.boundary')
disp (['patch_numbers: ', num2str(msh.boundary.patch_numbers)])
disp (['side_numbers:  ', num2str(msh.boundary.side_numbers)])
fprintf('\n')

disp (['Number of boundary patches, in space.boundary.npatch: ', num2str(space.boundary.npatch)]);
fprintf('\n')
disp ('Local to global numbering in space.boundary.gnum:')
for ii = 1:space.boundary.npatch
  disp(['PATCH ',num2str(ii),':  ', num2str(space.boundary.gnum{ii}.')])
end
fprintf('\n')

disp ('Boundary to volumetric numbering, in space.boundary.dofs')
disp (num2str(space.boundary.dofs.'))

fprintf('\n')
disp ('For the solution of problems in multipatch geometries, see also mp_solve_laplace.m')

%!test
%! nrb(1) = nrbdegelev (nrb4surf ([0 0], [1 0], [0 1], [1 1]), [1 1]);
%! nrb(2) = nrbtform (nrb(1), vectrans([1 0 0]));
%! [geometry, boundaries, interfaces, ~, boundary_interfaces] = mp_geo_load (nrb);
%! npatch = numel (geometry);
%! for iptc = 1:npatch
%!   knots = geometry(iptc).nurbs.knots;
%!   [qn, qw] = msh_set_quad_nodes (knots, msh_gauss_nodes (geometry(iptc).nurbs.order));
%!   local_meshes{iptc} = msh_cartesian (knots, qn, qw, geometry(iptc), 'der2', false);
%!   local_spaces{iptc}  = sp_nurbs (geometry(iptc).nurbs, local_meshes{iptc});
%! end
%! msh = msh_multipatch (local_meshes, boundaries);
%! space = sp_multipatch (local_spaces, msh, interfaces, boundary_interfaces);
%! assert (space.gnum{1}(:), [1 2 13 3 4 14 5 6 15].');
%! assert (space.gnum{2}(:), [13 7 8 14 9 10 15 11 12].');
%! assert (msh.boundary.patch_numbers, [1 1 1 2 2 2]);
%! assert (msh.boundary.side_numbers, [1 3 4 2 3 4]);
%! assert (space.boundary.npatch, 6)
%! assert (space.boundary.gnum{1}(:), [7 1 8].');
%! assert (space.boundary.gnum{2}(:), [7 2 9].');
%! assert (space.boundary.gnum{3}(:), [8 3 10].');
%! assert (space.boundary.gnum{4}(:), [11 4 12].');
%! assert (space.boundary.gnum{5}(:), [9 5 11].');
%! assert (space.boundary.gnum{6}(:), [10 6 12].');
%! assert (space.boundary.dofs(:), [3 2 6 10 7 11 1 5 13 15 8 12].');

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [z, sd, m] = cellzscore(x, dim, flag)

% [Z, SD] = CELLZSCORE(X, DIM, FLAG) computes the zscore, across all cells in x along
% the dimension dim, normalising by the total number of samples
%
% X should be an linear cell-array of matrices for which the size in at
% least one of the dimensions should be the same for all cells. If flag==1, the mean will
% be subtracted first (default behavior, but to save time on already demeaned data, it
% can be set to 0). SD is a vector containing the standard deviations, used for the normalisation.

nx = size(x);
if ~iscell(x) || length(nx)>2 || all(nx>1),
  error('incorrect input for cellstd');
end

if nargin<2,
  scx1 = cellfun('size', x, 1);
  scx2 = cellfun('size', x, 2);
  if     all(scx2==scx2(1)), dim = 2; %let second dimension prevail
  elseif all(scx1==scx1(1)), dim = 1;
  else   error('no dimension to compute mean for');
  end
elseif nargin==2,
  flag = 1;
end

if flag,
  m    = cellmean(x, dim);
  x    = cellvecadd(x, -m);
end

sd   = cellstd(x, dim, 0);
z    = cellvecmult(x, 1./sd);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% p = genpath_visible(d)
% save as genpath, but if any of the path is hidden, the directory is not
% added
function p = genpath_visible(d)

hiddenchar = '.';

p = '';
a = genpath(d);
if isempty(a),
  return;
end

% find the first character of each path
i0s = [0,find(a == pathsep)]+1;
if i0s(end) <= length(a),
  i0s(end+1) = length(a)+1;
end
n = length(i0s)-1;
off = length(d) + 1;

for i = 1:n,
  i0 = i0s(i);
  i1 = i0s(i+1)-2;
  if i0+off > i1,
    continue;
  end
  ishidden = a(i0+off) == hiddenchar || ...
    ~isempty(strfind(a(i0+off:i1),[filesep,hiddenchar]));
  if ~ishidden,
    p = [p,a(i0:i1),pathsep]; %#ok<AGROW>
    %fprintf('Added: %s\n',a(i0:i1));
    %else
    %fprintf('Hidden: %s\n',a(i0:i1));
  end
end

author: kristinbranson

repo: JAABA


text:
function L = symbchol()
%                                                          L = symbchol(X)
%
% SYMBCHOL Symbolic block sparse Cholesky factorization.
%   L = symbchol(X) returns a structure L that can be used
%   by the efficient block sparse Cholesky solver SPARCHOL.
%   The fields in L have the following meaning:
%
%   L.perm   - Multiple minimum degree ordering.
%
%   L.L      -  Sparse lower triangular matrix, has sparsity structure
%     of Cholesky factor of X(L.perm,L.perm).
%
%   L.xsuper - Supernode partition. Supernode jsup consists of
%     the nodes   L.xsuper(jsup) : L.xsuper(jsup)-1.
%
%   L.split  - Splitting of supernodes. Recommends to split supernode
%     in blocks of sizes   L.split(xsuper(jsup):L.xsuper(jsup)-1).
%
%   L.tmpsiz - Quantity used by SPARCHOL, to allocated enough working
%     storage.
%
%   L = symbchol(X,cachsz) optimizes L.split for a computer cache
%     of size CACHSZ * 1024 byte. Default cachsz = 512.
%
% See also sparchol, sparfwslv, sparbwslv, symbfact, symmmd, chol.

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

% ----------------------------------------
% Enter here the cache-size in KB, for shaping
% optimal dense blocks of floats.
% ----------------------------------------
global ADA_sedumi_
if ~issparse(ADA_sedumi_)
    error('X should be a sparse symmetric matrix')
end
cachsz = 512;
% ----------------------------------------
% Compute multiple minimum degree ordering.
% If the matrix is actually dense we don't bother.
% ----------------------------------------
if spars(ADA_sedumi_)<1
    perm = ordmmdmex(ADA_sedumi_);
    L = symfctmex(ADA_sedumi_,perm);
else
    L.perm=(1:size(ADA_sedumi_,1))';
    L.L=sparse(tril(ones(size(ADA_sedumi_))));
    L.xsuper=[1;size(ADA_sedumi_,1)+1];
end
% ----------------------------------------
% Symbolic Cholesky factorization structures, stored in L.
% ----------------------------------------
L.tmpsiz = choltmpsiz(L);
L.split = cholsplit(L,cachsz);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function fileNames = list_image_files(parameterString)

if nargin==0,
    %clc; clear;
    parameterString = '.';
end

temp = dir(parameterString);
fileNames = {};
[fileNames{1:length(temp),1}] = deal(temp.name);

% Find hidden files.
hiddenFiles = nan(size(fileNames));
for fileIndex = 1:length(fileNames),
    tempName = fileNames{fileIndex};
    hiddenFiles(fileIndex) = tempName(1) == '.';
end

% Find which entries in the listing are images.
imageFlags = nan(size(fileNames));
for fileIndex = 1:length(fileNames),
    [boo1, boo2, tempExtension] = fileparts(fileNames{fileIndex});
    tempExtension = lower(tempExtension);
    imageFlags(fileIndex) = ismember(tempExtension, {'.bmp', '.png', '.tif', '.jpg'});
end

% Find which entries in the listing are directories.
directoryFlags = {};
[directoryFlags{1:length(temp),1}] = deal(temp.isdir);
directoryFlags = cell2mat(directoryFlags);

tempFlags = (~hiddenFiles) & (~directoryFlags) & (imageFlags);
fileNames = fileNames(tempFlags);

repo: lrslibrary


text:
function state= util_warning(msg, id, varargin)
%UTIL_WARNING - Display warning message
%
%Description:
% Essentially this function is just a work-around such that one can
% use the increased functionality of the warning function of Matlab R>=13,
% while having programs that still run under Matlab R<13.
% Furthermore, this function allows to suppress warnings which follow each
% other within a short interval of time.
%
%Usage:
% STATE= util_warning(MSG, ID, <OPT>)
%
%Input:
% MSG:  (string) the warning message, or 'off' or 'query', see the
%       help of 'warning'.
% ID:   (string)
% OPT:  (struct(proplist) Property/value list or struct of optional parameters:
%   'Interval': Warnings that appear earlier than this interval of time [s]
%        after a warning (for the same ID) are suppressed. Default 10*60.
%   'File': Display this as the file that issued the warning.
%        Default: name of the calling function.
%
%
%Example 1 (within a matlab function):
% util_warning('outer model selection sucks', 'validation');
%
%Example 2:
% wstat= util_warning('off', 'validation');
% %% ... some code admittedly producing 'validation' warnings ...
% util_warning(wstat);

persistent id_list last_warning

st= dbstack;
if length(st)>1,
  default_File= st(2).name;
else
  default_File= '';
end
props= {'File'       default_File   'CHAR';
        'Interval'   10*60          '!DOUBLE[1]';
       };
opt= opt_proplistToStruct(varargin{:});
opt= opt_setDefaults(opt, props, 1);

misc_checkType(msg, 'CHAR');
misc_checkType(id, 'CHAR');

this_warning= clock;
id_idx= strmatch(msg, id_list, 'exact');
if isempty(id_idx),
  id_list= cat(2, id_list, {msg});
  last_warning= [last_warning; this_warning];
else
  if etime(this_warning, last_warning(id_idx,:)) < opt.Interval,
    return;
  end
  last_warning(id_idx,:)= this_warning;
end

if isempty(opt.File),
  ff= '';
else
  ff= [opt.File ': '];
end

a= sscanf(getfield(ver('MATLAB'), 'Release'), '(R%d)');
if a>=13,
  if isstruct(msg),
    warning(msg);
  elseif strcmp(msg, 'off'),
    if nargout>0,
      state= warning('query', ['BBCI:' id]);
    end
    warning('off', ['bbci:' id]);
  elseif strcmp(msg, 'query'),
    state= warning('query', ['BBCI:' id]);
  else
    warning(['bbci:' id], [ff msg]);
  end
else
  if ~isempty(msg) & isempty(strmatch(msg, {'off','query'})),
    state= warning([ff msg]);
  else
    state= [];
  end
end

author: bbci


text:
function add_ndt_paths_and_init_rand_generator

% This function adds all the appropriate directories for the NDT toolbox to work
%   and initializes the random number generator based on the clock (by default
%   Matlab uses the same seed everything it is started which can lead to unwanted behavior).

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

%     This code is part of the Neural Decoding Toolbox.
%     Copyright (C) 2011 by Ethan Meyers (emeyers@mit.edu)
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

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

% file the path where the toolbox is
toolbox_basedir_name = [fileparts(which('add_ndt_paths_and_init_rand_generator')) '/']

% add all the appropriate sub-directories
addpath([toolbox_basedir_name 'datasources/']);
addpath([toolbox_basedir_name 'feature_preprocessors/']);
addpath([toolbox_basedir_name 'classifiers/']);
addpath([toolbox_basedir_name 'cross_validators/']);
addpath([toolbox_basedir_name 'tools/']);
addpath([toolbox_basedir_name 'helper_functions/']);
addpath(genpath([toolbox_basedir_name 'external_libraries/']));  % include subdirectories of any external libraries

if isOctave
    addpath([toolbox_basedir_name 'octave_code/']);
    'Running Octave - adding Octave specific code to the search path'
end

% initialize the random number generator
rand('state',sum(100*clock));
disp('initializing the matlab random number generator to an aribrary clock value, i.e., the ''rand'' function should now work properly');

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ G ] = gsp_erdos_renyi( N,p,param )
%GSP_ERDOS_RENYI Create a random Erdos Renyi graph
%   Usage:  G = gsp_erdos_renyi( N,p,param );
%           G = gsp_erdos_renyi( N,p );
%
%   Input parameters:
%         N     : Number of nodes
%         p     : Probability of connection of a node with another
%         param : Structure of optional parameter
%   Output parameters:
%         G     : Graph structure.
%
%   'gsp_erdos_renyi( N,p,param )' initializes create a graph structure
%   containing the weighted adjacency matrix (G.W), the number of vertices
%   (G.N) for an Erdos Renyi graph. All edge weights are equal to 1.
%
%   The Erdos Renyi graph is constructed by connecting nodes randomly. Each
%   edge is included in the graph with probability p independent from every
%   other edge.
%
%   *param* a Matlab structure containing the following fields:
%
%     * *param.connected* : flag to force the graph to be connected. By
%       default, it is 0.
%
%     * *param.maxit* : is the maximum number of try to connect the graph.
%       By default, it is 10.
%
%     * *param.verbose* : 0 no log, 1 print main steps, 2 print all steps.
%       By default, it is 1.
%
%   Example::
%
%        G = gsp_erdos_renyi(100,0.05)

% Author : Nathanael Perraudin

% Optional input arguments
if nargin<3, param=struct; end

if ~isfield(param, 'connected'), param.connected=0; end
if ~isfield(param, 'maxit'), param.maxit=10; end
if ~isfield(param, 'verbose'), param.verbose=1 ; end

% Check if the user didn't put silly parameter
if param.connected
    if N*p<3
        if param.verbose
            fprintf(['The model has very low chance to create a '...
                ,'connected graph. Increase N, p or disable '...
                ,'param.connected \n']);
        end
    end
end

if p>1
   error('The probability p cannot be above 1');
end

if p<0
    error('The probability p cannot be negative');
end

% Create the graph structure
G=struct;
G.graph_type='erdos_renyi';
G.connected=param.connected;

% Create the graph
if param.connected
    bool=0;
    iter=1;
    if param.verbose>1
        fprintf('Entering the iteration process\n')
    end
    while ~bool

        G.W = sprandsym(N, p)>0;
        G.W(1:N+1:end) = 0;

        bool=gsp_check_connectivity(G.W);
        if iter==param.maxit
            warning('The graph is not strongly connected.');
            break;
        elseif ~bool
            if param.verbose>1
                fprintf('Iteration %i failed. Trying again \n',iter)
            end
            iter=iter+1;
        end
    end
    if param.verbose>1
        if bool
            fprintf('A connected graph has been created in %i iteration\n',iter);
        else
            fprintf('The iterative process has failed.\n');
        end
    end
else
   G.W = sprandsym(N,p)>0;
   G.W(1:N+1:end) = 0;
end

% try to create point to be associated with the points of the graph

% [x, y] = draw_dot(G.W);
%
%
% G.coords=[x',y'];
% G.limits=[-1e-4,1.01*max(x),-1e-4,1.01*max(y)];

G.type = 'erdos_renyi';

G = gsp_graph_default_parameters(G);

end

repo: gspbox


text:
filename='Bridge_quad_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'PROJECTED GRADIENT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.2;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%% Machine Learning Online Class - Exercise 2: Logistic Regression
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the second part
%  of the exercise which covers regularization with logistic regression.
%
%  You will need to complete the following functions in this exericse:
%
%     sigmoid.m
%     costFunction.m
%     predict.m
%     costFunctionReg.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% Load Data
%  The first two columns contains the X values and the third column
%  contains the label (y).

data = load('ex2data2.txt');
X = data(:, [1, 2]); y = data(:, 3);

plotData(X, y);

% Put some labels
hold on;

% Labels and Legend
xlabel('Microchip Test 1')
ylabel('Microchip Test 2')

% Specified in plot order
legend('y = 1', 'y = 0')
hold off;

%% =========== Part 1: Regularized Logistic Regression ============
%  In this part, you are given a dataset with data points that are not
%  linearly separable. However, you would still like to use logistic
%  regression to classify the data points.
%
%  To do so, you introduce more features to use -- in particular, you add
%  polynomial features to our data matrix (similar to polynomial
%  regression).
%

% Add Polynomial Features

% Note that mapFeature also adds a column of ones for us, so the intercept
% term is handled
X = mapFeature(X(:,1), X(:,2));

% Initialize fitting parameters
initial_theta = zeros(size(X, 2), 1);

% Set regularization parameter lambda to 1
lambda = 1;

% Compute and display initial cost and gradient for regularized logistic
% regression
[cost, grad] = costFunctionReg(initial_theta, X, y, lambda);

fprintf('Cost at initial theta (zeros): %f\n', cost);

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============= Part 2: Regularization and Accuracies =============
%  Optional Exercise:
%  In this part, you will get to try different values of lambda and
%  see how regularization affects the decision coundart
%
%  Try the following values of lambda (0, 1, 10, 100).
%
%  How does the decision boundary change when you vary lambda? How does
%  the training set accuracy vary?
%
% Initialize fitting parameters
initial_theta = zeros(size(X, 2), 1);

% Set regularization parameter lambda to 1 (you should vary this)
lambda = 1;
%lambda = 0;
%lambda = 10;
%lambda = 100;

% Set Options
options = optimset('GradObj', 'on', 'MaxIter', 400);

% Optimize
[theta, J, exit_flag] = ...
	fminunc(@(t)(costFunctionReg(t, X, y, lambda)), initial_theta, options);

% Plot Boundary
plotDecisionBoundary(theta, X, y);
hold on;
title(sprintf('lambda = %g', lambda))

% Labels and Legend
xlabel('Microchip Test 1')
ylabel('Microchip Test 2')

legend('y = 1', 'y = 0', 'Decision boundary')
hold off;

% Compute accuracy on our training set
p = predict(theta, X);

fprintf('Train Accuracy: %f\n', mean(double(p == y)) * 100);

repo: MachineLearning


text:
function seed = get_seed ( )

%*****************************************************************************80
%
%% GET_SEED returns a seed for the random number generator.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer SEED, a random seed value.
%
  I_MAX = 2147483647;

  time_array = clock;

  hour = time_array(4);
  minute = time_array(5);
  second = time_array(6);

  temp = ( second + 60 * ( minute + 60 * hour ) ) / ( 60.0 * 60.0 * 24.0 );

  if ( temp <= 0.0 )
    temp = temp + 1.0;
  end

  if ( 1.0 < temp )
    temp = temp - 1.0;
  end

  seed = 1 + floor ( I_MAX * temp );

  return
end

author: johannesgerer

repo: jburkardt-m


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 start the server on the KUKA iiwa controller
% Then run this script using Matlab

% Copyright: Mohammad SAFEEA, 13th of March 2018

% Important: Be careful when runnning the script, be sure that no human, nor obstacles
% are around the robot

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.5; % 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


text:
function hmObj = hlp_validateHeadModelObject(hmObj)
% helper function to check the head model object, and
% load from base workspace if necessary
% Author: Tim Mullen, SCCN/INC/UCSD 2013

if isempty(hmObj)
    hmObj = 'hmObj'; % default name of variable in base wkspace
end
if ischar(hmObj)
    % load variable from base workspace
    try
       hmObj = evalin('base',hmObj);
    catch err1
        % try to open it as a file names
        try
            hmObj = hlp_microcache('headmodels',@headModel.loadFromFile,env_translatepath(hmObj));
        catch err2
            if ~strcmp(err2.identifier,'MATLAB:UndefinedFunction')
                disp_once(['Note: The head model named "' hmObj '" was neither found in the base workspace nor could it be loaded as a file.']); end
            hmObj = '';
            return;
        end
    end
end

% check that the object is a valid head model object
if ~isobject(hmObj)
    warning('You did not supply a valid head model object');
    hmObj = '';
    return;
end
if ~isprop(hmObj,'leadFieldFile') || isempty(hmObj.leadFieldFile)
    warning('Head model object does not appear to be valid. The lead field matrix does not exist.');
    hmObj = '';
    return;
end


text:
function test_pull433

% MEM 2gb
% WALLTIME 00:20:00
% DEPENDENCY

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/pull433'));
load('SubjectUCI29_elec_tal_f.mat', 'elec_tal_f');
load('SubjectUCI29_hull_lh.mat', 'mesh');

% keepchannel
elec_tal_fr = elec_tal_f;
grids = {'LPG*', 'LTG*'};
for g = 1:numel(grids)
  cfg             = [];
  cfg.channel     = grids{g};
  cfg.keepchannel = 'yes';            % <-
  cfg.elec        = elec_tal_fr;
  cfg.method      = 'headshape';
  cfg.headshape   = mesh;
  cfg.warp        = 'dykstra2012';
  elec_tal_fr = ft_electroderealign(cfg);
end
assert(isequal(numel(elec_tal_fr.label), 152)); % realigned + unused is 152 elecs

% do not keepchannel
elec_tal_fr = elec_tal_f;
grids = {'LPG*'};
cfg             = [];
cfg.channel     = grids{1};
cfg.keepchannel = 'no';             % <-
cfg.elec        = elec_tal_fr;
cfg.method      = 'headshape';
cfg.headshape   = mesh;
cfg.warp        = 'dykstra2012';
elec_tal_fr = ft_electroderealign(cfg);
assert(isequal(numel(elec_tal_fr.label), 64)); % realigned LPG has 64 elecs


text:
% IM = pgmRead( FILENAME )
%
% Load a pgm image into a MatLab matrix.
%   This format is accessible from the XV image browsing utility.
%   Only works for 8bit gray images (raw or ascii)

% Hany Farid, Spring '96.  Modified by Eero Simoncelli, 6/96.

function im = pgmRead( fname );

[fid,msg] = fopen( fname, 'r' );

if (fid == -1)
  error(msg);
end

%%% First line contains ID string:
%%% "P1" = ascii bitmap, "P2" = ascii greymap,
%%% "P3" = ascii pixmap, "P4" = raw bitmap,
%%% "P5" = raw greymap, "P6" = raw pixmap
TheLine = fgetl(fid);
format  = TheLine;

if ~((format(1:2) == 'P2') | (format(1:2) == 'P5'))
  error('PGM file must be of type P2 or P5');
end

%%% Any number of comment lines
TheLine  = fgetl(fid);
while TheLine(1) == '#'
	TheLine = fgetl(fid);
end

%%% dimensions
sz = sscanf(TheLine,'%d',2);
xdim = sz(1);
ydim = sz(2);
sz = xdim * ydim;

%%% Maximum pixel value
TheLine  = fgetl(fid);
maxval = sscanf(TheLine, '%d',1);

%%im  = zeros(dim,1);
if (format(2) == '2')
  [im,count]  = fscanf(fid,'%d',sz);
else
  [im,count]  = fread(fid,sz,'uchar');
end

fclose(fid);

if (count == sz)
  im = reshape( im, xdim, ydim )';
else
  fprintf(1,'Warning: File ended early!');
  im = reshape( [im ; zeros(sz-count,1)], xdim, ydim)';
end

repo: SelfExSR


text:
function [x]=mesh
%------------------------------------------------
% GENERATE MATCHING DISTANCE X FOR DIFFERENT MESH
% -----------------------------------------------
global N L
dx=L/N;
x(1)=0;
x(2)=dx/2;
for i=3:N+1
    x(i)=x(i-1)+dx;
end
x(N+2)=L;

repo: mcbench-benchmarks


text:
function a = ge( x, y )

b = newcnstr( evalin( 'caller', 'cvx_problem', '[]' ), x, y, '>=' );
if nargout, a = b; end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function dbfprintf(varargin)
    if ~ZmapGlobal.Data.debug
        return
    end
    if isstring(varargin{1})
        fmtstr = char(varargin{1});
    else
        fmtstr = varargin{1};
    end

    if nargin>1
        fprintf(['\n---<strong> ZMAP DB: </strong>' fmtstr], varargin{2:end});
    else
        fprintf(['\n---<strong> ZMAP DB: </strong>' fmtstr]);
    end
end

author: CelsoReyes


text:
function [resid, resid_norm, Wxk,Hkt,centers,Vxtpix, Vxtpixbg] = computeresid(res,dpixc,peval)
% [resid, resid_norm, Wxk,Hkt,centers,Vxtpix, Vxtpixbg] =
% computeresid(res,dpixc,peval)
[Wxk,Hkt,centers,Vxtpix]=reshapeGaP(res.hvec,res.cxcy,peval);
Vxtpixbg=reshape(Wxk*Hkt,peval.nx,peval.ny,peval.nt)+peval.bg;
resid=(Vxtpixbg-dpixc);
resid_norm=resid./sqrt(Vxtpixbg);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function hb_to_st_test ( )

%*****************************************************************************80
%
%% HB_TO_ST_TEST tests the HB_TO_ST library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HB_TO_ST_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the HB_TO_ST library.\n' );

  hb_to_st_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HB_TO_ST_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function hb_to_st_test01 ( )

%*****************************************************************************80
%
%% HB_TO_ST_TEST01 tests HB_TO_ST.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
  input_file = 'rua_32.hb';
  output_file = 'rua_32.st';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HB_TO_ST_TEST01\n' );
  fprintf ( 1, '  We will read a Harwell-Boeing sparse matrix file\n' );
  fprintf ( 1, '  called "%s"\n', input_file );
  fprintf ( 1, '  and have HB_TO_ST convert it to an ST file.\n' );

  hb_to_st ( input_file, output_file );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function Q = alex2jaakko(A)
Q.W = A.A;
Q.X = A.S;
Q.mu = A.Mu;

if ~isempty(A.Av)
  D = rows(A.Av{1});
  M = length(A.Av);
  Q.CovW = zeros([D,D,M]);
  for m=1:M
    Q.CovW(:,:,m) = A.Av{m};
  end
end
if ~isempty(A.Sv)
  D = rows(A.Sv{1});
  N = length(A.Sv);
  Q.CovX = zeros([D,D,N]);
  for n=1:N
    Q.CovW(:,:,n) = A.Sv{n};
  end
end

if isempty(A.Muv)
  Q.v_mu = 0;
else
  Q.v_mu = A.Muv;
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function b_n = beamWeightsTorus2Spherical(N)
%BEAMWEIGHTSTORUS2SPHERICAL Generate beamweights for a raised torus
%
%   The tabulated beamweights below describe toroidal patterns of the form
%   D(theta) = |sin(\theta)|^N, using up to 4th-order components.
%   Note that N=2 and N=4 are exactly band-limited to 2nd and 4th-order
%   respectively, while the N=1 and N=3 extend to higher-orders with
%   decaying energy, and hence the weights below are just an approximation.
%   Because the pattern is axisymmetric only the N+1 coefficients of m=0
%   are returned.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% BEAMWEIGHTSTORUS2SPHERICAL.M - 11/7/2013
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch N
    case 1
        b_n = [2.7842    0.0000   -0.7781    0.0000   -0.1303];
    case 2
        b_n = [2.3633   -0.0000   -1.0569];
    case 3
        b_n = [2.0881   -0.0000   -1.1673    0.0000    0.1468];
    case 4
        b_n = [1.8906   -0.0000   -1.2079    0.0000    0.2701];
end

author: polarch

repo: Spherical-Array-Processing


text:
classdef SubCellConnecComputer < handle

    properties (GetAccess = public, SetAccess = private)
        connec
    end

    properties (Access = private)
        coord
        nElem
        isEdgeCutInElem
        levelSet
        vertexNodesInElem

        allNodesInElem

        firstCutEdge
        finalNode
        cutNodesInElem

        isSubCellInterior
        triangleSubMesher
    end

    methods (Access = public)

        function obj = SubCellConnecComputer(cParams)
            obj.init(cParams);
            obj.compute();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.firstCutEdge      = cParams.firstCutEdge;
            obj.vertexNodesInElem = cParams.backgroundConnec;
            obj.nElem             = cParams.nElem;
            obj.isEdgeCutInElem   = cParams.isEdgeCutInElem;
            obj.coord             = cParams.coord;
            obj.levelSet          = cParams.levelSet;
            obj.finalNode         = cParams.finalNode;

        end

        function compute(obj)
            obj.computeCutNodePerElem();
            obj.computeAllNodesInElem();
            obj.computeNodesInAllSubCells();
            obj.computeIsSubCellsInterior();
            obj.computeConnec();
        end

        function computeAllNodesInElem(obj)
            vertexNodes = obj.vertexNodesInElem;
            cutNodes    = obj.cutNodesInElem;
            allNodes    = [vertexNodes,cutNodes];
            obj.allNodesInElem = allNodes;
        end

        function computeNodesInAllSubCells(obj)
            s.nodesInElem     = obj.allNodesInElem;
            s.isEdgeCutInElem = obj.isEdgeCutInElem;
            s.nElem           = obj.nElem;
            s.coord           = obj.coord;
            tComputer = TriangleSubCellNodesComputer(s);
            tComputer.compute();
            obj.triangleSubMesher = tComputer;
        end

        function computeIsSubCellsInterior(obj)
            isTriInt  = obj.computeIsSubCellTriangleInterior();
            nSubCells = obj.triangleSubMesher.nSubCellsByElem;
            itIs = false(nSubCells,obj.nElem);
            itIs(1,isTriInt)  = true;
            itIs(2,~isTriInt) = true;
            itIs(3,~isTriInt) = true;
            obj.isSubCellInterior = itIs;
        end

        function itIs = computeIsSubCellTriangleInterior(obj)
            isoNode = obj.triangleSubMesher.isoNode;
            isoNodeIsFull = obj.levelSet(isoNode) < 0;
            itIs = isoNodeIsFull;
        end

        function computeConnec(obj)
            allConnec  = obj.triangleSubMesher.allSubCellsConnec;
            isInterior = obj.isSubCellInterior(:);
            obj.connec = allConnec(isInterior,:);
        end

        function cutNodePerElemen = computeCutNodePerElem(obj)
            firstCutEdgePerElem = obj.firstCutEdge;
            cutNodePerElemen = firstCutEdgePerElem + obj.finalNode;
            obj.cutNodesInElem = cutNodePerElemen;
        end

    end

end


text:
function [x,fval,exitflag,info,Opt] = opti_bintprog(f,A,b,Aeq,beq,x0,opts)
%OPTI_BINTPROG Solve a BIP using an OPTI BIP Solver (Matlab Overload)
%
%   [x,fval,exitflag,info] = opti_bintprog(f,A,b,Aeq,beq,x0) solves the
%   linear program min f'x where A,b are the inequality constraints,
%   Aeq,beq are the equality constraints and x0 is the initial guess. All
%   decision variables are binary only.
%
%   [x,fval,exitflag,info] = opti_bintprog(f,...,x0,opts) allows the user
%   to specify optiset options. This includes specifying a solver via the
%   'solver' field of optiset.
%
%   [x,...,info,Opt] = opti_bintprog(f,...) returns the internally built
%   OPTI object.

%   Copyright (C) 2011 Jonathan Currie (I2C2)

% Handle missing arguments
if nargin < 7, opts = optiset; end
if nargin < 6, x0 = []; end
if nargin < 5, beq = []; end
if nargin < 4, Aeq = []; end
if nargin < 3, error('You must supply at least 3 arguments to opti_bintprog'); end

%Build OPTI Object
Opt = opti('f',f,'ineq',A,b,'eq',Aeq,beq,'x0',x0,'int',repmat('B',size(f)),'options',opts);

%Solve
[x,fval,exitflag,info] = solve(Opt);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
classdef MaOEACSS < ALGORITHM
% <multi/many> <real/integer/label/binary/permutation>
% Many-objective evolutionary algorithms based on coordinated selection
% strategy
% t --- 0 --- Threshold value in environmental selection

%------------------------------- Reference --------------------------------
% Z. He and G. G. Yen, Many-objective evolutionary algorithms based on
% coordinated selection strategy, IEEE Transactions on Evolutionary
% Computation, 2017, 21(2): 220-233.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            t = Algorithm.ParameterSet(0);

            %% Generate random population
            Population = Problem.Initialization();
            Zmin       = min(Population.objs,[],1);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = MatingSelection(Population.objs,Zmin);
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                Zmin       = min([Zmin;Offspring.objs],[],1);
                Population = EnvironmentalSelection([Population,Offspring],Zmin,t,Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO