text:
% Usage: [xopt, fmin, retcode] = nlopt_minimize(algorithm, f, f_data, lb, ub,
%                                               xinit, stop)
%
% Minimizes a nonlinear multivariable function f(x, f_data{:}), where
% x is a row vector, returning the optimal x found (xopt) along with
% the minimum function value (fmin = f(xopt)) and a return code (retcode).
% A variety of local and global optimization algorithms can be used,
% as specified by the algorithm parameter described below.  lb and ub
% are row vectors giving the upper and lower bounds on x, xinit is
% a row vector giving the initial guess for x, and stop is a struct
% containing termination conditions (see below).
%
% This function is a front-end for the external routine nlopt_minimize
% in the free NLopt nonlinear-optimization library, which is a wrapper
% around a number of free/open-source optimization subroutines.  More
% details can be found on the NLopt web page (ab-initio.mit.edu/nlopt)
% and also under 'man nlopt_minimize' on Unix.
%
% f should be a handle (@) to a function of the form:
%
%    [val, gradient] = f(x, ...)
%
% where x is a row vector, val is the function value f(x), and gradient
% is a row vector giving the gradient of the function with respect to x.
% The gradient is only used for gradient-based optimization algorithms;
% some of the algorithms (below) are derivative-free and only require
% f to return val (its value).  f can take additional arguments (...)
% which are passed via the argument f_data: f_data is a cell array
% of the additional arguments to pass to f.  (Recall that cell arrays
% are specified by curly brackets { ... }.  For example, pass f_data={}
% for functions that require no additional arguments.)
%
% stop describes the termination criteria, and is a struct with a
% number of optional fields:
%     stop.ftol_rel = fractional tolerance on function value
%     stop.ftol_abs = absolute tolerance on function value
%     stop.xtol_rel = fractional tolerance on x
%     stop.xtol_abs = row vector of absolute tolerances on x components
%     stop.fmin_max = stop when f < fmin_max is found
%     stop.maxeval = maximum number of function evaluations
%     stop.maxtime = maximum run time in seconds
%     stop.verbose = > 0 indicates verbose output
% Minimization stops when any one of these conditions is met; any
% condition that is omitted from stop will be ignored.  WARNING:
% not all algorithms interpret the stopping criteria in exactly the
% same way, and in any case ftol/xtol specify only a crude estimate
% for the accuracy of the minimum function value/x.
%
% The algorithm should be one of the following constants (name and
% interpretation are the same as for the C function).  Names with
% _G*_ are global optimization, and names with _L*_ are local
% optimization.  Names with _*N_ are derivative-free, while names
% with _*D_ are gradient-based algorithms.  Algorithms:
%
% NLOPT_GD_MLSL_LDS, NLOPT_GD_MLSL, NLOPT_GD_STOGO, NLOPT_GD_STOGO_RAND,
% NLOPT_GN_CRS2_LM, NLOPT_GN_DIRECT_L, NLOPT_GN_DIRECT_L_NOSCAL,
% NLOPT_GN_DIRECT_L_RAND, NLOPT_GN_DIRECT_L_RAND_NOSCAL, NLOPT_GN_DIRECT,
% NLOPT_GN_DIRECT_NOSCAL, NLOPT_GN_ISRES, NLOPT_GN_MLSL_LDS, NLOPT_GN_MLSL,
% NLOPT_GN_ORIG_DIRECT_L, NLOPT_GN_ORIG_DIRECT, NLOPT_LD_AUGLAG_EQ,
% NLOPT_LD_AUGLAG, NLOPT_LD_LBFGS, NLOPT_LD_LBFGS_NOCEDAL, NLOPT_LD_MMA,
% NLOPT_LD_TNEWTON, NLOPT_LD_TNEWTON_PRECOND,
% NLOPT_LD_TNEWTON_PRECOND_RESTART, NLOPT_LD_TNEWTON_RESTART,
% NLOPT_LD_VAR1, NLOPT_LD_VAR2, NLOPT_LN_AUGLAG_EQ, NLOPT_LN_AUGLAG,
% NLOPT_LN_BOBYQA, NLOPT_LN_COBYLA, NLOPT_LN_NELDERMEAD,
% NLOPT_LN_NEWUOA_BOUND, NLOPT_LN_NEWUOA, NLOPT_LN_PRAXIS, NLOPT_LN_SBPLX
%
% For more information on individual algorithms, see their individual
% help pages (e.g. "help NLOPT_LN_SBPLX").
function [xopt, fmin, retcode] = nlopt_minimize(algorithm, f, f_data, lb, ub, xinit, stop)

  [xopt, fmin, retcode] = nlopt_minimize_constrained(algorithm, f, f_data, {}, {}, lb, ub, xinit, stop);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [t0, x0, u0] = shift(T, t0, x0, u,f)
st = x0;
con = u(1,:)';
f_value = f(st,con);
st = st+ (T*f_value);
x0 = full(st);

t0 = t0 + T;
u0 = [u(2:size(u,1),:);u(size(u,1),:)];
end

author: MMehrez


text:
function cordic_test ( )

%*****************************************************************************80
%
%% CORDIC_TEST tests the CORDIC library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 June 2007
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CORDIC_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the CORDIC library.\n' );

  cordic_test001 ( );
  cordic_test002 ( );
  cordic_test003 ( );
  cordic_test004 ( );
  cordic_test005 ( );
  cordic_test006 ( );
  cordic_test007 ( );
  cordic_test008 ( );
  cordic_test009 ( );
  cordic_test010 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CORDIC_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f = dtlz2a(x)

alpha = 1;

%% g function
sum = 0;
for i = 3:8
sum = sum + (x(:,i)-0.5)^2;
end
g = sum;

%% cost functions
f(:,1) = (1+g)*cos(((x(:,1)^alpha)*pi)/2)*cos(((x(:,2)^alpha)*pi)/2);
f(:,2) = (1+g)*cos(((x(:,1)^alpha)*pi)/2)*sin(((x(:,2)^alpha)*pi)/2);
f(:,3) = (1+g)*sin(((x(:,1)^alpha)*pi)/2);
end

author: Eric-Bradford

repo: TS-EMO


text:
function gui_components(rval_space,rval_time,max_pr,sizeA,A,Cn,options)

        figure;CC = plot_contours(A,Cn,options,0); title('Selected components','fontweight','bold','fontsize',14); close all;

        plt_cont = @(x) plot_contours(A,Cn,options,0,[],CC,[],x);

        d1 = options.d1;
        d2 = options.d2;

        figure;
            screensize = get(0,'Screensize' );
            fac = min(min((screensize(3:4)-200)./[2*d2,1.2*d1]),10);
            fig_size = round([100 100 fac*2*d2 fac*1.2*d1]);
            set(gcf, 'Position', fig_size);

            imagesc(Cn); axis equal; axis tight; axis off;
                set(gca,'Position',[0.01,0.1,0.5,0.8]);

        thr_space = @(x,y,z,w,v) find((rval_space >= x) & (rval_time >= y) & (sizeA >= z) & (sizeA < w) & (log(1-max_pr) < v));

        hsl1 = uicontrol('Style','slider','Min',0,'Max',1,...
                        'SliderStep',[0.01,0.01],'Value',options.space_thresh,'Tag','space_corr',...
                        'Position',[round(0.55*fig_size(3)) 50 20 round(0.85*fig_size(4))]);
        txt1 = uicontrol('Style','text','Position',[round(0.55*fig_size(3))-50 50+round(0.85*fig_size(4)) 100 50],'String','space threshold' );
        %ved1 = uicontrol('Style','edit','Position',[round(0.55*fig_size(3))-50 50+round(0.85*fig_size(4)) 100 50],'String',['val=',num2str(x1)] );

        hsl2 = uicontrol('Style','slider','Min',0,'Max',1,...
                        'SliderStep',[0.01,0.01],'Value',options.time_thresh,'Tag','time_corr',...
                        'Position',[round(0.6*fig_size(3)) 50 20 round(0.85*fig_size(4))]);
        txt2 = uicontrol('Style','text','Position',[round(0.6*fig_size(3))-50 50+round(0.85*fig_size(4)) 100 50],'String','time threshold' );

        hsl3 = uicontrol('Style','slider','Min',0,'Max',120,...
                        'SliderStep',[0.01,0.01],'Value',options.min_size_thr,'Tag','min_size',...
                        'Position',[round(0.65*fig_size(3)) 50 20 round(0.85*fig_size(4))]);
        txt3 = uicontrol('Style','text','Position',[round(0.65*fig_size(3))-50 50+round(0.85*fig_size(4)) 100 50],'String','minimum size' );

        hsl4 = uicontrol('Style','slider','Min',0,'Max',500,...
                        'SliderStep',[.01,.05],'Value',options.max_size_thr,'Tag','max_size',...
                        'Position',[round(0.7*fig_size(3)) 50 20 round(0.85*fig_size(4))]);

        txt4 = uicontrol('Style','text','Position',[round(0.7*fig_size(3))-50 50+round(0.85*fig_size(4)) 100 50],'String','maximum size' );

        hsl5 = uicontrol('Style','slider','Min',-10,'Max',0,...
                        'SliderStep',[0.01,0.01],'Value',log10(1-options.max_pr_thr),'Tag','max_pr',...
                        'Position',[round(0.75*fig_size(3)) 50 20 round(0.85*fig_size(4))]);

        txt5 = uicontrol('Style','text','Position',[round(0.75*fig_size(3))-50 50+round(0.85*fig_size(4)) 100 50],'String','log(1-max_pr)' );

        pb = uicontrol('Style','pushbutton','Callback',@close_Callback);

        set(hsl1,'Callback',@(hObject,eventdata) plt_cont(thr_space(get(hObject,'Value'),hsl2.Value,hsl3.Value,hsl4.Value,hsl5.Value)))
        set(hsl2,'Callback',@(hObject,eventdata) plt_cont(thr_space(hsl1.Value,get(hObject,'Value'),hsl3.Value,hsl4.Value,hsl5.Value)))
        set(hsl3,'Callback',@(hObject,eventdata) plt_cont(thr_space(hsl1.Value,hsl2.Value,get(hObject,'Value'),hsl4.Value,hsl5.Value)))
        set(hsl4,'Callback',@(hObject,eventdata) plt_cont(thr_space(hsl1.Value,hsl2.Value,hsl3.Value,get(hObject,'Value'),hsl5.Value)))
        set(hsl5,'Callback',@(hObject,eventdata) plt_cont(thr_space(hsl1.Value,hsl2.Value,hsl3.Value,hsl4.Value,get(hObject,'Value'))))

%     function close_Callback(hObject, eventdata)
%
%             close all;
%     end

end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
 function sr = de_ftab_fit_sprad(fit, sl)
%function sr = de_ftab_fit_sprad(fit, sl)
% for each s vector compute spectra radius of hessian of LS fit
% in
%	sl	cell{L}	[(Nd),L]
% out
%	sr	[(Nd),M]
% applies only to 'exp' fit

sll = ndgrid_jf('mat', sl{:}); % [(Nd),L]
Nd = size(sll); LL = Nd(end); Nd = Nd(1:end-1);
sll = reshape(sll, [], LL); % [*Nd,L]

if LL ~= 2, error 'only L=2 done', end
MM = fit.MM;

%sr = cell(fit.MM,1);
sr = zeros(prod(Nd), MM);
for mm=1:MM
	alf = fit.coef{mm}; % [ne,1]
	mac = fit.mac{mm}; % [ne,L]
	g0 = mac' * alf; % gradient of fit at s=0 (largest point)
	h0 = mac' * diag(alf) * mac - g0 * g0'
	norm(g0)
	norm(g0 * g0')
	norm(h0)
	for is=1:prod(Nd)
		ss = sll(is,:)'; % [L,1]
		q = alf .* exp(-mac * ss); % [ne,1]
		q = q / sum(q);
%		plot(q), drawnow
		g = mac' * q; % gradient of f_m(s)
		h = mac' * diag(q) * mac - g * g';
		sr(is,mm) = norm(g * g');
%		sr(is,mm) = norm(h);
	end
end
sr = reshape(sr, [Nd MM]);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function  AIMIDA  = readAIMIDA(fid)
%
% <<<<<<<<<< CLASSIFICATION: UNCLASSIFIED >>>>>>>>>>
%
% READAIMIDA Create structure for metadata fields.
%   AIMIDA = READAIMIDA(FID) returns a structure of fields and their
%   associated values.

% Written by Matt Donath, NGA, matthew.b.donath@nga.ic.gov
% References:
%  > MIL-STD-2500A, NATIONAL IMAGERY TRANSMISSION FORMAT VERSION 2.0
%  > MIL-STD-2500C, NATIONAL IMAGERY TRANSMISSION FORMAT VERSION 2.1
%  > STDI-0001,     NATIONAL SUPPORT DATA EXTENSIONS (SDE)(VERSION 1.3/CN2)
%                   FOR THE NATIONAL IMAGERY TRANSMISSION FORMAT (NITF)
%  > STDI-0002,     THE COMPENDIUM OF CONTROLLED EXTENSIONS (CE) FOR THE
%                   NATIONAL IMAGERY TRANSMISSION FORMAT (NITF) VERSION 2.1
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Create the metadata structure for AIMIDA.

AIMIDA.CETAG = fread(fid,6,'uint8=>char')';
AIMIDA.CEL = str2double(fread(fid,5,'uint8=>char')');

AIMIDA.MISSION_DATE = fread(fid,7,'uint8=>char')';
AIMIDA.MISSION_NO = fread(fid,4,'uint8=>char')';
AIMIDA.FLIGHT_NO = fread(fid,2,'uint8=>char')';
AIMIDA.OP_NUM = fread(fid,3,'uint8=>char')';
AIMIDA.START_SEGMENT = fread(fid,2,'uint8=>char')';
AIMIDA.REPRO_NUM = str2double(fread(fid,2,'uint8=>char')');
AIMIDA.REPLAY = fread(fid,3,'uint8=>char')';
AIMIDA.Reserved1 = fread(fid,1,'uint8=>char')';
AIMIDA.START_COLUMN = str2double(fread(fid,2,'uint8=>char')');
AIMIDA.START_ROW = str2double(fread(fid,5,'uint8=>char')');
AIMIDA.END_SEGMENT = fread(fid,2,'uint8=>char')';
AIMIDA.END_COLUMN = str2double(fread(fid,2,'uint8=>char')');
AIMIDA.END_ROW = str2double(fread(fid,5,'uint8=>char')');
AIMIDA.COUNTRY = fread(fid,2,'uint8=>char')';
AIMIDA.Reserved2 = fread(fid,4,'uint8=>char')';
AIMIDA.LOCATION = fread(fid,11,'uint8=>char')';
AIMIDA.TIME = fread(fid,5,'uint8=>char')';
AIMIDA.CREATION_DATE = fread(fid,7,'uint8=>char')';

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function pred = RunInference (factors)
% This function performs inference for a Markov network specified as a list
% of factors.
%
% Input:
%   factors: An array of struct factors, each containing 'var', 'card', and
%     'val' fields.
%
% Output:
%   pred: An array of predictions for every variable. In particular,
%     pred(i) is the predicted value for variable numbered i (as determined
%     by the 'var' fields in the input factors).
%
% Copyright (C) Daphne Koller, Stanford University, 2012

binaries = {'.\inference\doinference.exe', ...
            './inference/doinference-mac', ...
            './inference/doinference-linux'};

kFactorsFilename = 'factors.fg';
kStderrFilename = 'inf.log';
kInfBinary = binaries{[ispc ismac isunix]}; % NB: need ismac first so that if ismac and isunix are both 1, then mac is chosen
kInferenceType = 'map'; % choices are 'map' or 'pd'

factorsString = SerializeFactorsFg (factors);

fd = fopen(kFactorsFilename, 'wt');
fprintf (fd, '%s', factorsString);
fclose(fd);

if (isunix && ~ismac)
  command = [kInfBinary ' ' kFactorsFilename ' ' kInferenceType];
else
  command = [kInfBinary ' ' kFactorsFilename ' ' kInferenceType ' 2> ' kStderrFilename];
end

[retVal, output] = system(command);

if (retVal ~= 0)
    error('The doinference command failed. Look at the file %s to diagnose the cause', kStderrFilename);
end

pred = ParseOutput(output);

end

function pred = ParseOutput(output)

lines = strread(output, '%s', 'delimiter', sprintf('\n'));
lines(strcmp(lines, '')) = [];

numVars = str2double(lines{1});

if (numVars ~= length(lines) - 1)
    error('Error parsing output: %s', output);
end

pred = zeros(numVars, 1);
for i = 2:(numVars + 1)
    line = str2num(lines{i}); %#ok
    pred(i-1) = line(end);
end

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
  function [B Bgx Bgy Bgz] = makeB(ig, kg, varargin)
%|function [B Bgx Bgy Bgz] = makeB(ig, kg, [?])
%|
%| Create a Fatrix B operation of warping based on cubic B-spline
%|
%| in:
%|	'ig'		image geometry
%|	'kg'		knot geometry
%|
%| out:
%|	B		fatrix operation for warping
%|	Bgx		fatrix operation for warping
%|	Bgy		fatrix operation for warping
%|	Bgz		fatrix operation for warping
%|
%| Copyright August 2006, Se Young Chun and Jeff Fessler, University of Michigan

if nargin < 2, ir_usage, end

argB.power = 0;
argB.ig = ig;
argB.kg = kg;
[argB.kernelx argB.kernelgx] = makeKernel(3, kg.mx);
[argB.kernely argB.kernelgy] = makeKernel(3, kg.my);
if (kg.is3 == 1)
	[argB.kernelz argB.kernelgz] = makeKernel(3, kg.mz);
end

if (nargin == 2)
	B = Fatrix([ig.np kg.np], argB, 'forw', @makeB_forward, ...
		'back', @makeB_transpose, 'power', @makeB_power);
	Bgx = Fatrix([ig.np kg.np], argB, 'forw', @makeBgx_forward, ...
		'back', @makeBgx_transpose);
	Bgy = Fatrix([ig.np kg.np], argB, 'forw', @makeBgy_forward, ...
		'back', @makeBgy_transpose);
	if (kg.is3 == 1)
		Bgz = Fatrix([ig.np kg.np], argB, 'forw', @makeBgz_forward, ...
			'back', @makeBgz_transpose);
	end

else
	if (kg.is3 == 1)
		argB.WarpCell = varargin(1:3);
		Bgz = Fatrix([ig.np kg.np],argB,'forw',@makeBgz_forwardwarp, ...
			'back', @makeBgz_transposewarp);
	else
		argB.WarpCell = varargin(1:2);
	end
	B = Fatrix([ig.np kg.np], argB, 'forw', @makeB_forwardwarp, ...
		'back', @makeB_transposewarp);
	Bgx = Fatrix([ig.np kg.np], argB, 'forw', @makeBgx_forwardwarp, ...
		'back', @makeBgx_transposewarp);
	Bgy = Fatrix([ig.np kg.np], argB, 'forw', @makeBgy_forwardwarp, ...
		'back', @makeBgy_transposewarp);
end

%%%
function DFM = makeB_forward(argB, alpha)
alpha = double(reshape(alpha, argB.kg.dim));
if (argB.kg.is3 == 0) % 2D case
	DFM = BsplCo2ValZeroFilt(alpha, [argB.ig.nx, argB.ig.ny], ...
		[argB.kg.offset_x argB.kg.offset_y], ...
		[argB.kg.mx argB.kg.my], {argB.kernelx, argB.kernely});
else
	DFM = BsplCo2ValZeroFilt(alpha, [argB.ig.nx, argB.ig.ny argB.ig.nz], ...
		[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
		[argB.kg.mx argB.kg.my argB.kg.mz], ...
		{argB.kernelx, argB.kernely, argB.kernelz});
end
DFM = single(DFM(:));

%%%
function ALP = makeB_transpose(argB, imgco)
imgco = double(reshape(imgco, argB.ig.dim));
if (argB.kg.is3 == 0) % 2D case
	ALP = BsplCo2ValTranZeroFilt(imgco, [argB.kg.nx argB.kg.ny], ...
		[argB.kg.offset_x argB.kg.offset_y], ...
		[argB.kg.mx argB.kg.my], {argB.kernelx, argB.kernely});
else
	ALP = BsplCo2ValTranZeroFilt(imgco,[argB.kg.nx argB.kg.ny argB.kg.nz],...
		[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
		[argB.kg.mx argB.kg.my argB.kg.mz], ...
		{argB.kernelx, argB.kernely, argB.kernelz});
end
ALP = single(ALP(:));

%%%
function DFM = makeB_forwardwarp(argB, alpha)
alpha = single(reshape(alpha, argB.kg.dim));
if (argB.kg.is3 == 0) % 2D case
        DFM = BsplCo2ValZero(alpha, [argB.ig.nx, argB.ig.ny], ...
		[argB.kg.offset_x argB.kg.offset_y], ...
		[argB.kg.mx argB.kg.my], argB.WarpCell);
else
        DFM = BsplCo2ValZero(alpha, [argB.ig.nx, argB.ig.ny argB.ig.nz], ...
		[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
		[argB.kg.mx argB.kg.my argB.kg.mz], argB.WarpCell);
end
DFM = DFM(:);

%%%
function ALP = makeB_transposewarp(argB, imgco)
imgco = single(reshape(imgco, argB.ig.dim));
if (argB.kg.is3 == 0) % 2D case
	ALP = BsplCo2ValTranZero(imgco, [argB.kg.nx, argB.kg.ny], ...
		[argB.kg.offset_x argB.kg.offset_y], ...
		[argB.kg.mx argB.kg.my], argB.WarpCell);
else
	ALP = BsplCo2ValTranZero(imgco, [argB.kg.nx, argB.kg.ny argB.kg.nz], ...
		[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
		[argB.kg.mx argB.kg.my argB.kg.mz], argB.WarpCell);
end
ALP = ALP(:);

%%%
function DFM = makeBgx_forward(argB, alpha)
alpha = double(reshape(alpha, argB.kg.dim));
if (argB.kg.is3 == 0) % 2D case
	DFM = BsplCo2ValZeroFilt(alpha, [argB.ig.nx, argB.ig.ny], ...
		[argB.kg.offset_x argB.kg.offset_y], ...
		[argB.kg.mx argB.kg.my], {argB.kernelgx, argB.kernely});
else
	DFM = BsplCo2ValZeroFilt(alpha, [argB.ig.nx, argB.ig.ny argB.ig.nz], ...
		[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
		[argB.kg.mx argB.kg.my argB.kg.mz], ...
		{argB.kernelgx, argB.kernely, argB.kernelz});
end
DFM = single(DFM(:));

%%%
function ALP = makeBgx_transpose(argB, imgco)
imgco = double(reshape(imgco, argB.ig.dim));
if (argB.kg.is3 == 0) % 2D case
	ALP = BsplCo2ValTranZeroFilt(imgco, [argB.kg.nx argB.kg.ny], ...
		[argB.kg.offset_x argB.kg.offset_y], ...
		[argB.kg.mx argB.kg.my], ...
		{fliplr(argB.kernelgx), fliplr(argB.kernely)});
else
	ALP = BsplCo2ValTranZeroFilt(imgco,[argB.kg.nx argB.kg.ny argB.kg.nz],...
		[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
		[argB.kg.mx argB.kg.my argB.kg.mz], ...
		{fliplr(argB.kernelgx), fliplr(argB.kernely), ...
		fliplr(argB.kernelz)});
end
ALP = single(ALP(:));

%%%
function DFM = makeBgx_forwardwarp(argB, alpha)
alpha = single(reshape(alpha, argB.kg.dim));
if (argB.kg.is3 == 0) % 2D case
        DFM = BsplCo2GdXZero(alpha, [argB.ig.nx, argB.ig.ny], ...
		[argB.kg.offset_x argB.kg.offset_y], ...
		[argB.kg.mx argB.kg.my], argB.WarpCell);
else
        DFM = BsplCo2GdXZero(alpha, [argB.ig.nx, argB.ig.ny argB.ig.nz], ...
		[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
		[argB.kg.mx argB.kg.my argB.kg.mz], argB.WarpCell);
end
DFM = DFM(:);

%%%
function ALP = makeBgx_transposewarp(argB, imgco)
imgco = single(reshape(imgco, argB.ig.dim));
if (argB.kg.is3 == 0) % 2D case
	ALP = BsplCo2GdXTranZero(imgco, [argB.kg.nx, argB.kg.ny], ...
		[argB.kg.offset_x argB.kg.offset_y], ...
		[argB.kg.mx argB.kg.my], argB.WarpCell);
else
	ALP = BsplCo2GdXTranZero(imgco, [argB.kg.nx, argB.kg.ny argB.kg.nz], ...
		[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
		[argB.kg.mx argB.kg.my argB.kg.mz], argB.WarpCell);
end
ALP = ALP(:);

%%%
function DFM = makeBgy_forward(argB, alpha)
alpha = double(reshape(alpha, argB.kg.dim));
if (argB.kg.is3 == 0) % 2D case
	DFM = BsplCo2ValZeroFilt(alpha, [argB.ig.nx, argB.ig.ny], ...
		[argB.kg.offset_x argB.kg.offset_y], ...
		[argB.kg.mx argB.kg.my], {argB.kernelx, argB.kernelgy});
else
	DFM = BsplCo2ValZeroFilt(alpha, [argB.ig.nx, argB.ig.ny argB.ig.nz], ...
		[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
		[argB.kg.mx argB.kg.my argB.kg.mz], ...
		{argB.kernelx, argB.kernelgy, argB.kernelz});
end
DFM = single(DFM(:));

%%%
function ALP = makeBgy_transpose(argB, imgco)
imgco = double(reshape(imgco, argB.ig.dim));
if (argB.kg.is3 == 0) % 2D case
	ALP = BsplCo2ValTranZeroFilt(imgco, [argB.kg.nx argB.kg.ny], ...
		[argB.kg.offset_x argB.kg.offset_y], ...
		[argB.kg.mx argB.kg.my], ...
		{fliplr(argB.kernelx), fliplr(argB.kernelgy)});
else
	ALP = BsplCo2ValTranZeroFilt(imgco,[argB.kg.nx argB.kg.ny argB.kg.nz],...
		[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
		[argB.kg.mx argB.kg.my argB.kg.mz], ...
		{fliplr(argB.kernelx), fliplr(argB.kernelgy), ...
		fliplr(argB.kernelz)});
end
ALP = single(ALP(:));

%%%
function DFM = makeBgy_forwardwarp(argB, alpha)
alpha = single(reshape(alpha, argB.kg.dim));
if (argB.kg.is3 == 0) % 2D case
        DFM = BsplCo2GdYZero(alpha, [argB.ig.nx, argB.ig.ny], ...
		[argB.kg.offset_x argB.kg.offset_y], ...
		[argB.kg.mx argB.kg.my], argB.WarpCell);
else
        DFM = BsplCo2GdYZero(alpha, [argB.ig.nx, argB.ig.ny argB.ig.nz], ...
		[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
		[argB.kg.mx argB.kg.my argB.kg.mz], argB.WarpCell);
end
DFM = DFM(:);

%%%
function ALP = makeBgy_transposewarp(argB, imgco)
imgco = single(reshape(imgco, argB.ig.dim));
if (argB.kg.is3 == 0) % 2D case
	ALP = BsplCo2GdYTranZero(imgco, [argB.kg.nx, argB.kg.ny], ...
		[argB.kg.offset_x argB.kg.offset_y], ...
		[argB.kg.mx argB.kg.my], argB.WarpCell);
else
	ALP = BsplCo2GdYTranZero(imgco, [argB.kg.nx, argB.kg.ny argB.kg.nz], ...
		[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
		[argB.kg.mx argB.kg.my argB.kg.mz], argB.WarpCell);
end
ALP = ALP(:);

%%%
function DFM = makeBgz_forward(argB, alpha)
alpha = double(reshape(alpha, argB.kg.dim));
DFM = BsplCo2ValZeroFilt(alpha, [argB.ig.nx, argB.ig.ny argB.ig.nz], ...
	[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
	[argB.kg.mx argB.kg.my argB.kg.mz], ...
	{argB.kernelx, argB.kernely, argB.kernelgz});
DFM = single(DFM(:));

%%%
function ALP = makeBgz_transpose(argB, imgco)
imgco = double(reshape(imgco, argB.ig.dim));
ALP = BsplCo2ValTranZeroFilt(imgco,[argB.kg.nx argB.kg.ny argB.kg.nz],...
	[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
	[argB.kg.mx argB.kg.my argB.kg.mz], ...
	{fliplr(argB.kernelx), fliplr(argB.kernely), ...
	fliplr(argB.kernelgz)});
ALP = single(ALP(:));

%%%
function DFM = makeBgz_forwardwarp(argB, alpha)
alpha = single(reshape(alpha, argB.kg.dim));
DFM = BsplCo2GdZZero(alpha, [argB.ig.nx, argB.ig.ny argB.ig.nz], ...
	[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
	[argB.kg.mx argB.kg.my argB.kg.mz], argB.WarpCell);
DFM = DFM(:);

%%%
function ALP = makeBgz_transposewarp(argB, imgco)
imgco = single(reshape(imgco, argB.ig.dim));
ALP = BsplCo2GdZTranZero(imgco, [argB.kg.nx, argB.kg.ny argB.kg.nz], ...
	[argB.kg.offset_x argB.kg.offset_y argB.kg.offset_z], ...
	[argB.kg.mx argB.kg.my argB.kg.mz], argB.WarpCell);
ALP = ALP(:);

%%%
function ob = makeB_power(ob, sup)

if (sup == 1)

elseif ((sup == 2) && (ob.arg.power == 0))
	ob.arg.power = 1;
else
	error('only squares of each element are supported');
end
ob = ob;


text:
function varargout = ubcmatch(varargin)
% VL_UBCMATCH  Match SIFT features
%   MATCHES = VL_UBCMATCH(DESCR1, DESCR2) matches the two sets of SIFT
%   descriptors DESCR1 and DESCR2.
%
%   [MATCHES,SCORES] = VL_UBCMATCH(DESCR1, DESCR2) retuns the matches and
%   also the squared Euclidean distance between the matches.
%
%   The function uses the algorithm suggested by D. Lowe [1] to reject
%   matches that are too ambiguous.
%
%   VL_UBCMATCH(DESCR1, DESCR2, THRESH) uses the specified threshold
%   THRESH. A descriptor D1 is matched to a descriptor D2 only if the
%   distance d(D1,D2) multiplied by THRESH is not greater than the
%   distance of D1 to all other descriptors. The default value of
%   THRESH is 1.5.
%
%   The storage class of the descriptors can be either DOUBLE, FLOAT,
%   INT8 or UINT8. Usually integer classes are faster.
%
%   REFERENCES::
%     [1] D. G. Lowe, Distinctive image features from scale-invariant
%     keypoints. IJCV, vol. 2, no. 60, pp. 91-110, 2004.
%
%   See also: VL_HELP(), VL_SIFT().
[varargout{1:nargout}] = vl_ubcmatch(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function p = NormL1_primal(x,weights)

p = norm(x.*weights,1);

author: mpf


text:
% This is a demo for segmentation using local gaussian distribution (LGD)
% fitting energy
%
% Reference: <Li Wang, Lei He, Arabinda Mishra, Chunming Li.
% Active Contours Driven by Local Gaussian Distribution Fitting Energy.
% Signal Processing, 89(12), 2009,p. 2435-2447>
%
% Please DO NOT distribute this code to anybody.
% Copyright (c) by Li Wang
%
% Author:       Li Wang
% E-mail:       li_wang@med.unc.edu
% URL:          http://www.unc.edu/~liwa/
%
% 2010-01-02 PM

clc;clear all;close all;

Img=imread('3.bmp');
Img = double(Img(:,:,1));

NumIter = 1000; %iterations
timestep=0.1; %time step
mu=0.1/timestep;% level set regularization term, please refer to "Chunming Li and et al. Level Set Evolution Without Re-initialization: A New Variational Formulation, CVPR 2005"
sigma = 3;%size of kernel
epsilon = 1;
c0 = 2; % the constant value
lambda1=1.03;%outer weight, please refer to "Chunming Li and et al,  Minimization of Region-Scalable Fitting Energy for Image Segmentation, IEEE Trans. Image Processing, vol. 17 (10), pp. 1940-1949, 2008"
lambda2=1.0;%inner weight
%if lambda1>lambda2; tend to inflate
%if lambda1<lambda2; tend to deflate
nu = 0.0005*255*255;%length term
alf = 20;%data term weight

figure,imagesc(uint8(Img),[0 255]),colormap(gray),axis off;axis equal
[Height Wide] = size(Img);
[xx yy] = meshgrid(1:Wide,1:Height);
phi = (sqrt(((xx - 64).^2 + (yy - 76).^2 )) - 15);
phi = sign(phi).*c0;

Ksigma=fspecial('gaussian',round(2*sigma)*2 + 1,sigma); %  kernel
ONE=ones(size(Img));
KONE = imfilter(ONE,Ksigma,'replicate');
KI = imfilter(Img,Ksigma,'replicate');
KI2 = imfilter(Img.^2,Ksigma,'replicate');

figure,imagesc(uint8(Img),[0 255]),colormap(gray),axis off;axis equal,
hold on,[c,h] = contour(phi,[0 0],'r','linewidth',1); hold off
pause(0.5)

tic
for iter = 1:NumIter
    phi =evolution_LGD(Img,phi,epsilon,Ksigma,KONE,KI,KI2,mu,nu,lambda1,lambda2,timestep,alf);

    if(mod(iter,50) == 0)
        figure(2),
        imagesc(uint8(Img),[0 255]),colormap(gray),axis off;axis equal,title(num2str(iter))
        hold on,[c,h] = contour(phi,[0 0],'r','linewidth',1); hold off
        pause(0.02);
    end

end
toc

author: Sable

repo: mcbench-benchmarks


text:
function result = tetra_07 ( func, x, y, z )

%*****************************************************************************80
%
%% TETRA_07 approximates an integral inside a tetrahedron in 3D.
%
%  Integration region:
%
%    Points inside a tetrahedron whose four corners are given.
%
%  Discussion:
%
%    A 64 point 7-th degree conical product Gauss formula is used,
%    Stroud number T3:7-1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 November 2004
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Arthur H Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%    Stroud and Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966, pages 42-43.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied
%    function of three variables which is to be integrated,
%    of the form:
%      function value = func ( x, y, z )
%
%    Input, real X(4), Y(4), Z(4), the X, Y and Z coordinates of
%    the vertices.
%
%    Output, real RESULT, the approximate integral of the function.
%
  norder = 4;
  weight2 = [ ...
    0.1355069134, 0.2034645680, 0.1298475476, 0.0311809709 ];
  weight3 = [ ...
    0.1108884156, 0.1434587898, 0.0686338872, 0.0103522407 ];
  xtab2 = [ ...
    0.0571041961E+00, 0.2768430136E+00, 0.5835904324E+00, 0.8602401357E+00 ];
  xtab3 = [ ...
    0.0485005495E+00, 0.2386007376E+00, 0.5170472951E+00, 0.7958514179E+00 ];
%
%  Get the Gauss-Legendre weights and abscissas for [-1,1].
%
  [ xtab1, weight1 ] = legendre_set ( norder );
%
%  Adjust the rule for the interval [0,1].
%
  a = -1.0E+00;
  b = +1.0E+00;

  c =  0.0E+00;
  d =  1.0E+00;

  [ xtab1, weight1 ] = rule_adjust ( a, b, c, d, norder, xtab1, weight1 );
%
%  Carry out the quadrature.
%
  quad = 0.0;

  for i = 1 : norder
    for j = 1 : norder
      for k = 1 : norder
%
%  Compute the barycentric coordinates of the point in the unit triangle.
%
        t =                                                 xtab3(k);
        u =                        xtab2(j)   * ( 1.0E+00 - xtab3(k) );
        v = xtab1(i) * ( 1.0E+00 - xtab2(j) ) * ( 1.0E+00 - xtab3(k) );
        w = 1.0E+00 - t - u - v;
%
%  Compute the corresponding point in the triangle.
%
        xval = t * x(1) + u * x(2) + v * x(3) + w * x(4);
        yval = t * y(1) + u * y(2) + v * y(3) + w * y(4);
        zval = t * z(1) + u * z(2) + v * z(3) + w * z(4);
        quad = quad + 6.0 * weight1(i) * weight2(j) * weight3(k) ...
          * feval ( func, xval, yval, zval );

      end
    end
  end

  volume = tetra_volume ( x, y, z );
  result = quad * volume;

  return
end


text:
function [ loglik ] = LogLikelihood_gmm(X,Priors,Mu,Sigma,w)
%LOGLIKELIHOOD_GMM
%
%   input------------------------------------------------------------------
%
%       o X : (D x N), set of N data points of dimension  D
%
%       o w : (1 x N), set of weights if the data points are weighted
%
%   output ----------------------------------------------------------------
%
%       o loglik : (N x 1) , loglikelihood
%

N = size(X,2);
K = size(Priors,2);

Pxi = zeros(N,K);
for i=1:K
    Pxi(:,i) = gaussPDF(X, Mu(:,i), Sigma(:,:,i));
end

if exist('w','var') && ~isempty(w)
    %               (N x K)
    %             (N x 1)
    F            = (repmat(w(:),1,K) .* Pxi) * Priors';
else

    %Compute the log likelihood
    F            = Pxi*Priors';

end

F(F<realmin) = realmin;
F            = Pxi*Priors';
loglik       = sum(log(F));

end

author: epfl-lasa


text:
function varargout = miqps_matpower(varargin)
%MIQPS_MATPOWER  Deprecated, please use MIQPS_MASTER instead.

%   MATPOWER
%   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

[varargout{1:nargout}] = miqps_master(varargin{:});

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [pgr] = create_grid(pts, follow_meridians, trim_final_grid_to_shape)
    % Interactively define a grid
    %
    % CREATE_GRID will create a ZmapGrid interactively. Upon closing the window or choosing "Set"
    % the global ZmapData.Grid will be updated.
    %
    % CREATE_GRID(PTS) PTS is a polygon. only grid points with PTS will be shown
    %  GR is all data points, as [lon,lat; ...]
    %  ZGR is a ZmapGrid.
    % CREATE_GRID(PTS, FOLLOW_MERIDIANS) if FOLLOW_MERIDIANS is false, then distances are constant.
    % if FOLLOW_MERIDIANS is TRUE, then degrees of longitude are constant.
    %
    % Choose an origin point
    % Choose a point to define initial spacing
    % Use scroll Wheel to change spacing
    % Drag Point to change origin
    %
    % tests:
    %  create_grid('testpoly')
    %  create_grid('testworld');
    %
    % If FOLLOW_MERIDIANS is true, then grid gets smaller as pole is approached.
    % Unfortunately, routines like pcolor this is necessary for pcolor, when Y axis is a lon.
    % This could be avoided if coordinate system was transformed into X-Y.
    %
    % When FOLLOW_MERIDIANS is false, then longitudes drift in order to keep consistent sized boxes.

    % TODO: add edit fields that allow grid to be further modified
    % TODO: add SAVE and LOAD buttons
    % TODO: make this work with gridfun

    FOLLOW_MERIDIANS = exist('follow_meridians','var') && follow_meridians;
    USEPOLY=exist('pts','var') && ~isempty(pts) && ~isnan(pts(1));
    name='grid';
    changed=false;

    ZG=ZmapGlobal.Data;

    f=figure('Name','Grid Selection','Units','pixels','Position',[200 75 730 700]);

    % DISPLAY EVENTS
    ax=subplot(4,4,[1,11]);
    ax.Units='points';
    ax.Position=fix(ax.Position);
    plot(ax,ZG.primeCatalog.Longitude,ZG.primeCatalog.Latitude,'.',...
        'color',[.75 .75 .8],'DisplayName','events');
    ax.XLabel.String='longitude';
    ax.YLabel.String='latitude';
    ax.NextPlot='add';

    if exist('pts','var') && ischar(pts)
        switch pts
            case 'testpoly'
                pts =[... % SAMPLE POLYGON, FOR TESTING
                    7.3500   47.8007;...
                    8.3187   47.5009;...
                    8.8392   46.6736;...
                    8.4633   46.0321;...
                    7.4801   45.7263;...
                    6.5548   46.7095;...
                    7.3500   47.8007...
                    ];
            case 'testworld'
                USEPOLY=false
                clear pts
                axis([0 90 -20 80])
        end
    end

    %DISPLAY SURFACE FEATURES
    copyobj(ZG.features('borders'),ax);

    % DISPLAY POLYGON
    if USEPOLY
        plot(ax, pts(:,1),pts(:,2),'k:','LineWidth',2,'DisplayName','polygon');
    end

    legend(ax,'show')

    % SHOW INSTRUCTIONS
    t = uicontrol('style','text','units','pixels','position',[10 10 300 20],'String','temporary');

    % FEEDBACK
    tp = uicontrol('style','text','Units','pixels','Position',[10 50 300 20],'String','N Points: ???');
    d = uicontrol('style','text','Units','pixels','Position',[10 30 300 20],'String','Dist: ???? (deg) [???? (km)]');

    % ADDITIONAL CONTROLS
    uicontrol('style','pushbutton','Units','pixels','Position',[310 30 50 20],'String','SET','Callback',@set_grid);
    uicontrol('style','edit','Units','pixels','Position',[310 70 50 20],...
        'String',name,'Callback',@(~,~)update_plot());

    uicontrol('style','checkbox','Units','pixels','Position',[310 100 80 20],'String','Follow Meridians? (adjusts for lat)','Callback',@toggle_parallels);
    % SELECT FIXED POINT
    write_string(t,'Enter a fixed point for the grid');
    [fixed_x,fixed_y]=ginput(1);
    fp=plot(ax,fixed_x,fixed_y,'b+','LineWidth',2,'DisplayName','Origin');
    instruction_end(t);

    % SELECT OTHER POINT
    write_string(t,'Click at a distance that will define a grid');
    [x2,y2]=ginput(1);
    %tmph=plot(ax,x2,y2,'bo','LineWidth',2)
    instruction_end(t);

    dx = abs(x2-fixed_x);
    dy = abs(y2-fixed_y);

    % SELECT INITIAL GRID
    gpts_h2=plot(nan,nan,'gx','DisplayName','grid');

    update_plot();

    f.WindowScrollWheelFcn=@adjust_grid;
    f.WindowButtonDownFcn=@mouse_down;
    f.WindowButtonUpFcn=@mouse_up;
    f.DeleteFcn=@check_for_save;

    write_string(t,'Scroll the mouse wheel to scale')
    if nargout > 0
        waitfor(f);
        pgr =ZG.Grid;
    end

    function adjust_grid(~,ev)
        % adjust grid spacing when mouse wheel is scrolled
        scale=1.05;
        %disp(ev)
        if ev.VerticalScrollCount > 0
            dx=dx.*scale;
            dy=dy.*scale;
        elseif ev.VerticalScrollCount < 0
            dx=dx./scale;
            dy=dy./scale;
        end
        update_plot();
    end

    function set_grid(~,~)
        ZG.Grid=ZmapGrid(name,pgr.xs, pgr.ys, 'deg');
        if FOLLOW_MERIDIANS
            ZG.gridopt = GridOptions(deg2km(dx),deg2km(dy),[],'deg', FOLLOW_MERIDIANS);
        else
            ZG.gridopt = GridOptions(dx, dy, [], 'km', FOLLOW_MERIDIANS);
        end
        changed=false;
    end

    function toggle_parallels(src,~)
        FOLLOW_MERIDIANS=src.Value==1;
        update_plot();
    end

    function update_plot()
        changed=true;
        [gx,gy]=get_eq_grid(fixed_x,fixed_y,dx,dy);
        if USEPOLY
            ll=polygon_filter(pts(:,1),pts(:,2),gx(:),gy(:),'inside');
        else
            ll=true(size(gx(:)));
        end
        % plot points inside polygon, but grid still covers entire map.
        gpts_h2.XData=gx(ll);
        gpts_h2.YData=gy(ll);

        pgr.xs=gx;
        pgr.ys=gy;
        disp(pgr)

        % trim pgr to polygon before returning. Maybe I shouldn't!
        if USEPOLY && exist('trim_final_grid_to_shape','var') && trim_final_grid_to_shape
            ll=polygon_filter(pts(:,1),pts(:,2),pgr.xs,pgr.ys,'inside');
            pgr.xs(~ll)=nan;
            pgr.ys(~ll)=nan;
        end

        pgr=trim_nans(pgr);

        tp.String=sprintf('N Points: %d',sum(ll));
    end

    function [lonMat,latMat] = get_eq_grid(lon0,lat0,dLon,dLat)
        % GET_EQ_GRID
        % input is the origin point and arclength between points
        % output is 2 matrices (lon, lat)

        % base grid on a single distance, so that instead of separate dx & dy, we use dd
        dist_arc = max([...
            distance(lat0,lon0,lat0,lon0+dLon,'degrees'),...
            distance(lat0,lon0,lat0+dLat,lon0,'degrees')]);

        d.String=sprintf('Dist: %.3f (deg) [%.3f (km)]',dist_arc,deg2km(dist_arc));

        % use the axes limits (assumed degrees) to control size of grid
        ylims_deg = ylim(ax);
        xlims_deg = xlim(ax);

        % pick out latitude spacing. Our grid will have this many rows.
        lats = vector_including_origin(lat0, dist_arc, ylims_deg);
        lonMat=[];
        latMat=[];

        if FOLLOW_MERIDIANS
            % when following the meridian lines, the longitude span covered by
            % the arc-distance at lat0 (along the rhumb!) remains constant.
            % that is, dLon 45 from origin (0,0) will always be 45, regardless of latitude.
            [~,dLon]=reckon('rh',lat0,0,dist_arc,90);

            % resulting in a rectangular matrix where, on a globe lines will converge, but on a graph
            lonValues = vector_including_origin(lon0, dLon, xlims_deg);

            %creates a meshgrid of size numel(lonValues) x numel(lats)
            [lonMat,latMat]=meshgrid(lonValues,lats);

        else
            % when ignoring meridian lines, and aiming for an approximately constant distance,
            % the dLon at each latitude will differ.

            % number of degrees longitude covered by the arclength at each latitude
            [~,dLon_per_lat]=reckon('rh',lats,0,dist_arc,90);

            for n=1:numel(lats)
                theseLonValues = vector_including_origin(lon0, dLon_per_lat(n), xlims_deg);
                lonMat=[lonMat;theseLonValues(:)]; %#ok<AGROW>
                latMat=[latMat;repmat(lats(n),size(theseLonValues(:)))]; %#ok<AGROW>
            end

            [lonMat,latMat] = cols2matrix(lonMat,latMat);
            % each gridx & gridy are vectors.
        end
    end

    function v = vector_including_origin(orig_deg, delta_deg, lims_deg)
        v = unique([orig_deg : -delta_deg : min(lims_deg) , orig_deg : delta_deg :max(lims_deg)]);
    end

    function [xs, ys] = cols2matrix(lonCol,latCol)
        % COLS2MATRIX convert columns of lats & lons into a matrix.
        % this takes lots of stuff into account
        %

        % assign pgrid
        ugy=unique(latCol); % lats in matrix
        nrows=numel(ugy); % number of latitudes in matrix
        [~,example]=min(abs(latCol(:))); % latitude closest to equator will have most number of lons in matrix
        mostCommonY=latCol(example); % account for the abs possibly flipping signs
        base_lon_idx=find(lonCol(latCol==mostCommonY)==fixed_x); % longitudes that must line up
        ncols=sum(latCol(:)==mostCommonY); % most number of lons in matrix
        ys=repmat(ugy(:),1,ncols);
        xs=nan(nrows,ncols);
        for n=1:nrows
            thislat=ugy(n); % lat for this row
            idx_lons=(latCol==thislat); % mask of lons in this row
            these_lons=lonCol(idx_lons); % lons in this row
            row_length=numel(these_lons); % number of lons in this row

            main_lon_idx=find(these_lons==fixed_x); % offset of X in this row
            offset=base_lon_idx - main_lon_idx;
            xs(n,(1:row_length)+offset)=these_lons;
        end
    end

    function mouse_down(src,~)
        %mouse down, so make origin follow mouse around
        %disp(ev)
        %disp(ax.CurrentPoint)
        src.WindowButtonMotionFcn=@mouse_move;
        %disp('  ');
        src.Pointer='cross';
    end

    function mouse_up(src,~)
        %mouse up, so ignore position
        %disp(ev)
        %disp(ax.CurrentPoint);
        src.WindowButtonMotionFcn='';
        %disp('  ');
        update_plot()
        src.Pointer='arrow';
    end

    function mouse_move(~,~)
        fixed_x=ax.CurrentPoint(1,1);
        fixed_y=ax.CurrentPoint(1,2);
        fp.XData=fixed_x;
        fp.YData=fixed_y;
    end

    function check_for_save(~,~)
        if changed
            dosave=questdlg('Save Changes?','Grid Creation','Yes','No','Yes');
            if strcmpi(dosave,'Yes')
                set_grid([],[]);
            end
        end
    end

end

function write_string(h,str)
    for i=1:length(str)-1
        h.String=[str(1:i),'..'];
        pause(.01);
    end
    h.String=str;
end

function instruction_end(h)
    h.ForegroundColor=[0 .5 0];
    pause(.3);
    h.String='';
    h.ForegroundColor='k';
end

function pgr=trim_nans(pgr)
    % REMOVE GRID POINTS BEYOND POLYGON (rows & cols of all nans)
    nanrows=all(isnan(pgr.xs),2);
    nancols=all(isnan(pgr.xs),1);
    pgr.xs(nanrows,:)=[];pgr.xs(:,nancols)=[];
    pgr.ys(nanrows,:)=[];pgr.ys(:,nancols)=[];
end

%{
function c=mycontext()
    c=uicontextmenu('Tag','GridContext')

    uimenu(c,'Label','Select Rectangle');
    uimenu(c,'Label','Select Circle');
    uimenu(c,'Label','Select Polygon');
    uimenu(c,'Separator','on','Label','Create Polygon');
end
%}

author: CelsoReyes


text:
 function ob = power(ob, sup)
%function ob = power(ob, sup)
% array "power" method (A.^2)
% matrix power (mpower) would be for A^2

if isempty(ob.handle_power)
	error 'no power method for this object'
end

ob = ob.handle_power(ob, sup);
if isfield(ob, 'scale') && ob.scale ~= 1
	ob.scale = ob.scale .^ sup;
end


text:
function imdb = setupCaltech256(datasetDir, varargin)
% SETUPCALTECH256    Setup Caltech 256 and 101 datasets
%    This is similar to SETUPGENERIC(), with modifications to setup
%    Caltech-101 and Caltech-256 according to the standard
%    evaluation protocols. Specific options include:
%
%    Variant:: 'caltech256'
%      Either 'caltech101' or 'caltech256'.
%
%    AutoDownload:: true
%      Automatically download the data from the Internet if not
%      found at DATASETDIR.
%
%    See:: SETUPGENERIC().

% Author: Andrea Vedaldi

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

opts.lite = false ;
opts.numTrain = 30 ;
opts.numTest = 25 ;
opts.seed = 1 ;
opts.variant = 'caltech256' ;
opts.autoDownload = true ;
opts = vl_argparse(opts, varargin) ;

% Download and unpack
vl_xmkdir(datasetDir) ;
switch opts.variant
  case 'caltech101'
    name = '101_ObjectCategories' ;
    checkClassName = 'airplanes' ;
    url = 'http://www.vision.caltech.edu/Image_Datasets/Caltech101/101_ObjectCategories.tar.gz' ;
    numClasses = 102 ;
  case 'caltech256'
    name = '256_ObjectCategories' ;
    checkClassName = '001.ak47' ;
    url = 'http://www.vision.caltech.edu/Image_Datasets/Caltech256/256_ObjectCategories.tar' ;
    numClasses = 257 ;
  otherwise
    error('Uknown dataset variant ''%s''.', opts.variant) ;
end

if exist(fullfile(datasetDir, checkClassName), 'file')
  % ok
elseif exist(fullfile(datasetDir, name, checkClassName), 'file')
  datasetDir = fullfile(datasetDir, name) ;
elseif opts.autoDownload
  fprintf('Downloading %s data to ''%s''. This will take a while.', opts.variant, datasetDir) ;
  untar(url, datasetDir) ;
  datasetDir = fullfile(datasetDir, name) ;
else
  error('Could not find %s dataset in ''%s''', opts.variant, datasetDir) ;
end

% Read classes
imdb = setupGeneric(datasetDir, ...
  'numTrain', opts.numTrain, 'numVal', 0, 'numTest', opts.numTest,  ...
  'expectedNumClasses', numClasses, ...
  'seed', opts.seed, 'lite', opts.lite) ;

% In Caltech 256 background is not evaluated
switch opts.variant
  case 'caltech101'
  case 'caltech256'
    imdb.images.sets(classes == 257) = 0 ;
    ok = find(imdb.images.sets ~= 0) ;
    imdb.images.id = ids(ok) ;
    imdb.images.name = names(ok) ;
    imdb.images.set = sets(ok) ;
    imdb.images.class = classes(ok) ;
end

author: jianxiongxiao

repo: ProfXkit


text:
function [ r, z, rho, c, s ] = dchud ( r, ldr, p, x, z, ldz, nz, y, rho )

%*****************************************************************************80
%
%% DCHUD updates an augmented Cholesky decomposition.
%
%  Discussion:
%
%    DCHUD can also update the triangular part of an augmented QR
%    decomposition.
%
%    Specifically, given an upper triangular matrix R of order P, a row vector
%    X, a column vector Z, and a scalar Y, DCHUD determines a unitary matrix
%    U and a scalar ZETA such that
%
%           (R  Z)     (RR   ZZ )
%      U  * (    )  =  (        ),
%           (X  Y)     ( 0  ZETA)
%
%    where RR is upper triangular.
%
%    If R and Z have been obtained from the factorization of a least squares
%    problem, then RR and ZZ are the factors corresponding to the problem
%    with the observation (X,Y) appended.  In this case, if RHO is the
%    norm of the residual vector, then the norm of the residual vector of
%    the updated problem is sqrt ( RHO * RHO + ZETA * ZETA ).  DCHUD will
%    simultaneously update several triplets (Z, Y, RHO).
%
%    For a less terse description of what DCHUD does and how
%    it may be applied, see the LINPACK guide.
%
%    The matrix U is determined as the product U(P)*...*U(1),
%    where U(I) is a rotation in the (I,P+1) plane of the form
%
%      (     C(I)      S(I) )
%      (                    ).
%      (    -S(I)      C(I) )
%
%    The rotations are chosen so that C(I) is real.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real R(LDR,P), the upper triangular matrix to be
%    updated.  The part of R below the diagonal is not referenced.
%
%    Input, integer LDR, the leading dimension of the array R.
%    LDR must be at least equal to P.
%
%    Input, integer P, the order of the matrix R.
%
%    Input, real X(P), the row to be added to R.
%
%    Input, real Z(LDZ,NZ), contains NZ P-vectors to be updated with R.
%
%    Input, integer LDZ, the leading dimension of the array Z.
%    LDZ must be at least P.
%
%    Input, integer NZ, the number of vectors to be updated.  NZ may be
%    zero, in which case Z, Y, and RHO are not referenced.
%
%    Input, real Y(NZ), the scalars for updating the vectors Z.
%
%    Input, real RHO(NZ), the norms of the residual vectors to be updated.
%    If RHO(J) is negative, it is left unaltered.
%
%    Output, real R(LDR,P), the updated matrix.
%
%    Output, real Z(LDZ,NZ), the updated vectors.
%
%    Output, real RHO(NZ), the updated norms of the residual vectors.
%
%    Output, real C(P), S(P), the cosines and sines of the
%    transforming rotations.
%

%
%  Update R.
%
  for j = 1: p

    xj = x(j);
%
%  Apply the previous rotations.
%
    for i = 1 : j-1
      t = c(i) * r(i,j) + s(i) * xj;
      xj = c(i) * xj - s(i) * r(i,j);
      r(i,j) = t;
    end
%
%  Compute the next rotation.
%
    [ c(j), s(j), r(j,j), xj ] = drotg ( r(j,j), xj );

  end
%
%  If required, update Z and RHO.
%
  for j = 1 : nz

    zeta = y(j);

    for i = 1 : p
      t =    c(i) * z(i,j) + s(i) * zeta;
      zeta = c(i) * zeta   - s(i) * z(i,j);
      z(i,j) = t;
    end

    azeta = abs ( zeta );

    if ( azeta ~= 0.0 & 0.0 <= rho(j) )
      scale = azeta + rho(j);
      rho(j) = scale * sqrt ( ( azeta / scale )^2 + ( rho(j) / scale )^2 );
    end

  end

  return
end


text:
classdef pivprogress < handle

	%   Inspired by the parfor_progress script made by Daniel Terry:
	%   https://www.mathworks.com/matlabcentral/fileexchange/53773-parfor-progressbar

	% Public properties
	properties (SetAccess=protected, GetAccess=public)
		hdl;      % Waitbar figure object handle
		N;        % Total number of iterations expected before completion
	end

	properties (Dependent, GetAccess=public)
		percent;  % Percentage of completed iterations
	end

	% Internal properties
	properties (SetAccess=protected, GetAccess=protected, Hidden)
		ipcfile;  % Path to temporary file for inter-process communication
		htimer;   % Timer object that checks ipcfile for completed iterations
	end

	methods
		% Create a new progress bar with N_init iterations before completion.
		function this = pivprogress(N_init, handle)
			% Create a unique inter-process communication file.
			for i=1:10
				f = sprintf('%s%d.txt', mfilename, round(rand*1000));
				this.ipcfile = fullfile(tempdir, f);
				if ~exist(this.ipcfile,'file'), break; end
			end
			if exist(this.ipcfile,'file')
				error('Too many temporary files. Clear out tempdir.');
			end
			try
				setappdata(0,'par_start_time',tic)
				setappdata(0,'par_percent',0)
			catch
			end

			% Create a new waitbar
			this.N = N_init;
			if nargin == 2
				this.hdl = handle;% waitbar(0, varargin{:});
			else
				this.hdl = [];
			end
			% Create timer to periodically update the waitbar in the GUI thread.
			this.htimer = timer( 'ExecutionMode','fixedSpacing', 'Period',2, ...
				'BusyMode','drop', 'Name',mfilename, ...
				'TimerFcn',@(x,y)this.tupdate);
			start(this.htimer);
		end

		function delete(this)
			this.close();
		end

		function close(this)
			% Closer the progress bar and clean up internal state.
			% Stop the timer
			if isa(this.htimer,'timer') && isvalid(this.htimer)
				stop(this.htimer);
				pause(0.01);
				delete(this.htimer);
			end
			this.htimer = [];
			% Delete the IPC file.
			if exist(this.ipcfile,'file')
				delete(this.ipcfile);
			end
			this.hdl = [];
		end

		function percent = get.percent(this)
			% Calculate the fraction of completed iterations from IPC file.
			if ~exist(this.ipcfile, 'file')
				percent = 0;  % File may not exist before the first iteration
			else
				fid = fopen( this.ipcfile, 'r' );
				percent = sum(fscanf(fid, '%d')) / this.N;
				percent = max(0, min(1,percent) );
				fclose(fid);
			end
		end

		function iterate(this, Nitr)
			% Update the progress bar by Nitr iterations (or 1 if not specified).
			if nargin<2,  Nitr = 1;  end
			fid = fopen(this.ipcfile, 'a');
			fprintf(fid, '%d\n', Nitr);
			fclose(fid);
		end
	end %public methods

	methods (Access=protected, Hidden)
		function tupdate(this)
			% Check the IPC file and update the waitbar with progress.
			if ishandle(this.hdl)
				try
					par_start_time=getappdata(0,'par_start_time');
					last_par_percent=getappdata(0,'par_percent');
					elapsd = toc(par_start_time);
					if this.percent > 0.01
						remain_t = elapsd / this.percent *(1-this.percent);
					else
						remain_t=inf;
					end
					hrs=remain_t/60^2;
					mins=(hrs-floor(hrs))*60;
					secs=(mins-floor(mins))*60;
					hrs=floor(hrs);
					mins=floor(mins);
					secs=floor(secs);

					if elapsd < 5
						remain_string='Remaining time: calculating...';
					else
						remain_string=['Remaining time: ' sprintf('%2.2d', hrs) 'h ' sprintf('%2.2d', mins) 'm ' sprintf('%2.2d', secs) 's'];
					end
					if this.percent ~= last_par_percent
						if ~isinf(remain_t)
							set(this.hdl, 'string' , ['Total progress: ' num2str(floor(100*this.percent)) '%' sprintf('\n') remain_string ]); %#ok<SPRINTFN>
						end
						setappdata(0,'par_percent',this.percent)
					end
				catch
				end
				drawnow;
			elseif isempty(this.hdl)
				clc
				disp([num2str(floor(100*this.percent)) ' %']);
			else
				% Kill the timer if the waitbar is closed.
				close(this);
			end
			if exist('cancel_piv','file')
				close(this)
			end
		end
	end %private methods

end %classdef

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function ort_vec = get_ort_vec(lew_rgb,mask_local)
% compute the mean color w.r.t pixels covered by the binary saliency mask
    [ii,jj] = find(mask_local~=0);
    ort_vec = zeros(1,3); % RGB color vector
    for k = 1 : length(ii)
        ort_vec(1,1) = ort_vec(1,1) + lew_rgb(ii(k),jj(k),1);
        ort_vec(1,2) = ort_vec(1,2) + lew_rgb(ii(k),jj(k),2);
        ort_vec(1,3) = ort_vec(1,3) + lew_rgb(ii(k),jj(k),3);
    end
    ort_vec = ort_vec / length(ii);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
p1 = [3 3 0]';
p2 = [3 3 -2]';
axis([-5 5 -5 5]);
plot

author: petercorke

repo: machinevision-toolbox-matlab


text:
%% Spherical Functions
%
%%
% By a variable of type @S2Fun it is possible to represent an entire
% function on the two dimensional sphere. A typical example of such a
% function is the pole density function of a given ODF with respect to a
% fixed crystal direction.

% the famouse Santa Fe orientation distribution function
odf = SantaFe;

% the (100) pole density function
pdf = odf.calcPDF(Miller(1,0,0,odf.CS))

%%
% Since, the variable |pdf| stores all information about this function we
% may evaluate it for any direction |r|

% take a random direction
r = vector3d.rand;

% and evaluate the pdf at this direction
pdf.eval(r)

%%
% We may also plot the function in any spherical projection

plot(pdf)

%%
% or find its local maxima

[~,localMax] = max(pdf,'numLocal',12)

annotate(localMax)

%%
% A complete list of operations that can be performed with spherical
% functions can be found in section <S2FunOperations.html Operations>.
%

%% Representation of Spherical Functions
%
% In MTEX there exist different ways for representing spherical functions
% internally.
%
% || harmonic expansion || @S2FunHarmonic ||
% || finite elements || @S2FunTri ||
% || function handle || @S2FunHandle ||
% || Bingham distribution || @BinghamS2 ||
%
% All representations allow for the same operations which are specified for
% the abstact class @S2Fun. In particular it is possible
% to calculate with spherical functions as with ordinary numbers, i.e., you
% can add, multiply arbitrary functions, take the mean, integrate them or
% compute gradients.
%
%% Generalizations of Spherical Functions
%
% || spherical vector fields || @S2VectorField ||
% || spherical axis fields || @S2AxisField ||
% || radial spherical functions || @S2Kernel ||
% || symmetric spherical functions || @S2FunHarmonicSym ||
%

author: mtex-toolbox

repo: mtex

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


text:
function g =  getGroupOverlap(row, col)

    N = row*col;

    g = sparse(zeros(N,1));
    g = diag(g);

    %% build overlapping group

    % top let corner group

    g(1, 1 ) = 1;
    g(1, 2 ) = 1;
    g(1, 1+row ) = 1;
    g(1, 2+row ) = 1;

     % bottom let corner group

    g(row, row-1 ) = 1;
    g(row, row ) = 1;
    g(row, row-1+row ) = 1;
    g(row, row+row ) = 1;

    % top right corner group

    g((col-1)*row+1, (col-1)*row+1 ) = 1;
    g((col-1)*row+1, (col-2)*row+1 ) = 1;
    g((col-1)*row+1, (col-1)*row+2 ) = 1;
    g((col-1)*row+1, (col-2)*row+2 ) = 1;

    % bottom right corner group

    g((col-1)*row+row , (col-1)*row+row-1 ) = 1;
    g((col-1)*row+row , (col-2)*row+row-1 ) = 1;
    g((col-1)*row+row , (col-1)*row+row ) = 1;
    g((col-1)*row+row , (col-2)*row+row ) = 1;

    % boundary group

    for i=2:col-1
        % top row rgoup

        j = 1;

        g( (i-1)*row+j  , (i-2)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-2)*row+j+1 ) = 1;

        g( (i-1)*row+j  , (i-1)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j+1 ) = 1;

        g( (i-1)*row+j  , i*row+j   ) = 1;
        g( (i-1)*row+j  , i*row+j+1 ) = 1;

        % bottom row group

        j = row;

        g( (i-1)*row+j  , (i-2)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-2)*row+j   ) = 1;

        g( (i-1)*row+j  , (i-1)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j   ) = 1;

        g( (i-1)*row+j  , i*row+j-1 ) = 1;
        g( (i-1)*row+j  , i*row+j   ) = 1;

    end

    for j=2:row-1
          % left column rgoup

        i=1;

        g( (i-1)*row+j  , (i-1)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j+1 ) = 1;
        g( (i-1)*row+j  , i*row+j-1 ) = 1;
        g( (i-1)*row+j  , i*row+j   ) = 1;
        g( (i-1)*row+j  , i*row+j+1 ) = 1;

        % right  column group

        i=col;

        g( (i-1)*row+j  , (i-2)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-2)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-2)*row+j+1 ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j+1 ) = 1;

    end

    for i=2:col-1
        for j=2:row-1

            g( (i-1)*row+j  , (i-2)*row+j-1 ) = 1;
            g( (i-1)*row+j  , (i-2)*row+j   ) = 1;
            g( (i-1)*row+j  , (i-2)*row+j+1 ) = 1;
            g( (i-1)*row+j  , (i-1)*row+j-1 ) = 1;
            g( (i-1)*row+j  , (i-1)*row+j   ) = 1;
            g( (i-1)*row+j  , (i-1)*row+j+1 ) = 1;
            g( (i-1)*row+j  , i*row+j-1 ) = 1;
            g( (i-1)*row+j  , i*row+j   ) = 1;
            g( (i-1)*row+j  , i*row+j+1 ) = 1;

        end
    end

end

author: andrewssobral

repo: lrslibrary


text:
function sys = sdpvar(varargin)
%SDPVAR Create symbolic decision variable
%
%   You can create a sdpvar variable by:
%     X = SDPVAR(n)               Symmetric nxn matrix
%     X = SDPVAR(n,n)             Symmetric nxn matrix
%     X = SDPVAR(n,m)             Full nxm matrix (n~=m)
%
%   Definition of multiple scalars can be simplified
%     SDPVAR x y z w
%
%   The parametrizations supported are
%     X = SDPVAR(n,n,'full')      Full nxn matrix
%     X = SDPVAR(n,n,'symmetric') Symmetric nxn matrix
%     X = SDPVAR(n,n,'toeplitz')  Symmetric Toeplitz
%     X = SDPVAR(n,n,'hankel')    Symmetric Hankel
%     X = SDPVAR(n,n,'skew')      Skew-symmetric
%
%   The letters 'sy','f','ha', 't' and 'sk' are searched for in the third argument
%   hence sdpvar(n,n,'toeplitz') gives the same result as sdpvar(n,n,'t')
%
%   Only square Toeplitz and Hankel matries are supported
%
%   A scalar is defined as a 1x1 matrix
%
%   Higher-dimensional matrices are also supported, although this currently
%   is an experimental feature with limited use. The type flag applies to
%   the lowest level slice.
%
%     X = SDPVAR(n,n,n,'full')      Full nxnxn matrix
%
%   In addition to the matrix type, a fourth argument
%   can be used to obtain a complex matrix. All the
%   matrix types above apply to a complex matrix, and
%   in addition a Hermitian type is added
%
%     X = SDPVAR(n,n,'hermitian','complex') Complex Hermitian nxn matrix (X=X'=conj(X.'))
%
%   The other types are obtained as above
%     X = SDPVAR(n,n,'symmetric','complex') Complex symmetric nxn matrix (X=X.')
%     X = SDPVAR(n,n,'full','complex')      Complex full nxn matrix
%     ... and the same for Toeplitz, Hankel and skew-symmetric
%
%   See also @SDPVAR/SET, INTVAR, BINVAR, methods('sdpvar'), SEE

superiorto('sdpvar');
if nargin==0
    return
end

if isstruct(varargin{1})
    sys = class(varargin{1},'ncvar');
    return
end

% To speed up dualization, we keep track of primal SDP cones
% [0 0] :  Nothing known (cleared in some operator, or none-cone to start with)
% [1 0] :  Primal cone
% [1 1] :  Primal cone + DOUBLE
% [1 2 x] :  Primal cone + SDPVAR
% [-1 1] : -Primal cone + DOUBLE
% [-1 2 x] : -Primal cone + SDPVAR

conicinfo = [0 0];

if ischar(varargin{1})
    switch varargin{1}
        case 'clear'
            disp('Obsolete comand');
            return
        case 'nvars'
            sys = yalmip('nvars');%THIS IS OBSAOLETE AND SHOULD NOT BE USED
            return
        otherwise
            n = length(varargin);
            varnames = varargin;
            for k = 1:n
                varcmd{k}='(1,1)';
                lp=strfind(varargin{k},'(');
                rp=strfind(varargin{k},')');
                if isempty(lp) & isempty(rp)
                    if ~isvarname(varargin{k})
                        error('Not a valid variable name.')
                    end
                else
                    if (~isempty(lp))&(~isempty(rp))
                        if min(lp)<max(rp)
                            varnames{k} = varargin{k}(1:lp-1);
                            varcmd{k}=varargin{k}(lp:rp);
                        else
                            error('Not a valid variable name.')
                        end
                    else
                        error('Not a valid variable name.')
                    end
                end
            end
            for k = 1:n
                if isequal(varnames{k},'i') | isequal(varnames{k},'j')
                    if length(dbstack) == 1
                        assignin('caller',varnames{k},eval(['sdpvar' varcmd{k}]));
                    else
                        error(['Due to a bug in MATLAB, use ' varnames{k} ' = sdpvar' varcmd{k} ' instead.']);
                    end
                else
                    assignin('caller',varnames{k},eval(['ncvar' varcmd{k}]));
                end
            end
            return
    end
end

% *************************************************************************
% Maybe new NDSDPVAR syntax
% *************************************************************************
if nargin > 2 & isa(varargin{3},'double') & ~isempty(varargin{3})
    sys = ndsdpvar(varargin{:});
    return
end

% Supported matrix types
% - symm
% - full
% - skew
% - hank
% - toep
switch nargin
    case 1 %Bug in MATLAB 5.3!! sdpvar called from horzcat!!!!????
        if isempty(varargin{1})
            sys = varargin{1};
            return
        end
        if isa(varargin{1},'sdpvar')
            sys = varargin{1};
            sys.typeflag = 0;
            return
        end
        n = varargin{1};
        m = varargin{1};
        if sum(n.*m)==0
            sys = zeros(n,m);
            return
        end
        if (n==m)
            matrix_type = 'symm';
            nvar = sum(n.*(n+1)/2);
            conicinfo = [1 0];
        else
            matrix_type = 'full';
            nvar = sum(n.*m);
        end
    case 2
        n = varargin{1};
        m = varargin{2};
        if length(n)~=length(m)
            error('The dimensions must have the same lengths')
        end
        if sum(n.*m)==0
            sys = zeros(n,m);
            return
        end
        if (n==m)
            matrix_type = 'symm';
            nvar = sum(n.*(n+1)/2);
            conicinfo = [1 0];
        else
            matrix_type = 'full';
            nvar = sum(n.*m);
        end
    case {3,4}
        n = varargin{1};
        m = varargin{2};
        if sum(n.*m)==0
            sys = zeros(n,m);
            return
        end

        % Check for complex or real
        if (nargin == 4)
            if isempty(varargin{4})
                varargin{4} = 'real';
            else
                if ~ischar(varargin{4})
                    help sdpvar
                    error('Fourth argument should be ''complex'' or ''real''')
                end
            end
            index_cmrl = strmatch_octavesafe(varargin{4},{'real','complex'});
            if isempty(index_cmrl)
                error('Fourth argument should be ''complex'' or ''real''. See help above')
            end
        else
            if ~ischar(varargin{3})
                help sdpvar
                error('Third argument should be ''symmetric'', ''full'', ''hermitian'',...See help above')
            end
            index_cmrl = 1;
        end;

        if isempty(varargin{3})
            if n==m
                index_type = 7; %Default symmetric
            else
                index_type = 4;
            end
        else
            if ~isempty(strmatch_octavesafe(varargin{3},{'complex','real'}))
                % User had third argument as complex or real
                error(['Third argument should be ''symmetric'', ''full'', ''toeplitz''... Maybe you meant sdpvar(n,n,''full'',''' varargin{3} ''')'])
            end
            index_type = strmatch_octavesafe(varargin{3},{'toeplitz','hankel','symmetric','full','rhankel','skew','hermitian'});
        end

        if isempty(index_type)
            error(['Matrix type "' varargin{3} '" not supported'])
        else
            switch index_type+100*(index_cmrl-1)
                case 1
                    if n~=m
                        error('Toeplitz matrix must be square')
                    else
                        matrix_type = 'toep';
                        nvar = n;
                    end

                case 2
                    if n~=m
                        error('Hankel matrix must be square')
                    else
                        matrix_type = 'hank';
                        nvar = n;
                    end

                case 3
                    if n~=m
                        error('Symmetric matrix must be square')
                    else
                        matrix_type = 'symm';
                        nvar = sum(n.*(n+1)/2);
                    end

                case 4
                    matrix_type = 'full';
                    nvar = sum(n.*m);
                    if nvar==1
                        matrix_type = 'symm';
                    end

                case 5
                    if n~=m
                        error('Hankel matrix must be square')
                    else
                        matrix_type = 'rhankel';
                        nvar = 2*n-1;
                    end

                case 6
                    if n~=m
                        error('Skew symmetric matrix must be square')
                    else
                        matrix_type = 'skew';
                        nvar = (n*(n+1)/2)-n;
                    end

                case 7
                    if n~=m
                        error('Symmetric matrix must be square')
                    else
                        matrix_type = 'symm';
                        nvar = n*(n+1)/2;
                    end

                case 101
                    if n~=m
                        error('Toeplitz matrix must be square')
                    else
                        matrix_type = 'toep complex';
                        nvar = 2*n;
                    end

                case 102
                    if n~=m
                        error('Hankel matrix must be square')
                    else
                        matrix_type = 'hank complex';
                        nvar = (2*n);
                    end

                case 103
                    if n~=m
                        error('Symmetric matrix must be square')
                    else
                        matrix_type = 'symm complex';
                        nvar = 2*n*(n+1)/2;
                    end

                case 104
                    matrix_type = 'full complex';
                    nvar = 2*n*m;
                    if nvar==1
                        matrix_type = 'symm complex';
                    end

                case 105
                    if n~=m
                        error('Hankel matrix must be square')
                    else
                        matrix_type = 'rhankel complex';
                        nvar = 2*(2*n-1);
                    end

                case 106
                    if n~=m
                        error('Skew symmetric matrix must be square')
                    else
                        matrix_type = 'skew complex';
                        nvar = 2*((n*(n+1)/2)-n);
                    end

                case 107
                    if n~=m
                        error('Hermitian matrix must be square')
                    else
                        matrix_type = 'herm complex';
                        nvar = n*(n+1)/2+(n*(n+1)/2-n);
                    end

                otherwise
                    error('Bug! Report!');
            end

        end

    case 5 % Fast version for internal use
        sys.basis = varargin{5};
        sys.lmi_variables=varargin{4};
        sys.dim(1) = varargin{1};
        sys.dim(2) = varargin{2};
        sys.typeflag = 0;
        sys.savedata = [];
        sys.extra = [];
        sys.extra.expanded = [];
        sys.extra.createTime = definecreationtime;
        sys.originalbasis = [];
        sys.leftfactors{1} = [];
        sys.rightfactors{1} = [];
        sys.midfactors{1} = [];
        sys.conicinfo = 0;

        % Find zero-variables
        constants = find(sys.lmi_variables==0);
        if ~isempty(constants);
            sys.lmi_variables(constants)=[];
            sys.basis(:,1) = sys.basis(:,1) + sum(sys.basis(:,1+constants),2);
            sys.basis(:,1+constants)=[];
        end
        if isempty(sys.lmi_variables)
            sys = full(reshape(sys.basis(:,1),sys.dim(1),sys.dim(2)));
        else
            sys = class(sys,'sdpvar');
        end
        return
    case 6 % Fast version for internal use
        sys.basis = varargin{5};
        sys.lmi_variables=varargin{4};
        sys.dim(1) = varargin{1};
        sys.dim(2) = varargin{2};
        sys.typeflag = varargin{6};
        sys.savedata = [];
        sys.extra = [];
        sys.extra.expanded = [];
        sys.extra.createTime = definecreationtime;
        sys.conicinfo = 0;
        sys.originalbasis = [];
        sys.leftfactors{1} = [];
        sys.rightfactors{1} = [];
        sys.midfactors{1} = [];
        % Find zero-variables
        constants = find(sys.lmi_variables==0);
        if ~isempty(constants);
            sys.lmi_variables(constants)=[];
            sys.basis(:,1) = sys.basis(:,1) + sum(sys.basis(:,1+constants),2);
            sys.basis(:,1+constants)=[];
        end
        if isempty(sys.lmi_variables)
            sys = full(reshape(sys.basis(:,1),sys.dim(1),sys.dim(2)));
        else
            sys = class(sys,'ncvar');
        end
        return
    case 7 % Fast version for internal use
        sys.basis = varargin{5};
        sys.lmi_variables=varargin{4};
        sys.dim(1) = varargin{1};
        sys.dim(2) = varargin{2};
        sys.typeflag = varargin{6};
        sys.savedata = [];
        sys.extra = varargin{7};
        sys.extra.expanded = [];
        sys.extra.createTime = definecreationtime;
        sys.conicinfo = varargin{7};
        % Find zero-variables
        constants = find(sys.lmi_variables==0);
        if ~isempty(constants);
            sys.lmi_variables(constants)=[];
            sys.basis(:,1) = sys.basis(:,1) + sum(sys.basis(:,1+constants),2);
            sys.basis(:,1+constants)=[];
        end
        if isempty(sys.lmi_variables)
            sys = full(reshape(sys.basis(:,1),sys.dim(1),sys.dim(2)));
        else
            sys = class(sys,'sdpvar');
        end
        return

    otherwise
        error('Wrong number of arguments in sdpvar creation');
end

if isempty(n) | isempty(m)
    error('Size must be integer valued')
end;
if ~((abs((n-ceil(n)))+ abs((m-ceil(m))))==0)
    error('Size must be integer valued')
end

nonCommutingTable = yalmip('nonCommutingTable');
[monomtable,variabletype] = yalmip('monomtable');
lmi_variables = (1:nvar)+max(size(nonCommutingTable,1),size(monomtable,1));

for blk = 1:length(n)
    switch matrix_type

        case 'full'
            basis{blk} = [spalloc(n(blk)*m(blk),1,0) speye(n(blk)*m(blk))];%speye(nvar)];

        case 'full complex'
            basis = [spalloc(n*m,1,0) speye(nvar/2) speye(nvar/2)*sqrt(-1)];

        case 'symm'
            if 0
                basis = spalloc(n^2,1+nvar,n^2);
                l = 2;
                an_empty = spalloc(n,n,2);
                for i=1:n
                    temp = an_empty;
                    temp(i,i)=1;
                    basis(:,l)=temp(:);
                    l = l+1;
                    for j=i+1:n,
                        temp = an_empty;
                        temp(i,j)=1;
                        temp(j,i)=1;
                        basis(:,l)=temp(:);
                        l = l+1;
                    end
                end
            else
                % Hrm...fast but completely f*d up

                Y = reshape(1:n(blk)^2,n(blk),n(blk));
                Y = tril(Y);
                Y = (Y+Y')-diag(sparse(diag(Y)));
                [uu,oo,pp] = unique(Y(:));
                if 1
                    basis{blk} = sparse(1:n(blk)^2,pp+1,1);
                else
                    basis{blk} = lazybasis(n^2,1+(n*(n+1)/2),1:n(blk)^2,pp+1,ones(n(blk)^2,1));
                end

            end

        case 'symm complex'
            basis = spalloc(n^2,1+nvar,2);
            l = 2;
            an_empty = spalloc(n,n,2);
            for i=1:n
                temp = an_empty;
                temp(i,i)=1;
                basis(:,l)=temp(:);
                l = l+1;
                for j=i+1:n,
                    temp = an_empty;
                    temp(i,j)=1;
                    temp(j,i)=1;
                    basis(:,l)=temp(:);
                    l = l+1;
                end
            end
            for i=1:n
                temp = an_empty;
                temp(i,i)=sqrt(-1);
                basis(:,l)=temp(:);
                l = l+1;
                for j=i+1:n,
                    temp = an_empty;
                    temp(i,j)=sqrt(-1);
                    temp(j,i)=sqrt(-1);
                    basis(:,l)=temp(:);
                    l = l+1;
                end
            end

        case 'herm complex'
            basis = spalloc(n^2,1+nvar,2);
            l = 2;
            an_empty = spalloc(n,n,2);
            for i=1:n
                temp = an_empty;
                temp(i,i)=1;
                basis(:,l)=temp(:);
                l = l+1;
                for j=i+1:n,
                    temp = an_empty;
                    temp(i,j)=1;
                    temp(j,i)=1;
                    basis(:,l)=temp(:);
                    l = l+1;
                end
            end
            for i=1:n
                for j=i+1:n,
                    temp = an_empty;
                    temp(i,j)=sqrt(-1);
                    temp(j,i)=-sqrt(-1);
                    basis(:,l)=temp(:);
                    l = l+1;
                end
            end

        case 'skew'
            basis = spalloc(n^2,1+nvar,2);
            l = 2;
            an_empty = spalloc(n,n,2);
            for i=1:n
                for j=i+1:n,
                    temp = an_empty;
                    temp(i,j)=1;
                    temp(j,i)=-1;
                    basis(:,l)=temp(:);
                    l = l+1;
                end
            end

        case 'skew complex'
            basis = spalloc(n^2,1+nvar,2);
            l = 2;
            an_empty = spalloc(n,n,2);
            for i=1:n
                for j=i+1:n,
                    temp = an_empty;
                    temp(i,j)=1;
                    temp(j,i)=-1;
                    basis(:,l)=temp(:);
                    l = l+1;
                end
            end
            for i=1:n
                for j=i+1:n,
                    temp = an_empty;
                    temp(i,j)=sqrt(-1);
                    temp(j,i)=-sqrt(-1);
                    basis(:,l)=temp(:);
                    l = l+1;
                end
            end

        case 'toep'
            basis = spalloc(n^2,1+nvar,2);
            an_empty = spalloc(n,1,1);
            for i=1:n,
                v = an_empty;
                v(i)=1;
                temp = sparse(toeplitz(v));
                basis(:,i+1) = temp(:);
            end

            % Notice, complex Toeplitz not Hermitian
        case 'toep complex'
            basis = spalloc(n^2,1+nvar,2);
            an_empty = spalloc(n,1,1);
            for i=1:n,
                v = an_empty;
                v(i)=1;
                temp = sparse(toeplitz(v));
                basis(:,i+1) = temp(:);
            end
            for i=1:n,
                v = an_empty;
                v(i)=sqrt(-1);
                temp = sparse(toeplitz(v));
                basis(:,n+i+1) = temp(:);
            end

        case 'hank'
            basis = spalloc(n^2,1+nvar,2);
            an_empty = spalloc(n,1,1);
            for i=1:n,
                v = an_empty;
                v(i)=1;
                temp = sparse(hankel(v));
                basis(:,i+1) = temp(:);
            end

        case 'hank complex'
            basis = spalloc(n^2,1+nvar,2);
            an_empty = spalloc(n,1,1);
            for i=1:n,
                v = an_empty;
                v(i)=1;
                temp = sparse(hankel(v));
                basis(:,i+1) = temp(:);
            end
            for i=1:n,
                v = an_empty;
                v(i)=sqrt(-1);
                temp = sparse(hankel(v));
                basis(:,n+i+1) = temp(:);
            end

        case 'rhankel'
            basis = spalloc(n^2,1+nvar,2);
            an_empty = spalloc(2*n-1,1,1);
            for i=1:nvar,
                v = an_empty;
                v(i)=1;
                temp = sparse(hankel(v(1:n),[v(n);v(n+1:2*n-1)]));
                basis(:,i+1) = temp(:);
            end

        case 'rhankel complex'
            basis = spalloc(n^2,1+nvar,2);
            an_empty = spalloc(2*n-1,1,1);
            for i=1:nvar/2,
                v = an_empty;
                v(i)=1;
                temp = sparse(hankel(v(1:n),[v(n);v(n+1:2*n-1)]));
                basis(:,i+1) = temp(:);
            end
            for i=1:nvar/2,
                v = an_empty;
                v(i)=sqrt(-1);
                temp = sparse(hankel(v(1:n),[v(n);v(n+1:2*n-1)]));
                basis(:,nvar/2+i+1) = temp(:);
            end

        otherwise
            error('Bug! Report')
    end

end

% Update noncommuting table and monomtables
nonCommutingTable(lmi_variables,1) = nan;
nonCommutingTable(lmi_variables,2) = lmi_variables;
yalmip('nonCommutingTable',nonCommutingTable);
variabletype(lmi_variables) = 0;
monomtable(lmi_variables(end),lmi_variables(end)) = 0;
yalmip('setmonomtable',monomtable,variabletype);
%if strcmp(matrix_type,'NonHermitian')
%    yalmip('setNonHermitianNonCommuting',lmi_variables);
%end

% Create an object
if isa(basis,'cell')
    top = 1;
    for blk = 1:length(n)
        sys{blk}.basis=basis{blk};
        nn = size(sys{blk}.basis,2)-1;
        sys{blk}.lmi_variables = lmi_variables(top:top+nn-1);
        top = top + nn;
        sys{blk}.dim(1) = n(blk);
        sys{blk}.dim(2) = m(blk);
        sys{blk}.typeflag = 0;
        sys{blk}.savedata = [];
        sys{blk}.extra = [];
        sys{blk}.extra.expanded = [];
        sys{blk}.extra.createTime = definecreationtime;
        sys{blk}.conicinfo = conicinfo;
        sys{blk}.originalbasis = [];
        sys{blk}.leftfactors{1} = [];
        sys{blk}.rightfactors{1} = [];
        sys{blk}.midfactors{1} = [];
        sys{blk} = class(sys{blk},'ncvar');
    end
    if length(n)==1
        sys = sys{1};
    end
else
    sys.basis=basis;
    sys.lmi_variables = lmi_variables;
    sys.dim(1) = n;
    sys.dim(2) = m;
    sys.typeflag = 0;
    sys.savedata = [];
    sys.extra = [];
    sys.extra.expanded = [];
    sys.extra.createTime = definecreationtime;
    sys.conicinfo = conicinfo;
    sys.originalbasis = [];
    sys.leftfactors{1} = [];
    sys.rightfactors{1} = [];
    sys.midfactors{1} = [];
    sys = class(sys,'ncvar');
    if ~isreal(basis)
        % Add internal information about complex pairs
        complex_elements = find(any(imag(basis),2));
        complex_pairs = [];
        for i = 1:length(complex_elements)
            complex_pairs = [complex_pairs;lmi_variables(find(basis(complex_elements(i),:))-1)];
        end
        complex_pairs = uniquesafe(complex_pairs,'rows');
        yalmip('addcomplexpair',complex_pairs);
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Compute the frequency range of the mel filterbank
% Inputs:
%   start_freq:     the lowerest linear frequency included in the calculation
%   linear_samp:    sampling frequency
%   N_mel:          number of mel banks used
% Output:
%   bin_lower_edge:    the lower edge linear frequency of each mel bank
%   bin_upper_edge:    the upper edge linear frequency of each mel bank

function [bin_lower_edge, bin_upper_edge] = mel_bank_range(start_freq, linear_samp, N_mel)

% get the mel version of start frequency and linear sampling rate
start_mel = linear2mel(start_freq);
samp_mel = linear2mel(linear_samp/2);

% slow implementation
if 0
    for i=1:N_mel
        % find the lower edge linear frequency of the bank
        tmp_mel = (i-1)*(samp_mel-start_mel) / (N_mel+1);
        bin_lower_edge(i) = mel2linear(tmp_mel+start_mel);
        % find the upper edge linear frequency of the bank
        tmp_mel = (i+1)*(samp_mel-start_mel) / (N_mel+1);
        bin_upper_edge(i) =  mel2linear(tmp_mel+start_mel);
    end
else
    % fast implementation
    i=1:N_mel;
    % find the lower edge linear frequency of the bank
    tmp_mel = (i-1)*(samp_mel-start_mel) / (N_mel+1);
    bin_lower_edge = mel2linear(tmp_mel+start_mel);
    % find the upper edge linear frequency of the bank
    tmp_mel = (i+1)*(samp_mel-start_mel) / (N_mel+1);
    bin_upper_edge =  mel2linear(tmp_mel+start_mel);
end

author: singaxiong

repo: SignalGraph


text:
function [formula, nH, charge] = getFormulaAndChargeFromInChI(inchi)

% [formula,charge] = getFormulaAndChargeFromInChI(inchi)
%
% INPUT
% inchi.......Nonstandard IUPAC InChI for a particular pseudoisomer of a
%             metabolite
%
% OUTPUTS
% formula....The chemical formula for the input pseudoisomer
% charge.....The charge on the input pseudoisomer

layers = regexp(inchi,'/','split');
f1 = layers{2}; % Fully protonated formula

p = {};
if ~isempty(strmatch('p',layers))
    p = layers{strmatch('p',layers)}; % nH to add or subtract
end

q = {};
if ~isempty(strmatch('q',layers))
    q = layers{strmatch('q',layers)}; % charge
end

if ~isempty(q)
   charge = str2double(q(2:end));
else
    charge = 0;
end

f1_nH = numAtomsOfElementInFormula(f1,'H'); % nH in fully protonated formula
if ~isempty(p)
    nH = f1_nH + str2double(p(2:end)); % nH in pseudoisomer formula
else
    nH = f1_nH;
end

formula = regexprep(f1,'H[a-z]*[0-9]*',''); % Remove all H from fully protonated formula
if nH == 1
    formula = [formula 'H'];
elseif nH > 1
    formula = [formula 'H' num2str(nH)]; % Add appropriate nH back in to create pseudoisomer formula
elseif nH < 0
    error('Negative number of H in formula.') % Should never get here
end

% In case there is Hg in formula
f1_nHg = numAtomsOfElementInFormula(f1,'Hg');
if f1_nHg == 1
    formula = [formula 'Hg'];
elseif f1_nHg > 1
    formula = [formula 'Hg' num2str(f1_nHg)];
end

end

author: opencobra

repo: cobratoolbox


text:
% function [datapeaks datavalleys] = findpeak(plotswitch,datain)
% This program finds peaks
% Arguments:
%   ploswitch: 0 or 1 : disable or enable plots (for debugging)
%   datain: 1D variable
% Dependencies:
%
% Copyright (c) 2011, Arun Ramakrishnan
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in the
%       documentation and/or other materials provided with the distribution.
%     * Neither the name of the <organization> nor the
%       names of its contributors may be used to endorse or promote products
%       derived from this software without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
% DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
% DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
% (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
% LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
% ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

function [datapeaks datavalleys] = findinflections(plotswitch,datain)
% datain = gconv;
datad1=[diff(datain)];
mp = find(datad1>0);
mn = find(datad1<=0);

datapeaks = [];
for pindex = 1:length(mp)
    [maxvalue, maxindex] = max(datain(mp(pindex):mn(min(find(mn>mp(pindex))))));
    if ~isempty(maxindex)
        datapeaks = cat(1,datapeaks, floor(mean(maxindex))+mp(pindex)-1);
    end
end
datapeaks= unique(datapeaks);

datavalleys = [];
for pindex = 1:length(mn)
    [minvalue, minindex] = min(datain(mn(pindex):mp(min(find(mp>mn(pindex))))));
    if ~isempty(minindex)
        datavalleys = cat(1,datavalleys, floor(mean(minindex))+mn(pindex)-1);
    end
end
datavalleys= unique(datavalleys);

if (plotswitch)
    plot(1:length(datain),datain,'k-',mp,datain(mp),'b.',mn,datain(mn),'r.',datapeaks,datain(datapeaks),'g.',datavalleys,datain(datavalleys),'c.');
end


text:
function IND = determine_search_location(A,method,params)

% determine the search location for updating each spatial component
% INPUTS:
% A:        current estimate of spatial component ( d x nr sparse matrix)
% method:   method to be used of determining search locations
%     available methods:
%       'ellipse': draw an ellipse centered at the center of mass with
%           axes the first two principal components and expand it to
%           determine the search location  (default)
%       'dilate' : grow the spatial footprint by dilating the current
%           footprint
%      if a different argument is passed then search location covers the whole field of view
% params:   hyper-parameter struct for the various methods (see default settings below for description)
%
% OUTPUT:
% IND:      binary d x nr matrix. IND(i,j) = 1 if pixel i is included in the search location of component j
%
% Written by Eftychios A. Pnevmatikakis, Simons Foundation
%   with input from Weijian Yang, Columbia University

if nargin < 3
    params = [];
end

if nargin < 2 || isempty(method)
    method = 'ellipse';         % default method
end

[d,nr] = size(A);

if ~isfield(params,'d1'); d1 = sqrt(d); params.d1 = d1; else d1 = params.d1; end          % # of rows
if ~isfield(params,'d2'); d2 = sqrt(d); params.d2 = d2; else d2 = params.d2; end          % # of columns
if ~isfield(params,'d3'); d3 = 1; params.d3 = d3; else d3 = params.d3; end
if ~isfield(params,'min_size') || isempty(params.min_size); min_size = 3; else min_size = params.min_size; end    % minimum size of ellipse axis
if ~isfield(params,'max_size') || isempty(params.max_size); max_size = 8; else max_size = params.max_size; end    % maximum size of ellipse axis
if ~isfield(params,'dist')  || isempty(params.dist); dist = 3; else dist = params.dist; end                              % expansion factor of ellipse
if ~isfield(params,'bSiz')  || isempty(params.dist); bSiz = 4; else bSiz = params.bSiz; end                              % expansion factor of ellipse
if ~isfield(params,'se')  || isempty(params.se);  % morphological element (for 'dilate')
    if d3 == 1; expandCore = strel('disk',bSiz,0); else expandCore = strel(ones(bSiz,bSiz,2)); end
else
    expandCore = params.se;
end

if strcmpi(method,'ellipse'); method = 'ellipse';
elseif strcmpi(method,'dilate'); method = 'dilate';
else fprintf('Method not recongnized. Search location equals the entire field of view. \n');
end

IND = false(d,nr);
% find all zero components
ind_empty = (sum(A, 1)==0);
if any(ind_empty)
    A(1, ind_empty) = 1;
end
switch method
    case 'ellipse'
        Coor.x = kron(ones(d2*d3,1),(1:d1)');
        Coor.y = kron(ones(d3,1),kron((1:d2)',ones(d1,1)));
        Coor.z = kron((1:d3)',ones(d1*d2,1));
        if ~(dist==Inf)             % determine search area for each neuron
           %cm = zeros(nr,2);        % vector for center of mass
           cm = com(A,d1,d2,d3);
           if d3 == 1
               cm = [cm,ones(nr,1)];
           end
           Vr = cell(nr,1);
           IND = zeros(d,nr);       % indicator for distance
           %cm(:,1) = Coor.x'*A(:,1:nr)./sum(A(:,1:nr));
           %cm(:,2) = Coor.y'*A(:,1:nr)./sum(A(:,1:nr));          % center of mass for each components
           parfor i = 1:nr            % calculation of variance for each component and construction of ellipses
               if d3 == 1
                   Vr{i} = ([Coor.x - cm(i,1), Coor.y - cm(i,2)]'*spdiags(A(:,i),0,d,d)*[Coor.x - cm(i,1), Coor.y - cm(i,2)])/sum(A(:,i));
                   [V,D] = eig(Vr{i});
                   cor = [Coor.x - cm(i,1),Coor.y - cm(i,2)];
               else
                   Vr{i} = ([Coor.x - cm(i,1), Coor.y - cm(i,2), Coor.z - cm(i,3)]'*spdiags(A(:,i),0,d,d)*[Coor.x - cm(i,1), Coor.y - cm(i,2), Coor.z - cm(i,3)])/sum(A(:,i));
                   [V,D] = eig(Vr{i});
                   cor = [Coor.x - cm(i,1),Coor.y - cm(i,2),Coor.z - cm(i,3)];
               end
               d11 = min(max_size^2,max(min_size^2,D(1,1)));
               d22 = min(max_size^2,max(min_size^2,D(2,2)));
               if d3 == 1
                   IND(:,i) = sqrt((cor*V(:,1)).^2/d11 + (cor*V(:,2)).^2/d22)<=dist;
               else
                   d33 = min((max_size/2)^2,max((min_size/2)^2,D(3,3)));
                   IND(:,i) = sqrt((cor*V(:,1)).^2/d11 + (cor*V(:,2)).^2/d22 + (cor*V(:,3)).^2/d33)<=dist;       % search indeces for each component
               end
           end
        else
            IND = true(d,nr);
        end
    case 'dilate'
        A = threshold_components(A,params);
        parfor i = 1:nr
            A_temp = imdilate(reshape(full(A(:,i)),d1,d2,d3),expandCore);
            IND(:,i) = A_temp(:)>0;
        end
    otherwise
        IND = true(d,nr);
end
if any(ind_empty)
    IND(:, ind_empty) = false;
end

author: zhoupc


text:
function aux_FMD(params, hParentFigure)
% function aux_FMD(params, hParentFigure);
%-------------------------------------------
% Plot determination of Mc for a grid point using  calc_McCdf_plot
%
% Incoming variables:
% params        : all variables
% hParentFigure : Handle of the parent figure
%
% J.Woessner, woessner@seismo.ifg.ethz.ch
% last update: 23.01.03

% Get the axes handle of the plotwindow
axes(sv_result('GetAxesHandle', hParentFigure, [], guidata(hParentFigure)));
hold on;
% Select a point in the plot window with the mouse
[fX, fY] = ginput(1);
disp(['X: ' num2str(fX) ' Y: ' num2str(fY)]);
% Plot a small circle at the chosen place
plot(fX,fY,'ok');

% Get closest gridnode for the chosen point on the map
[fXGridNode fYGridNode,  nNodeGridPoint] = calc_ClosestGridNode(params.mPolygon, fX, fY);
plot(fXGridNode, fYGridNode, '*r');
hold off;

% Get the data for the grid node
mNodeCatalog_ = params.mCatalog(params.caNodeIndices{nNodeGridPoint}, :);

% Start calculation
%[mResult, fMls, fMc, fMu, fSigma, mDatPredBest, vPredBest] = plot_McEMR(mNodeCatalog_, 0.1)
[mResult, fMls, fMc, fMu, fSigma, mDatPredBest, vPredBest] =plot_McCdfnormal(mNodeCatalog_, 0.1)
%plot_McCdf2(mNodeCatalog_, 0.1);


text:
filename='Cantileverbeam_Tetrahedra_Linear_Structured_SYM';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1];
constraint = {'volume'};
optimizer = 'MMA';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'INCREASING LAST STEP';
showBC = true;

nsteps = 10;
Vfrac_final = 0.1;
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:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "Hill+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the target point positions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function targets = update_Target_Point_Positions(dt,current_time,targets)

IDs = targets(:,1);                 % Stores Lag-Pt IDs in col vector
%xPts= targets(:,2);                % Original x-Values of x-Target Pts.
yPts= targets(:,3);                 % Original y-Values of y-Target Pts.
%kStiffs = targets(:,4);             % Stores Target Stiffnesses

N_target = length(targets(:,1));    %Gives total number of target pts!

% NEEDED FOR UPDATE FILE FROM CORAL GEOMETRY! %
N_half = 90;                        % # of pts. on ONE tentacle

    %frames per second in data
     F = 30.0;
     %pixels per dm in data
     S = 1733.333333333333333333;
     %time in frames (mod one pulse)
     t =  F*mod(current_time,2.30);
     tp = F*2.25;
     tend = F*2.30;
     %initial offset from center of the domain in dm
     xoffset = 0.019;
     yoffset = 0.25;
     xoffset2grid1 = 0.2;
     xoffset2grid2 = 0.4;

     qc11 = -8.018203211249667e-13;
     qc12 = 1.706310971452060e-10;
     qc13 = -1.371526883228119e-8;
     qc14 = 5.116506115621489e-7;
     qc15 = -8.613944795521634e-6;
     qc16 = 4.951095471423292e-5;
     qc17 = -1.237254488643602e-4;
     qc21 = -1.093988700839764e-10;
     qc22 = 2.370737151797367e-8;
     qc23 = -1.946684756263760e-6;
     qc24 = 7.447999061932255e-5;
     qc25 = -1.285757196003425e-3;
     qc26 = 7.184871148722453e-3;
     qc27 = -2.357284935085283e-2;
     qc31 = -3.111491615328405e-9;
     qc32 = 6.922752279229525e-7;
     qc33 = -5.951659336196739e-5;
     qc34 = 2.461541908755476e-3;
     qc35 = -4.830999119980220e-2;
     qc36 = 3.293063238297988e-1;
     qc37 = -1.893567487439901;
     qc41 = 2.028150369106308e-8;
     qc42 = -4.188949936239104e-6;
     qc43 = 3.237380275082298e-4;
     qc44 = -1.158763275070524e-2;
     qc45 = 1.908244836752949e-1;
     qc46 = -1.142372536532416;
     qc47 = 3.584246435337805;

    if (t<=tp)

        qpoly1=qc17+qc16*t+qc15*t*t+qc14*t*t*t+qc13*t*t*t*t+qc12*t*t*t*t*t+qc11*t*t*t*t*t*t;
        qpoly2=qc27+qc26*t+qc25*t*t+qc24*t*t*t+qc23*t*t*t*t+qc22*t*t*t*t*t+qc21*t*t*t*t*t*t;
        qpoly3=qc37+qc36*t+qc35*t*t+qc34*t*t*t+qc33*t*t*t*t+qc32*t*t*t*t*t+qc31*t*t*t*t*t*t;
        qpoly4=qc47+qc46*t+qc45*t*t+qc44*t*t*t+qc43*t*t*t*t+qc42*t*t*t*t*t+qc41*t*t*t*t*t*t;

    else
        qpolytf1=qc17+qc16*tp+qc15*tp*tp+qc14*tp*tp*tp+qc13*tp*tp*tp*tp+qc12*tp*tp*tp*tp*tp+qc11*tp*tp*tp*tp*tp*tp;
        qpolytf2=qc27+qc26*tp+qc25*tp*tp+qc24*tp*tp*tp+qc23*tp*tp*tp*tp+qc22*tp*tp*tp*tp*tp+qc21*tp*tp*tp*tp*tp*tp;
        qpolytf3=qc37+qc36*tp+qc35*tp*tp+qc34*tp*tp*tp+qc33*tp*tp*tp*tp+qc32*tp*tp*tp*tp*tp+qc31*tp*tp*tp*tp*tp*tp;
        qpolytf4=qc47+qc46*tp+qc45*tp*tp+qc44*tp*tp*tp+qc43*tp*tp*tp*tp+qc42*tp*tp*tp*tp*tp+qc41*tp*tp*tp*tp*tp*tp;

        qpoly1 = qpolytf1*((tend-tp)-(t-tp))/(tend-tp)+qc17*(t-tp)/(tend-tp);
        qpoly2 = qpolytf2*((tend-tp)-(t-tp))/(tend-tp)+qc27*(t-tp)/(tend-tp);
        qpoly3 = qpolytf3*((tend-tp)-(t-tp))/(tend-tp)+qc37*(t-tp)/(tend-tp);
        qpoly4 = qpolytf4*((tend-tp)-(t-tp))/(tend-tp)+qc47*(t-tp)/(tend-tp);
    end

for i=1:N_target                    % Loops over all target points!

    shifted_Xt1 = yPts(i)-yoffset;

    X_target = xoffset+qpoly4/S-qpoly3*shifted_Xt1+qpoly2*S*shifted_Xt1*shifted_Xt1-qpoly1*S*S*shifted_Xt1*shifted_Xt1*shifted_Xt1;

    if i<= N_target/2
        if i<=N_half
            X_target = -X_target+ xoffset2grid1;
        else
            X_target = X_target + xoffset2grid1;
        end
    else
        if i<=N_target/2+N_half
            X_target = -X_target+ xoffset2grid2;
        else
            X_target = X_target + xoffset2grid2;
        end
    end
    targets(IDs(i),2) = X_target;      % Store new xVals
    %targets(IDs(i),3) = yPts(IDs(i)); % Store new yVals

end

author: nickabattista


text:
function [marginals, marginalsT] = project_joint_onto_marginals(joint, hnodes, ns)

ss = length(ns);
jointT = dpot(hnodes, ns(hnodes), joint);
marginalsT = cell(1, ss);
marginals = cell(1,ss);
for i=hnodes(:)'
  marginalsT{i} = marginalize_pot(jointT, i);
  m = pot_to_marginal(marginalsT{i});
  marginals{i} = m.T(:);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt