text:
%% RPCA demo using TFOCS
% This is a demo of Robust PCA (RPCA) using <http://tfocs.stanford.edu TFOCS>. Since 2009,
% there has been much interest in this specific RPCA formulation
% (RPCA can refer to many different formulations; we will state
% our formulation precisely below), but most solvers
% assume equality constraints. One great feature of TFOCS
% is how easy it is to prototype new formulations, so we show
% here a formulation using inequality constraints.
%
% The basic (equality constraint) RPCA formulation is:
% ( see <http://perception.csl.uiuc.edu/matrix-rank/references.html here>
% for references )
%
% $$ \min_{S,L} \|L\|_* + \lambda \|S\|_{\ell_1} \quad\textrm{subject to}\;L+S = X $$
%
% The idea is that some object $X$ can be broken into two components $X=L+S$,
% where $L$ is low-rank and $S$ is sparse.  To tease out this separation,
% RPCA penalizes $\|L\|_*$ (the nuclear norm of $L$), which is a proxy for
% the rank of $L$.  We also penalizes $\|S\|_{\ell_1}$ which is just the sum
% of the absolute value of all the entries of the matrix $S$; this is a proxy
% for the number of non-zero entries of $S$.
%
% The parameter $\lambda$ is important, as it controls how much weight
% to put on $S$ relative to $L$.  This parameter is not overly sensitive,
% and we picked the value used in this demo by a few quick rounds of trial-and-error.
%
% _Code by Stephen Becker, May 2011_

%% Inequality constraints
% In this demo, $X$ is a video. In particular, each column of $X$ is the set
% of $m n$ pixels from a $m \times n$ pixel frame of the video, reshaped
% so that it is now a vector instead of a matrix.
%
% The particular video we use is not compressed using a video codec, so
% each frame of the video is stored. This leads to huge files, so to save
% a bit of space, the value of each pixel is stored as an 8 bit number,
% which means we can think of it as a value from $0,1,2,\ldots,255$.  Now,
% we may imagine that there is some "true" video where each pixel
% is a _real_ number between $[0,255]$, and that we see a quantized version of it.
% It would be nice to apply RPCA to this true version instead of applying
% it to the quantized version.  In particular, it is unlikely that the quantized
% version can be split nicely into a low-rank and sparse component.
% To account for this quantization, we will ask that $S+L$ is not exactly equal
% to $X$, but rather that $S+L$ agrees with $X$ up to the precision
% of the quantization.  The quantization can induce at most an error of .5 in
% the pixel value (e.g. true value is 5.6, which is rounded to 6.0, so an error of .4;
% or a true value of 6.4, which is also rounded to 6.0, so also an error of 0.4).
% A nice way to capture this is via the $\ell_\infty$ norm.
%
% So, the inequality constrained version of RPCA uses the same objective
% function, but instead of the constraints $L+S=X$, the constraints are
%
% $$ \|L+S-X\|_{\ell_\infty} \le 0.5. $$

%% Numerical demo
% We demonstrate this on a video clip taken from a surveillance
% camera in a subway station.  Not only is there a background
% and a foreground (i.e. people), but there is an escalator
% with periodic motion.  Conventional background subtraction
% algorithms would have difficulty with the escalator, but
% this RPCA formulation easily captures it as part of the low-rank
%   structure.
% The clip is taken from the data at this website (after
% a bit of processing to convert it to grayscale):
% http://perception.i2r.a-star.edu.sg/bk_model/bk_index.html

% Load the data:
load escalator_data % contains X (data), m and n (height and width)
X = double(X);

% addpath ~/Dropbox/TFOCS/   % add TFOCS to your path, so change this to suit your computer

%%
nFrames     = size(X,2);

lambda  = 1e-2;

opts = [];
opts.stopCrit       = 4;
opts.printEvery     = 1;
opts.tol            = 1e-4;

opts.maxIts         = 25;

opts.errFcn{1}      = @(f,d,p) norm(p{1}+p{2}-X,'fro')/norm(X,'fro');

largescale      = false;

for inequality_constraints = 0:1

    if inequality_constraints
        % if we already have equality constraint solution,
        %   it would make sense to "warm-start":
%         x0      = { LL_0, SS_0 };
        % but it's more fair to start all over:
        x0      = { X, zeros(size(X))   };
        z0      = [];
    else
        x0      = { X, zeros(size(X))   };
        z0      = [];
    end

    obj    = { prox_nuclear(1,largescale), prox_l1(lambda) };
    affine = { 1, 1, -X };

    mu = 1e-4;
    if inequality_constraints
        epsilon  = 0.5;
        dualProx = prox_l1(epsilon);
    else
        dualProx = proj_Rn;
    end

    tic
    % call the TFOCS solver:
    [x,out,optsOut] = tfocs_SCD( obj, affine, dualProx, mu, x0, z0, opts);
    toc

    % save the variables
    LL = x{1};
    SS = x{2};
    if ~inequality_constraints
        z0      = out.dual;
        LL_0    = LL;
        SS_0    = SS;
    end

end % end loop over "inequality_constriants" variable

%% show all together in movie format
% If you run this in your own computer, you can see the movie. On the webpage,
%   we have a youtube version of the video.
% The top row is using equality constraints, and the bottom row
% is using inequality constraints.
%  The first column of both rows is the same (i.e. it is the original image).
mat  = @(x) reshape( x, m, n );
figure();
colormap( 'Gray' );
k = 1;
for k = 1:nFrames

    imagesc( [mat(X(:,k)), mat(LL_0(:,k)),  mat(SS_0(:,k)); ...
              mat(X(:,k)), mat(LL(:,k)),    mat(SS(:,k))  ] );

    axis off
    axis image

    drawnow;
    pause(.05);

    if k == round(nFrames/2)
        snapnow; % Take a single still snapshot for publishing the m file to html format
    end
end

%% Is there a difference between the two versions?
% Compare the equality constrained version and the inequality
%   constrained version.  To do this, we can check whether
%   the "L" components are really low rank
%   and whether the "S" components are really sparse.
% Even if the two versions appear visually similar, the variables
%   may behave quite differently with respect to low-rankness
%   and sparsity.

fprintf('"S" from equality constrained version has  \t%.1f%% nonzero entries\n',...
    100*nnz(SS_0)/numel(SS_0) );
fprintf('"S" from inequality constrained version has  \t%.1f%% nonzero entries\n',...
    100*nnz(SS)/numel(SS) );

s  = svd(LL);    % inequality constraints
s0 = svd(LL_0);  % equality constraints

fprintf('"L" from equality constrained version has numerical rank\t %d (of %d possible)\n',...
    sum( s0>1e-6), min( m*n, nFrames)  );
fprintf('"L" from inequality constrained version has numerical rank\t %d (of %d possible)\n',...
    sum( s > 1e-6), min( m*n, nFrames)  );

figure();
semilogy( s0 ,'o-')
hold all;
semilogy( s ,'o-')
legend('Equality constraints','Inequality constraints');
xlabel('sorted singular value location');
ylabel('value of singular value');
title('Comparison of RPCA using equality and inequality constraints');

%%


text:
%% (Internal) Display a message framed to a string
%
%     This function format a string for fprintf like functions in order to
%     present a message string2display framed with # characters
%
%   str_out = disp_string_framed(style, string2display )
%
% Arguments:
%
%      + style: if numeric, could be the stream identifier for fprintf
%      function. If char, can be the "style" argument for cprintf function.
%
%      + string2display: The string to frame.
%
% Output:
%
%      + str_out: resulting string
%
% Example:
%
%
% See also disp_option_enumeration, disp_string_title
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function str_out = disp_string_framed(style, string2display )

    str_out = [];

    if( nargin == 1 || isempty(string2display) )
        string2display = style;
        style = 'text';
        fid = 1;
    end

    if( isempty(style) )
        style = 'text';
        fid = 1;
    elseif( isnumeric(style) )
        fid = style;
        style = 'text';
    else
        if( ischar(style) )
            fid = 1;
        else
            style = 'text';
            fid = 1;
        end
    end

    str_aux = sprintf('# %s #', string2display);
    str_aux1 = repmat('#', 1, length(str_aux));
    if( fid == 0)
        string2display = strrep(string2display, '\', '\\');
        str_out = sprintf('\n%s\n# %s #\n%s\n\n', str_aux1, string2display, str_aux1);
    elseif( fid == 1)
        fprintf(fid, '\n%s\n# ', str_aux1);
        string2display = strrep(string2display, '\', '\\');
        cprintf(style, [string2display ' ']);
        fprintf(fid, '#\n%s\n\n', str_aux1);
    else
        fprintf(fid, '\n%s\n# %s #\n%s\n\n', str_aux1, string2display, str_aux1);
    end

end

author: marianux

repo: ecg-kit

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


text:
function value = halham_leap_check ( dim_num, leap )

%*****************************************************************************80
%
%% HALHAM_LEAP_CHECK checks LEAP for a Halton or Hammersley sequence.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer LEAP(1:DIM_NUM), the leap vector.
%    Each entry should be 1 or greater.  Only the integer
%    part of a base is used.
%
%    Output, logical VALUE, is true if LEAP is legal.
%
  if ( any ( leap(1:dim_num) < 1 ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HALHAM_LEAP_CHECK - Fatal error!\n' );
    fprintf ( 1, '  At least one of the input leap entries is <= 1!\n' );
    value = 0;
  else
    value = 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [limits state] = hlp_adaptMinMaxLimits(varargin)
% adapt min and max limits using an exponential window moving average. The
% moving average is calculated in place based on new data values and
% previous min/max limits.

g = arg_define(varargin, ...
        arg_norep({'values'},[],[],'data values. can be scalar, vector or matrix'), ...
        arg_nogui({'state'},[],[],'current state of the adaptation'), ...
        arg_sub({'adaptOpts'},{},@hlp_scaleLimits,'Adaptation options'), ...
        arg({'reset'},false,[],'Reset adaptation state') ...
        );

state = g.state;

if g.reset || isempty(state)
    % initialize state
    state.lastMin      = min(g.values(:));
    state.lastMax      = max(g.values(:));
    state.numRunsSoFar = 0;
end

if isfield(g.adaptOpts,'arg_direct')
    g.adaptOpts = rmfield(g.adaptOpts,'arg_direct');
end

% adapt limits using exponential window moving average
[state.lastMin state.lastMax] = hlp_scaleLimits(g.values,state.lastMin, ...
                                state.lastMax,state.numRunsSoFar,g.adaptOpts);
state.numRunsSoFar = state.numRunsSoFar + 1;

limits = [state.lastMin, state.lastMax];

repo: REST


text:
% This program illustrates various examples presented in
% 'Advanced Mathematics and Mechanics Applications Using MATLAB'
% 3rd Ed., CRC Press, 2003, by Howard Wilson, Louis Turcotte and
% David Halpern. The program is executed by typing rundynam.
% Then select the individually numbered examples. The separate
% program modules are listed below.
%
% rundynam   main program which calls the other modules
% beamwave   vibrating beam with initial deflection
% beszeros   zeros of integer order Bessel functions
% circwave   circular membrane with an oscillating force
% elipmodes  elliptic membrane modes by Moler's method
% forcmove   moving force on a semi-infinite string
% rectwave   rectangular membrane with an oscillating force
% runchain   falling elastic chain with one end fixed
% runpenl    pendulum with large amplitude oscillations
% runtop     spinning and nutating top motion
% smdplot    spring-mass-damper with harmonic applied force
% stringmo   vibrating string shaken at one end
% strngwave  vibrating string with initial deflection
% trusvibs   natural frequencies of a pin-connected truss

author: Sable

repo: mcbench-benchmarks


text:
function pd = prtEvalPdAtPf(classifier,dataSet,pfDesired,nFolds)
% prtEvalPdAtPf   Returns the probability of false alarm at a desired
%   probability of detection on the receiver operating curve.
%
%   PD = prtEvalPdAtPf(prtClassifier, prtDataSet, pfDesired) returns the
%   probabilty of detection PD.  prtDataSet must be a labeled, binary
%   prtDataSetStandard object. prtClassifier must be a prtClass object.
%   pfDesired is the desired probability of false alarm and must be between
%   0 and 1.
%
%   pfDesired = prtEvalPdAtPf(prtClassifier, prtDataSet, pfDesired, nFolds)
%   returns the probabilty of false alarm pfDesired on the receiver
%   operating curve with K-fold cross-validation. prtDataSet must be a
%   labeled, binary prtDataSetStandard object. prtClassifier must be a
%   prtClass object. pfDesired is the desired probability of detection and
%   must be between 0 and 1. nFolds is the number of folds in the K-fold
%   cross-validation.
%
%   pfDesired = prtEvalPdAtPf(prtClassifier, prtDataSet, pfDesired,
%   xValInds) same as above, but use crossValidation with specified indices
%   instead of random folds.
%
%   Example:
%   dataSet = prtDataGenSpiral;
%   classifier = prtClassDlrt;
%   pf = prtEvalPdAtPf(classifier, dataSet,.01)
%
%   See Also: prtEvalAuc, prtEvalPfAtPd, prtEvalPercentCorrect,
%   prtEvalMinCost

if nargin < 4 || isempty(nFolds)
    nFolds = 1;
end
results = prtUtilEvalParseAndRun(classifier,dataSet,nFolds);

[pf,pd] = prtScoreRoc(results.getObservations,dataSet.getTargets);

[pf,sortInd] = sort(pf(:),'ascend');
pd = pd(sortInd);

ind = find(pf >= pfDesired);
if ~isempty(ind)
    pd = pd(ind(1));
else
    pd = nan;
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% ir_ct_roi_split1.m
% "frequency split" approach to ROI recon of Lin Fu et al. from Fully 3D 2015.
% Low frequency part of ROI sinogram comes from reprojecting FBP ROI image,
% whereas high frequency part comes from original sinogram.
% ROI reconstructed using PWLS-OS-LALM
% 2015-06-17 Jeff Fessler, University of Michigan

if ~isvar('Af'), printm 'setup geometry'
	f.down = 4;
	igf = image_geom('nx', 512, 'fov', 50, 'down', f.down);
	igf.mask = igf.circ > 0;
	sg = sino_geom('ge1', 'units', 'cm', 'strip_width', 'd', ...
		'down', f.down);

	igr = igf;
	igr.mask = igr.circ(11, 11, -10, 4) > 0;
	im(igf.mask + igr.mask)

	% system objects
	if has_mex_jf
		Af = Gtomo2_dscmex(sg, igf); % full
		Ar = Gtomo2_dscmex(sg, igr); % roi
	else
		Af = Gtomo_nufft_new(sg, igf);
		Ar = Gtomo_nufft_new(sg, igr);
	end
end

if ~isvar('xtrue'), printm 'xtrue, sinogram'
	% read image
    xtrue256 = ir_get_data('ncat_256_slice_140_ct_x100.fld');
	xtrue256 = single(xtrue256) / 200 * 0.4; % convert to 1/cm units

	if 1 % more realistic sinogram from finer image, avoid "inverse crime"
		ig_big = image_geom('nx', 512, 'fov', igf.fov, 'down', 2);
		if has_mex_jf
			Abig = Gtomo2_dscmex(sg, ig_big);
		else
			Abig = Gtomo_nufft_new(sg, ig_big);
		end
		sino_true = Abig * xtrue256;
	end
	xtrue = downsample2(xtrue256, f.down/2);

	im plc 2 2
	clim = [0 0.4];
	im(1, xtrue, 'x true', clim), cbar
	xlabelf('units: 1 / %s', sg.units)
	im(2, sino_true, 'sino true'), cbar
	im(3, xtrue + igf.mask + igr.mask)

	clear ddir ig_big Abig
drawnow
end

mask2 = conv2(single(igr.mask), ones(9)/9^2, 'same') > 0.999; % for roi rmse
im(3, igf.mask + igr.mask + mask2, 'ROIs')
xl = @(x) xlabelf('RMSE = %.3f / %s', rms(col(x - xtrue)), sg.units);
xr = @(x) xlabelf('RMSE = %.3f / %s', rms(x(mask2) - xtrue(mask2)), sg.units);

if ~isvar('sino'), printm 'noisy sinogram'
	rng(0)
	% transmission data:
	I0 = 1e5;
	yi = poisson(I0 * exp(-sino_true), 0, 'factor', 0.1);
	if any(yi(:) == 0)
		warn('%d of %d values are 0 in sinogram!', ...
			sum(yi(:)==0), length(yi(:)));
	end
	sino = log(I0 ./ max(yi,1)); % noisy fan-beam sinogram
	im(4, sino, 'sino noisy'), cbar
drawnow
%	ir_savefig ir_ct_roi_split1a
end

if ~isvar('fbp'), printm 'fbp 2d fan-beam reconstruction'
	tmp = fbp2(sg, igf);
	fbp = fbp2(sino, tmp, 'window', 'hanning,0.75');
	im(2, fbp, 'FBP Hanning', clim), cbar
	xl(fbp)
prompt
end

if ~isvar('kappa'), printm 'kappa: try to make resolution approximately uniform'
	wi = yi; % will give 0 weight to any ray where yi=0!
	kappa = sqrt( div0(Af' * wi, Af' * ones(size(wi))) );
	im(3, kappa, 'kappa'), cbar
prompt
end

% use local psf to help select beta
if ~isvar('R'), printm 'R'
	f.l2b = 10; % maybe a bit too big, but ok for now
	f.delta = 0.001;
%	f.pot_arg = {'lange3', f.delta}; % todo: why not as sharp as hyper3?
	f.pot_arg = {'hyper3', f.delta};
	R = Reg1(kappa, 'beta', 2^f.l2b, 'pot_arg', f.pot_arg);
	Rr = Reg1(kappa .* igr.mask, 'beta', 2^f.l2b, 'pot_arg', f.pot_arg);
%	qpwls_psf(A, R, 1, igf.mask, Gdiag(wi), 'loop', 1); % use this to choose beta
end

f.niter = 6;
f.nblock = 41; % 41 subsets

% OS-LALM
if ~isvar('xlalmf'), printm 'iterative reconstruction - lalm full'
	Ab = Gblock(Af, f.nblock);
	xlalmf = ir_pwls_os_lalm(fbp(igf.mask), Ab, sino, R, 'wi', wi, ...
		'isave', 'last', 'niter', f.niter);
	xlalmf = igf.embed(xlalmf);
	im(4, xlalmf(:,:,end), 'LALM full', clim), cbar
	xl(xlalmf(:,:,end))
%	ir_savefig ir_ct_roi_split1b
end

if ~isvar('fbpr'), printm 'fbp 2d fan-beam reconstruction - roi'
	tmp = fbp2(sg, igr);
	fbpr = fbp2(sino, tmp);
	im(1, xtrue .* igr.mask, 'True roi', clim), cbar
	im(2, fbpr, 'FBP roi', clim), cbar
	xr(fbpr)
prompt
end

if ~isvar('Hhi'), printm 'filters'
	nf = 2 * sg.nb;
	u = [-nf/2:nf/2-1]'/nf;
	cut = 0.1;
	Hhi = min((u/cut).^2, 1);
	Hlo = 1 - Hhi;
	plot(u, Hlo, '-o')
	clear u
end

if ~isvar('sino_roi'), printm 'sino_roi'
	sino_f = fft(sino, nf, 1);
	sino_hi = ifft(sino_f .* repmat(ifftshift(Hhi), [1 sg.na]), [], 1);
	sino_hi = sino_hi(1:sg.nb, :);
	im(1, fftshift(sino_f, 1))
	im(2, sino_hi, 'Hi-pass sino')

	sino_Ar = Ar * fbpr;
	sino_Arf = fft(sino_Ar, nf, 1);
	im(1, sino_Ar, 'Reproj FBP')
        sino_lo = ifft(sino_Arf .* repmat(ifftshift(Hlo), [1 sg.na]), [], 1);
	sino_lo = sino_lo(1:sg.nb, :);
	im(3, sino_lo, 'Low-pass sino')

	sino_roi = sino_lo + sino_hi;
	im(4, sino_roi, 'Synth ROI sino')
	clear sino_f sino_Arf
%	ir_savefig ir_ct_roi_split1c
prompt
end

% OS-LALM
if ~isvar('xlalmr'), printm 'iterative reconstruction - lalm roi'
	Ab = Gblock(Ar, f.nblock);
	xlalmr = ir_pwls_os_lalm(fbpr(igr.mask), Ab, sino_roi, Rr, 'wi', wi, ...
		'isave', 'last', 'niter', 2*f.niter);
	xlalmr = igr.embed(xlalmr);
end

if 1 % pics
	im(4, xlalmr(:,:,end), 'PWLS-OS-LALM roi', clim), cbar
	xr(xlalmr(:,:,end))

	im(1, xtrue .* igr.mask, 'True roi', clim), cbar
	im(2, fbpr, 'FBP Ramp roi', clim), cbar
	xr(fbpr)
	im(3, fbp .* igr.mask, 'FBP Hanning roi', clim), cbar
	xr(fbp)
%	ir_savefig ir_ct_roi_split1d
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%kkimagpart 'Output = Imaginary Component of (Complex) Input'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kimagpart.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'input data object'
% OutputFile: oi 'Output', required: 'resulting output data object'
%
% Example: oi = kkimagpart(i, {'i','';'oi',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kimagpart - Output = Imaginary Component of (Complex) Input
%
%  DESCRIPTION
% The "Imaginary Part" operator returns the imaginary component of each
% complex data point in the input data object, \fBInput".  If the input
% data is not complex, the imaginary component is assumed to be zero.
%
% The output object data type will be the same as the
% input (regarding single or double floats).
%
%  "Map Data" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_1input
%
%  "Validity Mask" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "Explicit Location and Time Data" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "Failure Modes"
% .cI $DATAMANIP/repos/shared/man/sections/fail_1input
%
% Executing "Imaginary Part" runs the program \fIkcmplx2real\fP with the -oi
% output option.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::kcmplx2real, DATAMANIP::kcmplx, DATAMANIP::kreal2cmplx
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% TFSA PLOTTING FACILITIES
%
% Plotting routines used to visualise data. The following types of
% graphs can be employed:
%
%  (i)   Plot
%  (ii)  Image
%  (iii) Pseudo-Colour
%  (iv)  Waterfall (1 and 2)
%  (v)   Mesh
%  (vi)  Surf
%  (vii) Contour
%  (viii)TFSAPL
%
% All plots except (i) can be used to plot time-frequency
% distributions. The time domain and frequency spectra one
% dimensional plots can be added to the figure along the two axes.
% The plots can be configured by specifying various options in the
% figure. The 'Expand' button can be used to reveal more options. On
% the actual plot figure the appearance of the plots can be altered
% from a popup menu which is activated by right clicking.
%
% EXPLANATION OF VARIOUS PLOTTING ROUTINES
%
%  (i) Plot
%
%      Used for one dimensional plots of signals, IF estimation,
%      wavelet coefficients, frequency spectra etc.
%
%  (ii) Image
%
%      Used to plot two dimensional matrices such as time-frequency
%      distributions and scalograms. Displays the matrix as an image.
%
%  (iii) Pseudo-Colour
%
%      Pseudo-colour plots a colour grid for each element of the input
%      matrix. Also used to plot two dimensional matrices such as
%      time-frequency distributions and scalograms.
%
%  (iv) Waterfall 1 and 2
%
%      Again also used to plot two dimensional matrices such as
%      time-frequency distributions and scalograms. This uses a
%      waterfall style plot of the matrix, using different viewing
%      angles which distinguish Waterfall 1 and Waterfall 2.
%
%  (v) Mesh
%
%      Same as the waterfall plots except that column lines are
%      drawn to produce a 'mesh' surface.
%
%  (vi) Surf
%
%      Equivalent to the mesh plot except that a 'surface' is now
%      defined over the mesh.
%
%  (vii) Contour
%
%      Contour plot of the two dimensional matrix where the
%      different coloured lines represent varying degrees of magnitudes.
%
%  (viii) TFSAL
%
%      TFSA graphing function specifically optimised for viewing
%      time-frequency distributions.
%
%
%
%
%  See Also:  tfsapl

author: HeLiangHIT


text:
function Vo = spm_file_split(V, odir)
% Convert a 4D volume file into a series of 3D volume files
% FORMAT Vo = spm_file_split(V, odir)
% V         - filename or spm_vol struct
% odir      - output directory [default: same as input]
%
% Vo        - spm_vol struct array of output files
%__________________________________________________________________________
% Copyright (C) 2009-2018 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_file_split.m 7391 2018-08-13 09:55:17Z guillaume $

if ~nargin
    [V, sts] = spm_select(1,'nifti','Select a 4D volume file to split');
    if ~sts, return; end
end
if ischar(V)
    V = spm_vol(spm_file(V,'number',''));
end

[p,n,e] = spm_fileparts(V(1).fname);
if nargin < 2
    if isempty(p), p = pwd; end
    odir = p;
end

Voo = cell(numel(V),1);
spm_progress_bar('Init',numel(V),'Splitting 4D Volume','Volumes Complete');
for i=1:numel(V)
    Voo{i}         = fullfile(odir,sprintf('%s_%05d%s',n,i,e));
    ni             = nifti;
    ni.dat         = file_array(Voo{i},V(i).dim(1:3),V(i).dt,0,V(i).pinfo(1),V(i).pinfo(2));
    ni.mat         = V(i).mat;
    ni.mat0        = V(i).mat;
    ni.mat_intent  = V(i).private.mat_intent;
    ni.mat0_intent = V(i).private.mat0_intent;
    ni.descrip     = [V(i).descrip sprintf(' - %d',i)];
    create(ni);
    ni.dat(:,:,:)  = V(i).private.dat(:,:,:,i);
    spm_progress_bar('Set',i);
end
spm_progress_bar('Clear');

if nargout
    Vo = spm_vol(char(Voo));
end

author: spm


text:
%  INTERNAL FUNCTION: recast the exit flag from fsolve and lsqnonlin
%
%  ::
%
%    flag=exitflag(flag,x1,fx1,TolFun)
%
%  Args:
%
%     - **flag** [scalar]: exit flag from fsolve or lsqnonlin
%     - **x1** [scalar|vector|matrix]: optimum of the function to set to zero
%       (f(x)=0)
%     - **fx1** [scalar]: value or norm of f(x1)
%     - **TolFun** [{sqrt(eps)}|scalar]: scalar or vector of model objects
%
%  Returns:
%     :
%
%     - **flag** [scalar]: recomputed exitflag
%
%

author: jmaih


text:
function showgraphmatrix(G,node)
%% SHOWGRAPHMATRIX displays a planar graph
%
%    showgraphmatrix(G,node) displays a planar undirected graph represented
%    by a sparse matrix G.
%
%   See also showgraph, showmesh, findedge
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

[i,j] = find(triu(G,1));
showgraph(node,[i j]);

repo: ifem


text:
function vw = spatialBlurTSeries(vw, scanList, kernelSize_mm, newTypeName)
% vw = spatialBlurTSeries(vw, scanList, kernelSize_mm, newTypeName)
%
% Performs spatial filtering on the time series data.
% Useful only in the inplane view.
%
% If 'dialog' is entered for the scanList argument, will pop up a dialog
% to get the parameters.
%
% Essentially an evil function that should be used only for fun. It's
% main use is for standardizing the sampling resolution across different
% sessions where the true voxel size changed for some reason.
%
% INPUTS:
% vw : mrVista view structure : default selectedInplane
% scanList : scans to filter : default all
% kernelSize_mm : FWHM of the kernel in mm : default 3mm. If a 3 element
% vector then sizes should be [x,y,z]: NOTE! Z (down the slices) is first
% newTypeName: default 'Blurred'
% ARW: 051804: Wrote it.
%
% Example:
%   spatialBlurTSeries(INPLANE{1},[1:2],3,'Blurred3mmAverages');
% .. then recompute the coranal in the new data type "Blurred3mmAverages"
%
% ras: 2005: added a dialog to get params. Don't think it's completely
% evil.

mrGlobals;

if (~exist('scanList','var') || isempty(scanList)), scanList = er_selectScans(vw); end
if (~exist('kernelSize_mm','var') || isempty(kernelSize_mm)), kernelSize_mm = 5; end
if (~exist('newTypeName','var') || isempty(newTypeName)), newTypeName='Blurred'; end

if isequal(lower(scanList), 'dialog')
    % pop up a dialog (or two)
    scanList = er_selectScans(vw, 'Select Scans to Spatially Blur');

    switch viewGet(vw, 'viewType')
        case 'Inplane'
            prompt={'Scan List:', 'Gaussian Kernel Size (mm):', ...
                'Name of New Data Type'};
            name='Input for Peaks function';
            defaults = {num2str(scanList), '3', 'Blurred3mm'};
        case 'Gray'
            prompt={'Scan List:', 'NumIterations Sigma:', ...
                'Name of New Data Type'};
            name='Input for Peaks function';
            defaults = {num2str(scanList), '9  0.5', 'dhkGraySmooth'};
    end
    resp = inputdlg(prompt,name, 1, defaults);
    scanList        = str2num(resp{1}); %#ok<*ST2NM>
    kernelSize_mm   = str2num(resp{2});
    newTypeName     = resp{3};
end

if ~existDataType(newTypeName), addDataType(newTypeName); end

origDataType=vw.curDataType;
srcDtName = dataTYPES(origDataType).name;

% if ~strcmp(vw.viewType,'Inplane')
%     error('This function only operates on inplane data');
% end

verbose = prefsVerboseCheck;

switch viewGet(vw, 'viewType')

    case 'Inplane'

        % Work out the size of the kernel we need in voxels
        % First find the effective voxel size:
        effectiveVoxSize_mm=mrSESSION.functionals(1).effectiveResolution(:);
        kernelSize_mm=kernelSize_mm(:);
        if (length(kernelSize_mm)~=3)
            kernelSize_mm=ones(3,1)*kernelSize_mm(1);
        end
        kernelRatio=kernelSize_mm./effectiveVoxSize_mm;
        support=kernelRatio*2;

        % Because we are going to have 'z' as the first dimension, rotate both
        % these vectors
        kernelRatio=kernelRatio([3 1 2]);
        support=support([3 1 2]);
        fprintf('\nKernel ratio=%d %d %d\n',kernelRatio(1),kernelRatio(2),kernelRatio(3));

        spatialFilter=gauss3d(support,kernelRatio);
        spatialFilter=spatialFilter./sum(spatialFilter(:));

        hiddenView = initHiddenInplane;

        % Pre-define a large matrix to hold all the data for a single scan
        % This will be nVoxels*nVoxels*nSlices*nTRs
        scanParams= dataTYPES(origDataType).scanParams(1);
        nx=scanParams.cropSize(1);
        ny=scanParams.cropSize(2);
        nSlices=length(scanParams.slices);

    case 'Gray'
        hiddenView = initHiddenGray;
end

hiddenView = selectDataType(hiddenView,existDataType(newTypeName));

for thisScan=1:length(scanList)

    % Initialize the new scan in dataTYPES.
    newScanNum = viewGet(hiddenView, 'num scans')+1;
    initScan(vw, newTypeName, newScanNum, {srcDtName scanList(thisScan)});

    switch viewGet(vw, 'viewType')
        case 'Inplane'
            nTR = viewGet(vw, 'NumFrames', scanList(thisScan));
            if verbose, disp('Allocating large matrix'); end
            dataArray=zeros(nTR,nx,ny,nSlices);

            % Get the tSeries directory for this dataType
            % (make the directory if it doesn't already exist).
            tseriesdir = tSeriesDir(hiddenView);

            % Make the Scan subdirectory for the new tSeries (if it doesn't exist)
            scandir = fullfile(tseriesdir,['Scan',num2str(newScanNum)]);
            if ~exist(scandir,'dir')
                mkdir(tseriesdir,['Scan',num2str(newScanNum)]);
            end

            % Load in the full data set for each scan
            % (nVoxelsx*nVoxelsy*nSlices*nTR)
            % Then loop through the TRs performing a spatial blur (convn) on the volume
            % at each time point.
            disp('Loading');
            for thisSlice=1:nSlices
                thisSliceData=loadtSeries(vw,scanList(thisScan),thisSlice); % Data come in as nTR*nVoxels
                dataArray(:,:,:,thisSlice)=reshape(thisSliceData,nTR,nx,ny);
            end

            % Now loop over TRs doing the convolution
            disp('Convolving');
            for thisTR=1:nTR
                dataArray(thisTR,:,:,:)=convn(dataArray(thisTR,:,:,:),spatialFilter,'same');
            end

            % Now loop over slices again, saving out the data
            disp('Saving');

            dataArray = permute(dataArray, [1 2 4 3]);

            savetSeries(dataArray,hiddenView,thisScan);

        case 'Gray'
            tSeries = loadtSeries(vw,scanList(thisScan));
            iterlambda = kernelSize_mm;
            tSeries = dhkGraySmooth(vw,tSeries,iterlambda);
            %This should be fine since we this is a 'gray' view
            savetSeries(tSeries, hiddenView,thisScan,1);
    end

    fprintf('\nDone scan %d\n',thisScan);

end % next scan

% Loop through the open views, switch their curDataType appropriately,
% and update the dataType popups
INPLANE = resetDataTypes(INPLANE);
VOLUME  = resetDataTypes(VOLUME);
FLAT    = resetDataTypes(FLAT);

% in the provided view, select the new data type
vw = selectDataType(vw, newTypeName);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% shear_poisson2bulk
% Below is a demonstration of the features of the |shear_poisson2bulk| function

%%
clear; close all; clc;

%% Syntax
% |k=shear_poisson2bulk(mu,v);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function a = clement2_inverse ( n, x, y )

%*****************************************************************************80
%
%% CLEMENT2_INVERSE returns the inverse of the CLEMENT2 matrix.
%
%  Example:
%
%    N = 6, X and Y arbitrary:
%
%     0                1/Y1 0         -X2/(Y1*Y3) 0   X2*X4/(Y1*Y3*Y5)
%     1/X1             0    0          0          0    0
%     0                0    0          1/Y3       0   -X4/(Y3*Y5)
%    -Y2/(X1*X3)       0    1/X3       0          0    0
%     0                0    0          0          0    1/Y5
%     Y2*Y4/(X1*X3*X5) 0   -Y4/(X3*X5) 0          1/X5 0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Clement,
%    A class of triple-diagonal matrices for test purposes,
%    SIAM Review,
%    Volume 1, 1959, pages 50-52.
%
%  Parameters:
%
%    Input, integer N, the order of A.  N must not be odd%
%
%    Input, real X(N-1), Y(N-1), the first super and
%    subdiagonals of the matrix A.  None of the entries
%    of X or Y may be zero.
%
%    Output, real A(N,N), the matrix.
%
  if ( mod ( n, 2 ) == 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CLEMENT2_INVERSE - Fatal error!\n' );
    fprintf ( 1,'  The matrix is singular for odd N.\n' );
    error ( 'CLEMENT2_INVERSE - Fatal error!' );
  end

  for i = 1 : n-1

    if ( x(i) == 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CLEMENT2_INVERSE - Fatal error!\n' );
      fprintf ( 1, '  The matrix is singular\n' );
      fprintf ( 1, '  X(I) = 0 for I = %d\n', i );
      error ( 'CLEMENT2_INVERSE - Fatal error!' );
    elseif ( y(i) == 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CLEMENT2_INVERSE - Fatal error!\n' );
      fprintf ( 1, '  The matrix is singular\n' );
      fprintf ( 1, '  Y(I) = 0 for I = %d\n', i );
      error ( 'CLEMENT2_INVERSE - Fatal error!' );
    end

  end

  a = zeros ( n, n );

  for i = 1 : n

    if ( mod ( i, 2 ) == 1 )

      for j = i : 2 : n-1

        if ( j == i )
          prod1 = 1.0 / y(j);
          prod2 = 1.0 / x(j);
        else
          prod1 = - prod1 * x(j-1) / y(j);
          prod2 = - prod2 * y(j-1) / x(j);
        end

        a(i,j+1) = prod1;
        a(j+1,i) = prod2;

      end

    end

  end

  return
end


text:
function c = correlation_pentaspherical ( n, rho, rho0 )

%*****************************************************************************80
%
%% CORRELATION_PENTASPHERICAL evaluates the pentaspherical correlation function.
%
%  Discussion:
%
%    This correlation is based on the volume of overlap of two spheres
%    of radius RHO0 and separation RHO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Petter Abrahamsen,
%    A Review of Gaussian Random Fields and Correlation Functions,
%    Norwegian Computing Center, 1997.
%
%  Parameters:
%
%    Input, integer N, the number of arguments.
%
%    Input, real RHO(N,1), the arguments.
%
%    Input, real RHO0, the correlation length.
%
%    Output, real C(N,1), the correlations.
%
  rho = rho ( : );

  rhohat = min ( abs ( rho ) / rho0, 1.0 );

  c = 1.0 - 1.875 * rhohat + 1.25 * rhohat.^3 - 0.375 * rhohat.^5;

  return
end


text:
function lmi = lmi_input(lmi, umax, phi)
%LMI_INPUT

% reference:
%	K. Tanaka, H. O. Wang
%	Fuzzy Control Systems Design and Analysis, (2002)
%	page 66,69 (theorem 11,13)

% TODO: known x(0), output constraint

R = size(lmi.A, 1);
X = lmi.X;
M = lmi.M;

% constraints on the control value

% phi^2 I < X
lmi.F = lmi.F + set(phi^2 * eye(lmi.n) < X, 'phi^2 I < X');

% [X, Mr'; Mr, mu^2 I] > 0
for r = 1:R
	lmi.F = lmi.F + set([X M{r}'; M{r} umax^2*eye(lmi.m)] > 0, sprintf('type3 lmi %d', r));
end


text:
%% KEOLLOGG Problem
%
% KELLOGG solves a diffusion equation with jump coefficients with AFEM.
%
% KELLOGG solves the problem within maxN number of vertices. The
% input argument theta is a parameter used in the marking step.
%
% The KELLOGG command, if no input arguments, use maxN = 5e3 and theta = 0.5.
%
% EXAMPLE
%
%    Kellogg
%
% See also  crack, Lshape
%
% TODO: rewrite M-lint
%
% Created by Chen-Song Zhang. Modified by Long Chen.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

close all
%% Problem setting
% $$-\nabla\cdot(d\nabla u) = 0  \quad \Omega = (-1,1)\times (-1,1)$$
%
% $$u = g_D \quad \partial \Omega$$
%
% The diffusion constant is discontinous; see the figure below. We set a2  =
% 1; a1 = 161.4476387975881 and choose boundary condition g_D such that the
% exact solution is $z = r^{0.1}\mu(\theta)$ in the poloar coordinate, where
% the formula of mu can be found in exactu function.

[x,y] = meshgrid(-1:1:1,-1:1:1);
z = 0*x;
surf(x,y,z,'linewidth',2); view(2);
axis equal; axis tight;
text(0.5,0.5,'a1','FontSize',12,'FontWeight','bold');
text(-0.5,-0.5,'a1','FontSize',12,'FontWeight','bold');
text(-0.5,0.5,'a2','FontSize',12,'FontWeight','bold');
text(0.5,-0.5,'a2','FontSize',12,'FontWeight','bold');

%% Parameters
maxN = 3e3;     theta = 0.2;    maxIt = 1000;
N = zeros(maxIt,1);     uIuhErr = zeros(maxIt,1);
errH1 = zeros(maxIt,1);

%%  Generate an initial mesh
[node,elem] = squaremesh([-1 1 -1 1], 0.25);
bdFlag = setboundary(node,elem,'Dirichlet');

%% Set up PDE data
pde = Kelloggdata;

%%  Adaptive Finite Element Method
% *SOLVE* -> *ESTIMATE* -> *MARK* -> *REFINE*

for k = 1:maxIt
    % Step 1: SOLVE
    [u,Du,eqn] = Poisson(node,elem,pde);
    % Plot mesh and solution
    figure(1);  showresult(node,elem,u,[27,26]);
    % Step 2: ESTIMATE
   eta = estimaterecovery(node,elem,u);            % recovery type
%     eta = estimateresidual(node,elem,u,pde);    % residual type
    % Record error and number of vertices
    uI = pde.exactu(node);
    uIuhErr(k) = sqrt((uI-u)'*eqn.A*(uI-u));
    errH1(k) = getH1error(node,elem,@pde.Du,u);
    N(k) = size(node,1);
    if (N(k)>maxN), break; end
    % Step 3: MARK
    markedElem = mark(elem,eta,theta);
    % Step 4: REFINE
    [node,elem,bdFlag] = bisect(node,elem,markedElem,bdFlag);
end

%%  Plot convergent rates in energy norm
N = N(1:k);
uIuhErr = uIuhErr(1:k);
errH1 = errH1(1:k);
figure(2)
showrate2(N,uIuhErr,10,'-*','||Du_I-Du_h||',...
          N,errH1,10,'k-.','||Du-Du_h||');

repo: Aero-matlab


text:
function varargout = column_vector(varargin)
%COLUMN_VECTOR returns the given vectors as column vectors
%
%   Usage: [x1,x2,...] = column_vector(x1,x2,...)
%
%   Input parameters:
%       x1,x2,...  - one or more vectors
%
%   Output parameters:
%       x1,x2,...  - input vectors as column vectors
%
%   See also: row_vector

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ===================================
isargvector(varargin{:});

%% ===== Computation =====================================================
for ii = 1:nargin
    if size(varargin{ii},2)>size(varargin{ii},1)
        varargout{ii} = varargin{ii}';
    else
        varargout{ii} = varargin{ii};
    end
end

author: sfstoolbox

repo: sfs-matlab


text:
function prob_test043 ( )

%*****************************************************************************80
%
%% TEST043 tests DERANGED_MEAN, DERANGED_VARIANCE, DERANGED_SAMPLE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST043\n' );
  fprintf ( 1, '  For the Deranged PDF:\n' );
  fprintf ( 1, '  DERANGED_MEAN computes the mean.\n' );
  fprintf ( 1, '  DERANGED_VARIANCE computes the variance.\n' );
  fprintf ( 1, '  DERANGED_SAMPLE samples.\n' );

  a = 7;

  check = deranged_check ( a );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST043 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = deranged_mean ( a );
  variance = deranged_variance ( a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %14f\n', a );
  fprintf ( 1, '  PDF mean =        %14f\n', mean );
  fprintf ( 1, '  PDF variance =    %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = deranged_sample ( a, seed );
  end

  mean = i4vec_mean ( nsample, x );
  variance = i4vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %6d\n', xmax );
  fprintf ( 1, '  Sample minimum =  %6d\n', xmin );

  return
end


text:
function [Maxima,MaxPos,Minima,MinPos]=MinimaMaxima3D(Input,Robust,LookInBoundaries,numbermax,numbermin)
% V 1.0 Dec 13, 07
% Author Sam Pichardo.
% This  function finds the local minima and maxima in a 3D Cartesian data.
% It's assumed that the data is uniformly distributed.
% The minima and maxima are calculated using a multi-directional derivation.
%
% Use:
%
%  [Maxima,MaxPos,Minima,MinPos]=MinimaMaxima3D(Input,[Robust],[LookInBoundaries],[numbermax],[numbermin])
%
% where Input is the 3D data and Robust (optional and with a default value
% of 1) indicates if the multi-directional derivation should include the
% diagonal derivations.
%
% Input has to have a size larger or equal than [3 x 3 x 3]
%
% If Robust=1, the total number of derivations taken into account are 26: 6
% for all surrounding elements colliding each of the faces of the unit cube;
% 10 for all the surrounding elements in diagonal.
%
% If Robust =0, then only the 6 elements of the colliding faces are considered
%
% The function returns in Maxima and MaxPos, respectively,
% the values (numbermax) and subindexes (numbermax x 3) of local maxima
% and position in Input. Maxima (and the subindexes) are sorted in
% descending order.
% Similar situation for Minima and MinimaPos witn a numbermin elements but
% with the execption of being sorted in ascending order.
%
% IMPORTANT: if numbermin or numbermax are not specified, ALL the minima
% or maxima will be returned. This can be a useless for highly
% oscillating data
%
% LookInBoundaries (default value of 0) specifies if a search of the minima/maxima should be
% done in the boundaries of the matrix. This situation depends on the
% the desire application. When it is not activated, the algorithm WILL NOT
% FIND ANY MINIMA/MAXIMA on the 6 layers of the boundaries.
% When it is activated, the finding minima and maxima on the boundaries is done by
% replicating the extra layer as the layer 2 (or layer N-1, depending of the boundary)
% By example (and using a 2D matrix for simplicity reasons):
% For the matrix
% [ 4 1 3 7
%   5 7 8 8
%   9 9 9 9
%   5 6 7 9]
%
% the calculation of the partial derivate following the -x direction will be done by substrascting
% [ 5 7 8 8
%   4 1 3 7
%   5 7 8 8
%   9 9 9 9]
% to the input. And so on for the other dimensions.
% Like this, the value "1" at the coordinate (1,2) will be detected as a
% minima. Same situation for the value "5" at the coordinate (4,1)

if nargin <1
    test=load('temp.mat');
    pf=test.uresTot(test.EvalLims(2,1):test.EvalLims(2,2));
    pf=reshape(pf,length(test.EvalCoord{2}.Ry),length(test.EvalCoord{2}.Rx),length(test.EvalCoord{2}.Rz));
    Input = abs(pf)*1.5e6;
    clear test;
    clear pf;
    Robust =1;
end

Asize=size(Input);

if length(Asize)<3
    error('MinimaMaxima3D can only works with 3D matrices ');
end

if (Asize(1)<3 || Asize(2)<3 || Asize(3)<3)
    error('MinimaMaxima3D can only works with matrices with dimensions equal or larger to [3x3x3]');
end

if ~isreal(Input)
    warning('ATTENTION, complex values detected!!, using abs(Input)');
    Input=abs(Input);
end

if ~exist('Robust','var')
    Robust=1;
end

if ~exist('LookInBoundaries','var')
    LookInBoundaries=0;
end

if ~exist('numbermax','var')
    numbermax=0;
end

if ~exist('numbermin','var')
    numbermin=0;
end

[xx_base,yy_base,zz_base]=ndgrid(1:Asize(1),1:Asize(2),1:Asize(3));

IndBase=sub2ind(Asize,xx_base(:),yy_base(:),zz_base(:));

if Robust ~= 0
    Numbder_dd=26;
else
    Numbder_dd=6;
end

if LookInBoundaries==0
    lx=1:Asize(1);
    lx_p1=[2:Asize(1),Asize(1)];
    lx_m1=[1,1:Asize(1)-1];
    ly=1:Asize(2);
    ly_p1=[2:Asize(2),Asize(2)];
    ly_m1=[1,1:Asize(2)-1];
    lz=1:Asize(3);
    lz_p1=[2:Asize(3),Asize(3)];
    lz_m1=[1,1:Asize(3)-1];
else
    lx=1:Asize(1);
    lx_p1=[2:Asize(1),Asize(1)-1]; %We replicate the layer N-1 as the layer N+1
    lx_m1=[2,1:Asize(1)-1]; %We replicate the layer 2 as the layer -1
    ly=1:Asize(2);
    ly_p1=[2:Asize(2),Asize(2)-1]; %We replicate the layer N-1 as the layer N+1
    ly_m1=[2,1:Asize(2)-1]; %We replicate the layer 2 as the layer -1
    lz=1:Asize(3);
    lz_p1=[2:Asize(3),Asize(3)-1]; %We replicate the layer N-1 as the layer N+1
    lz_m1=[2,1:Asize(3)-1];%We replicate the layer 2 as the layer -1
end

for n_dd=1:Numbder_dd
    switch n_dd
        case 1
            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1)
            [xx,yy,zz]=ndgrid(lx_p1,ly,lz);

        case 2
            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1)
            [xx,yy,zz]=ndgrid(lx_m1,ly,lz);

        case 3
            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(y)-elem(y+1)
            [xx,yy,zz]=ndgrid(lx,ly_p1,lz);

        case 4
            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(y)-elem(y-1)
            [xx,yy,zz]=ndgrid(lx,ly_m1,lz);

        case 5
            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(z)-elem(z+1)
            [xx,yy,zz]=ndgrid(lx,ly,lz_p1);

         case 6
            %%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(z)-elem(z-1)
            [xx,yy,zz]=ndgrid(lx,ly,lz_m1);
        case 7
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y+1)
            [xx,yy,zz]=ndgrid(lx_p1,ly_p1,lz);
        case 8
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y-1)
            [xx,yy,zz]=ndgrid(lx_p1,ly_m1,lz);
        case 9
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y-1)
            [xx,yy,zz]=ndgrid(lx_m1,ly_m1,lz);
        case 10
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y+1)
            [xx,yy,zz]=ndgrid(lx_m1,ly_p1,lz);
        case 11
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,z+1)
            [xx,yy,zz]=ndgrid(lx_p1,ly,lz_p1);
        case 12
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,z-1)
            [xx,yy,zz]=ndgrid(lx_p1,ly,lz_m1);
        case 13
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,z-1)
            [xx,yy,zz]=ndgrid(lx_m1,ly,lz_m1);
        case 14
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,z+1)
            [xx,yy,zz]=ndgrid(lx_m1,ly,lz_p1);
        case 15
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(y+1,z+1)
            [xx,yy,zz]=ndgrid(lx,ly_p1,lz_p1);
        case 16
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(y+1,z-1)
            [xx,yy,zz]=ndgrid(lx,ly_p1,lz_m1);
        case 17
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(y-1,z-1)
            [xx,yy,zz]=ndgrid(lx,ly_m1,lz_m1);
        case 18
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(y-1,z+1)
            [xx,yy,zz]=ndgrid(lx,ly_m1,lz_p1);
         case 19
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y+1,z+1)
            [xx,yy,zz]=ndgrid(lx_p1,ly_p1,lz_p1);
         case 20
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y+1,z-1)
            [xx,yy,zz]=ndgrid(lx_p1,ly_p1,lz_m1);
         case 21
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y-1,z+1)
            [xx,yy,zz]=ndgrid(lx_p1,ly_m1,lz_p1);
         case 22
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y-1,z-1)
            [xx,yy,zz]=ndgrid(lx_p1,ly_m1,lz_m1);
         case 23
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y+1,z+1)
            [xx,yy,zz]=ndgrid(lx_m1,ly_p1,lz_p1);
         case 24
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y+1,z-1)
            [xx,yy,zz]=ndgrid(lx_m1,ly_p1,lz_m1);
         case 25
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y-1,z+1)
            [xx,yy,zz]=ndgrid(lx_m1,ly_m1,lz_p1);
         case 26
           %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y-1,z-1)
            [xx,yy,zz]=ndgrid(lx_m1,ly_m1,lz_m1);

    end

    Ind_dd=sub2ind(Asize,xx(:),yy(:),zz(:));

    part_deriv = Input(IndBase)-Input(Ind_dd);

    if n_dd >1
        MatMinMax= (sign_Prev_deriv==sign(part_deriv)).*MatMinMax;
    else
        MatMinMax=sign(part_deriv);
    end

    sign_Prev_deriv=sign(part_deriv);
end

%Well , now the easy part, all values MatMinMax ==1 are local maximum and
%the values MatMinMax ==-1 are minimun

AllMaxima=find(MatMinMax==1);
AllMinima=find(MatMinMax==-1);

if numbermax ==0
    nmax=length(AllMaxima);
else
    nmax=numbermax;
end
nmax=min([nmax,length(AllMaxima)]);
smax=1:nmax;

if numbermin ==0
    nmin=length(AllMinima);
else
    nmin=numbermin;
end

nmin=min([nmin,length(AllMinima)]);

smin=1:nmin;

[Maxima,IndMax]=sort(Input(AllMaxima),'descend');
Maxima=Maxima(smax);
IndMax=AllMaxima(IndMax(smax));

MaxPos=zeros(nmax,3);
[MaxPos(:,1),MaxPos(:,2),MaxPos(:,3)]=ind2sub(Asize,IndMax);

[Minima,IndMin]=sort(Input(AllMinima));
Minima=Minima(smin);
IndMin=AllMinima(IndMin(smin));

MinPos=zeros(nmin,3);
[MinPos(:,1),MinPos(:,2),MinPos(:,3)]=ind2sub(Asize,IndMin);


text:
%TRANSL2 SE(2) translational homogeneous transform
%
% Create a translational SE(2) matrix::
%
% T = TRANSL2(X, Y) is an SE(2) homogeneous transform (3x3) representing a
% pure translation.
%
% T = TRANSL2(P) is a homogeneous transform representing a translation or
% point P=[X,Y]. P (Mx2) represents a sequence and T (3x3xM) is a
% sequence of homogenous transforms such that T(:,:,i) corresponds to the
% i'th row of P.
%
% Extract the translational part of an SE(2) matrix::
%
% P = TRANSL2(T) is the translational part of a homogeneous transform as a
% 2-element column vector.  T (3x3xM) is a homogeneous transform
% sequence and the rows of P (Mx2) are the translational component of the
% corresponding transform in the sequence.
%
% Notes::
% - Somewhat unusually, this function performs a function and its inverse.  An
%   historical anomaly.
%
% See also SE2.t, ROT2, ISHOMOG2, TRPLOT2, TRANSL.

%## 2d homogeneous translation

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

function T = transl2(x, y)
    if nargin == 1
        if ishomog2(x)
            if ndims(x) == 3
                % transl(T)  -> P, trajectory case
                T = squeeze(x(1:2,3,:))';
            else
                % transl(T)  -> P
                T = x(1:2,3);
            end
        elseif all(size(x) == [3 3])
            T = x(1:2,3);
        elseif length(x) == 2
            % transl(P) -> T
            t = x(:);
            T =    [eye(2)          t(:);
                0   0   1];
        else
            % transl(P) -> T, trajectory case
            n = size(x,1);
            T = repmat(eye(3,3), [1 1 n]);
            T(1:2,3,:) = x';
        end
    elseif nargin == 2
        % transl(x,y) -> T
        t = [x; y];
        T =    [ eye(2) t; 0 0 1];
    end
end

% one less function to upload for Cody/LTI assessments

function h = ishomog2(tr, rtest)
    d = size(tr);
    if ndims(tr) >= 2
        h =  all(d(1:2) == [3 3]);

        if h && nargin > 1
            h = abs(det(tr(1:2,1:2)) - 1) < eps;
        end
    else
        h = false;
    end
end

author: petercorke

repo: spatialmath-matlab


text:
function show_curve(b)
% function show_curve(b)

if ~strcmpi(class(b),'bezier')
  error('Expecting a bezier object');
end

t=0:0.01:1;
t2=t.^2;
t3=t2.*t;
x=b.ax*t3+b.bx*t2+b.cx*t+b.x0;
y=b.ay*t3+b.by*t2+b.cy*t+b.y0;

plot(x,y,'b',[b.x0],[b.y0],'r*',...
     [b.x1],[b.y1],'ro',...
     [b.x0,b.x1],[b.y0,b.y1],'r',...
     [b.x2],[b.y2],'ro',...
     [b.x3],[b.y3],'r*',...
     [b.x2,b.x3],[b.y2,b.y3],'r');

text([b.x0],[b.y0],'(x_0,y_0)');
text([b.x1],[b.y1],'(x_1,y_1)');
text([b.x2],[b.y2],'(x_2,y_2)');
text([b.x3],[b.y3],'(x_3,y_3)');


text:
function [vec,rotMat]=TEME2GCRS(x,Jul1,Jul2,deltaTTUT1,xpyp,dXdY,LOD)
%%TEME2GCRS Convert from the True Equator Mean Equinox (TEME) of date
%           coordinate system to the Geocentric Celestrial Reference System
%           (GCRS), a type of Earth-Centered Inertial (ECI) coordinate
%           system. The TEME system is non-standard and is generally
%           only used in the Specialized General Perturbations 4 (SGP4)
%           orbit propagation algorithm.
%
%INPUTS: x The NXnumVec collection of vectors in TEME coordinates to
%          convert. N can be 3, or 6. If the vectors are 3D, then they are
%          position. 6D vectors are assumed to be position and velocity,
%          whereby the angular velocity of the Earth's rotation is taken
%          into account using a non-relativistic formula.
% Jul1, Jul2 Two parts of a Julian date given in terrestrial time (TT).
%          The units of the date are days. The full date is the sum of both
%          terms. The date is broken into two parts to provide more bits of
%          precision. It does not matter how the date is split.
% deltaTTUT1 An optional parameter specifying the difference between TT and
%          UT1 in seconds. This information can be obtained from
%          http://www.iers.org/nn_11474/IERS/EN/DataProducts/EarthOrientationData/eop.html?__nnn=true
%          or
%          http://www.usno.navy.mil/USNO/earth-orientation/eo-products
%          If this parameter is omitted or if an empty matrix is passed,
%          then the value provided by the function getEOP will be used
%          instead.
%     xpyp xpyp=[xp;yp] are the polar motion coordinates in radians
%          including the effects of tides and librations. If this parameter
%          is omitted or if an empty matrix is passed, the value from the
%          function getEOP will be used.
%     dXdY dXdY=[dX;dY] are the celestial pole offsets with respect to the
%          IAU 2006/2000A precession/nutation model in radians If this
%          parameter is omitted or if an empty matrix is passed, the value
%          from the function getEOP will be used.
%      LOD The difference between the length of the day using terrestrial
%          time, international atomic time, or UTC without leap seconds and
%          the length of the day in UT1. This is an instantaneous parameter
%          (in seconds) proportional to the rotation rate of the Earth.
%          This is only needed if more than just position components are
%          being converted.
%
%OUTPUTS: vec A 3XN or 6XN matrix of vectors converted from TEME
%             coordinates to GCRS coordinates.
%      rotMat The 3X3 rotation matrix used for the conversion of the
%             positions. That is, vec(1:3)=rotMat*x(1:3)
%
%This function just calls the functions TEME2ITRS and ITRS2GCRS, reusing
%the same Earth orientation parameters across calls.
%
%Different celestial coordinate systems are compared in [1].
%
%REFERENCES:
%[1] D. F. Crouse, "An Overview of Major Terrestrial, Celestial, and
%    Temporal Coordinate Systems for Target Tracking," Formal Report, Naval
%    Research Laboratory, no. NRL/FR/5344--16-10,279, 10 Aug. 2016, 173
%    pages.
%
%January 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%If more EOPs are needed.
if(nargin<7&&size(x,1)==6||nargin<6||isempty(dXdY)||isempty(xpyp)||isempty(deltaTTUT1))
    [JulUTC1,JulUTC2]=TT2UTC(Jul1,Jul2);
    [xpypget,dXdYget,~,deltaTTUT1get,LODget]=getEOP(JulUTC1,JulUTC2);
end

if(nargin<7||isempty(LOD))
   LOD=LODget;
end

if(nargin<6||isempty(dXdY))
    dXdY=dXdYget;
end

if(nargin<5||isempty(xpyp))
    xpyp=xpypget;
end

if(nargin<4||isempty(deltaTTUT1))
    deltaTTUT1=deltaTTUT1get;
end

[x,rotMat1]=TEME2ITRS(x,Jul1,Jul2,deltaTTUT1,xpyp,LOD);
[vec,rotMat2]=ITRS2GCRS(x,Jul1,Jul2,deltaTTUT1,xpyp,dXdY,LOD);

%The combined rotation matrix.
rotMat=rotMat2*rotMat1;
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef GeneralizedHoughBallard < handle
    %GENERALIZEDHOUGHBALLARD  Generalized Hough transform
    %
    % Finds arbitrary template in the grayscale image using
    % Generalized Hough Transform.
    %
    % Detects position only without translation and rotation.
    %
    % ## References
    % > Ballard, D.H. (1981). Generalizing the Hough transform to detect
    % > arbitrary shapes. Pattern Recognition 13 (2): 111-122.
    %
    % See also: cv.GeneralizedHoughGuil
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % Canny low threshold. default 50
        CannyLowThresh
        % Canny high threshold. default 100
        CannyHighThresh
        % Minimum distance between the centers of the detected objects.
        % default 1.0
        MinDist
        % Inverse ratio of the accumulator resolution to the image resolution.
        % default 1.0
        Dp
        % Maximal size of inner buffers. default 1000
        MaxBufferSize

        % R-Table levels. default 360
        Levels
        % The accumulator threshold for the template centers at the detection
        % stage. The smaller it is, the more false positions may be detected.
        % default 100
        VotesThreshold
    end

    %% GeneralizedHoughBallard
    methods
        function this = GeneralizedHoughBallard()
            %GENERALIZEDHOUGHBALLARD  Constructor
            %
            %     obj = cv.GeneralizedHoughBallard()
            %
            % See also: cv.GeneralizedHoughBallard.detect
            %
            this.id = GeneralizedHoughBallard_(0, 'new');
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.GeneralizedHoughBallard
            %
            if isempty(this.id), return; end
            GeneralizedHoughBallard_(this.id, 'delete');
        end
    end

    %% Algorithm
    methods (Hidden)
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.GeneralizedHoughBallard.empty,
            %  cv.GeneralizedHoughBallard.load
            %
            GeneralizedHoughBallard_(this.id, 'clear');
        end

        function load(this, filename, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.GeneralizedHoughBallard.save
            %
            GeneralizedHoughBallard_(this.id, 'load', filename, varargin{:});
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.GeneralizedHoughBallard.load
            %
            GeneralizedHoughBallard_(this.id, 'save', filename);
        end

        function b = empty(this)
            %EMPTY  Returns true if the algorithm is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.GeneralizedHoughBallard.clear,
            %  cv.GeneralizedHoughBallard.load
            %
            b = GeneralizedHoughBallard_(this.id, 'empty');
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.GeneralizedHoughBallard.save,
            %  cv.GeneralizedHoughBallard.load
            %
            name = GeneralizedHoughBallard_(this.id, 'getDefaultName');
        end
    end

    %% GeneralizedHough
    methods
        function [positions,votes] = detect(this, varargin)
            %DETECT  Find template on image
            %
            %     [positions,votes] = hough.detect(image)
            %     [positions,votes] = hough.detect(edges, dx, dy)
            %
            % ## Input
            % * __image__ input image, 8-bit 1-channel image
            % * __edges__ image edges
            % * __dx__ image x-derivative of the same size as `edges` and
            %   single-precision floating type.
            % * __dy__ image y-derivate of the same size and type as `dx`.
            %
            % ## Output
            % * __positions__ Cell array of 4-element vectors, each of the
            %   form: `[posx, posy, scale, angle]`
            % * __votes__ Cell array of 3-element vectors, of the same length
            %   as `positions`.
            %
            [positions,votes] = GeneralizedHoughBallard_(this.id, 'detect', varargin{:});
        end

        function setTemplate(this, varargin)
            %SETTEMPLATE  Set template to search
            %
            %     hough.setTemplate(templ)
            %     hough.setTemplate(edges, dx, dy)
            %     hough.setTemplate(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __templ__ template, 8-bit 1-channel image
            % * __edges__ template edges
            % * __dx__ template x-derivative of the same size as `edges` and
            %   single-precision floating type.
            % * __dy__ template y-derivate of the same size and type as `dx`.
            %
            % ## Options
            % * __Center__ Template center `[x,y]`. The default `[-1,-1]`
            %   will use `[size(templ,2) size(templ,1)]./2` as center.
            %
            % In the first variant of the function (with the `templ` input),
            % the `edges` are internally calculated using the cv.Canny filter,
            % and the derivatives `dx` and `dy` using cv.Sobel operator.
            %
            GeneralizedHoughBallard_(this.id, 'setTemplate', varargin{:});
        end
    end

    %% Getters/Setters
    methods
        function value = get.CannyHighThresh(this)
            value = GeneralizedHoughBallard_(this.id, 'get', 'CannyHighThresh');
        end
        function set.CannyHighThresh(this, value)
            GeneralizedHoughBallard_(this.id, 'set', 'CannyHighThresh', value);
        end

        function value = get.CannyLowThresh(this)
            value = GeneralizedHoughBallard_(this.id, 'get', 'CannyLowThresh');
        end
        function set.CannyLowThresh(this, value)
            GeneralizedHoughBallard_(this.id, 'set', 'CannyLowThresh', value);
        end

        function value = get.Dp(this)
            value = GeneralizedHoughBallard_(this.id, 'get', 'Dp');
        end
        function set.Dp(this, value)
            GeneralizedHoughBallard_(this.id, 'set', 'Dp', value);
        end

        function value = get.MaxBufferSize(this)
            value = GeneralizedHoughBallard_(this.id, 'get', 'MaxBufferSize');
        end
        function set.MaxBufferSize(this, value)
            GeneralizedHoughBallard_(this.id, 'set', 'MaxBufferSize', value);
        end

        function value = get.MinDist(this)
            value = GeneralizedHoughBallard_(this.id, 'get', 'MinDist');
        end
        function set.MinDist(this, value)
            GeneralizedHoughBallard_(this.id, 'set', 'MinDist', value);
        end

        function value = get.Levels(this)
            value = GeneralizedHoughBallard_(this.id, 'get', 'Levels');
        end
        function set.Levels(this, value)
            GeneralizedHoughBallard_(this.id, 'set', 'Levels', value);
        end

        function value = get.VotesThreshold(this)
            value = GeneralizedHoughBallard_(this.id, 'get', 'VotesThreshold');
        end
        function set.VotesThreshold(this, value)
            GeneralizedHoughBallard_(this.id, 'set', 'VotesThreshold', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function submit()
  addpath('./lib');

  conf.assignmentSlug = 'anomaly-detection-and-recommender-systems';
  conf.itemName = 'Anomaly Detection and Recommender Systems';
  conf.partArrays = { ...
    { ...
      '1', ...
      { 'estimateGaussian.m' }, ...
      'Estimate Gaussian Parameters', ...
    }, ...
    { ...
      '2', ...
      { 'selectThreshold.m' }, ...
      'Select Threshold', ...
    }, ...
    { ...
      '3', ...
      { 'cofiCostFunc.m' }, ...
      'Collaborative Filtering Cost', ...
    }, ...
    { ...
      '4', ...
      { 'cofiCostFunc.m' }, ...
      'Collaborative Filtering Gradient', ...
    }, ...
    { ...
      '5', ...
      { 'cofiCostFunc.m' }, ...
      'Regularized Cost', ...
    }, ...
    { ...
      '6', ...
      { 'cofiCostFunc.m' }, ...
      'Regularized Gradient', ...
    }, ...
  };
  conf.output = @output;

  submitWithConfiguration(conf);
end

function out = output(partId, auxstring)
  % Random Test Cases
  n_u = 3; n_m = 4; n = 5;
  X = reshape(sin(1:n_m*n), n_m, n);
  Theta = reshape(cos(1:n_u*n), n_u, n);
  Y = reshape(sin(1:2:2*n_m*n_u), n_m, n_u);
  R = Y > 0.5;
  pval = [abs(Y(:)) ; 0.001; 1];
  yval = [R(:) ; 1; 0];
  params = [X(:); Theta(:)];
  if partId == '1'
    [mu sigma2] = estimateGaussian(X);
    out = sprintf('%0.5f ', [mu(:); sigma2(:)]);
  elseif partId == '2'
    [bestEpsilon bestF1] = selectThreshold(yval, pval);
    out = sprintf('%0.5f ', [bestEpsilon(:); bestF1(:)]);
  elseif partId == '3'
    [J] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                       n, 0);
    out = sprintf('%0.5f ', J(:));
  elseif partId == '4'
    [J, grad] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                             n, 0);
    out = sprintf('%0.5f ', grad(:));
  elseif partId == '5'
    [J] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                       n, 1.5);
    out = sprintf('%0.5f ', J(:));
  elseif partId == '6'
    [J, grad] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                             n, 1.5);
    out = sprintf('%0.5f ', grad(:));
  end
end

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford


text:
function data = loadubjson(fname,varargin)
%
% data=loadubjson(fname,opt)
%    or
% data=loadubjson(fname,'param1',value1,'param2',value2,...)
%
% parse a JSON (JavaScript Object Notation) file or string
%
% authors:Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
%            date: 2013/08/01
%
% $Id: loadubjson.m 417 2014-01-21 22:34:49Z fangq $
%
% input:
%      fname: input file name, if fname contains "{}" or "[]", fname
%             will be interpreted as a UBJSON string
%      opt: a struct to store parsing options, opt can be replaced by
%           a list of ('param',value) pairs. The param string is equivallent
%           to a field in opt.
%
% output:
%      dat: a cell array, where {...} blocks are converted into cell arrays,
%           and [...] are converted to arrays
%
% license:
%     BSD license, see LICENSE_BSD.txt files for details
%
% -- this function is part of jsonlab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

global pos inStr len  esc index_esc len_esc isoct arraytoken

if(regexp(fname,'[\{\}\]\[]','once'))
   string=fname;
elseif(exist(fname,'file'))
   fid = fopen(fname,'rb');
   string = fread(fid,inf,'uint8=>char')';
   fclose(fid);
else
   error('input file does not exist');
end

pos = 1; len = length(string); inStr = string;
isoct=exist('OCTAVE_VERSION');
arraytoken=find(inStr=='[' | inStr==']' | inStr=='"');
jstr=regexprep(inStr,'\\\\','  ');
escquote=regexp(jstr,'\\"');
arraytoken=sort([arraytoken escquote]);

% String delimiters and escape chars identified to improve speed:
esc = find(inStr=='"' | inStr=='\' ); % comparable to: regexp(inStr, '["\\]');
index_esc = 1; len_esc = length(esc);

opt=varargin2struct(varargin{:});
jsoncount=1;
while pos <= len
    switch(next_char)
        case '{'
            data{jsoncount} = parse_object(opt);
        case '['
            data{jsoncount} = parse_array(opt);
        otherwise
            error_pos('Outer level structure must be an object or an array');
    end
    jsoncount=jsoncount+1;
end % while

jsoncount=length(data);
if(jsoncount==1 && iscell(data))
    data=data{1};
end

if(~isempty(data))
      if(isstruct(data)) % data can be a struct array
          data=jstruct2array(data);
      elseif(iscell(data))
          data=jcell2array(data);
      end
end

%%
function newdata=parse_collection(id,data,obj)

if(jsoncount>0 && exist('data','var'))
    if(~iscell(data))
       newdata=cell(1);
       newdata{1}=data;
       data=newdata;
    end
end

%%
function newdata=jcell2array(data)
len=length(data);
newdata=data;
for i=1:len
      if(isstruct(data{i}))
          newdata{i}=jstruct2array(data{i});
      elseif(iscell(data{i}))
          newdata{i}=jcell2array(data{i});
      end
end

%%-------------------------------------------------------------------------
function newdata=jstruct2array(data)
fn=fieldnames(data);
newdata=data;
len=length(data);
for i=1:length(fn) % depth-first
    for j=1:len
        if(isstruct(getfield(data(j),fn{i})))
            newdata(j)=setfield(newdata(j),fn{i},jstruct2array(getfield(data(j),fn{i})));
        end
    end
end
if(~isempty(strmatch('x0x5F_ArrayType_',fn)) && ~isempty(strmatch('x0x5F_ArrayData_',fn)))
  newdata=cell(len,1);
  for j=1:len
    ndata=cast(data(j).x0x5F_ArrayData_,data(j).x0x5F_ArrayType_);
    iscpx=0;
    if(~isempty(strmatch('x0x5F_ArrayIsComplex_',fn)))
        if(data(j).x0x5F_ArrayIsComplex_)
           iscpx=1;
        end
    end
    if(~isempty(strmatch('x0x5F_ArrayIsSparse_',fn)))
        if(data(j).x0x5F_ArrayIsSparse_)
            if(~isempty(strmatch('x0x5F_ArraySize_',fn)))
                dim=double(data(j).x0x5F_ArraySize_);
                if(iscpx && size(ndata,2)==4-any(dim==1))
                    ndata(:,end-1)=complex(ndata(:,end-1),ndata(:,end));
                end
                if isempty(ndata)
                    % All-zeros sparse
                    ndata=sparse(dim(1),prod(dim(2:end)));
                elseif dim(1)==1
                    % Sparse row vector
                    ndata=sparse(1,ndata(:,1),ndata(:,2),dim(1),prod(dim(2:end)));
                elseif dim(2)==1
                    % Sparse column vector
                    ndata=sparse(ndata(:,1),1,ndata(:,2),dim(1),prod(dim(2:end)));
                else
                    % Generic sparse array.
                    ndata=sparse(ndata(:,1),ndata(:,2),ndata(:,3),dim(1),prod(dim(2:end)));
                end
            else
                if(iscpx && size(ndata,2)==4)
                    ndata(:,3)=complex(ndata(:,3),ndata(:,4));
                end
                ndata=sparse(ndata(:,1),ndata(:,2),ndata(:,3));
            end
        end
    elseif(~isempty(strmatch('x0x5F_ArraySize_',fn)))
        if(iscpx && size(ndata,2)==2)
             ndata=complex(ndata(:,1),ndata(:,2));
        end
        ndata=reshape(ndata(:),data(j).x0x5F_ArraySize_);
    end
    newdata{j}=ndata;
  end
  if(len==1)
      newdata=newdata{1};
  end
end

%%-------------------------------------------------------------------------
function object = parse_object(varargin)
    parse_char('{');
    object = [];
    type='';
    count=-1;
    if(next_char == '$')
        type=inStr(pos+1); % TODO
        pos=pos+2;
    end
    if(next_char == '#')
        pos=pos+1;
        count=double(parse_number());
    end
    if next_char ~= '}'
        num=0;
        while 1
            str = parseStr(varargin{:});
            if isempty(str)
                error_pos('Name of value at position %d cannot be empty');
            end
            %parse_char(':');
            val = parse_value(varargin{:});
            num=num+1;
            eval( sprintf( 'object.%s  = val;', valid_field(str) ) );
            if next_char == '}' || (count>=0 && num>=count)
                break;
            end
            %parse_char(',');
        end
    end
    if(count==-1)
        parse_char('}');
    end

%%-------------------------------------------------------------------------
function [cid,len]=elem_info(type)
id=strfind('iUIlLdD',type);
dataclass={'int8','uint8','int16','int32','int64','single','double'};
bytelen=[1,1,2,4,8,4,8];
if(id>0)
    cid=dataclass{id};
    len=bytelen(id);
else
    error_pos('unsupported type at position %d');
end
%%-------------------------------------------------------------------------

function [data adv]=parse_block(type,count,varargin)
global pos inStr isoct
[cid,len]=elem_info(type);
if(isoct)
    data=typecast(int8(inStr(pos:pos+len*count-1)),cid);
else
    data=typecast(uint8(inStr(pos:pos+len*count-1)),cid);
end
adv=double(len*count);

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

function object = parse_array(varargin) % JSON array is written in row-major order
global pos inStr isoct
    parse_char('[');
    object = cell(0, 1);
    dim=[];
    type='';
    count=-1;
    if(next_char == '$')
        type=inStr(pos+1);
        pos=pos+2;
    end
    if(next_char == '#')
        pos=pos+1;
        if(next_char=='[')
            dim=parse_array(varargin{:});
            count=prod(double(dim));
        else
            count=double(parse_number());
        end
    end
    if(~isempty(type))
        if(count>=0)
            [object adv]=parse_block(type,count,varargin{:});
            if(~isempty(dim))
                object=reshape(object,dim);
            end
            pos=pos+adv;
            return;
        else
            endpos=matching_bracket(inStr,pos);
            [cid,len]=elem_info(type);
            count=(endpos-pos)/len;
            [object adv]=parse_block(type,count,varargin{:});
            pos=pos+adv;
            parse_char(']');
            return;
        end
    end
    if next_char ~= ']'
         while 1
            val = parse_value(varargin{:});
            object{end+1} = val;
            if next_char == ']'
                break;
            end
            %parse_char(',');
         end
    end
    if(jsonopt('SimplifyCell',0,varargin{:})==1)
      try
        oldobj=object;
        object=cell2mat(object')';
        if(iscell(oldobj) && isstruct(object) && numel(object)>1 && jsonopt('SimplifyCellArray',1,varargin{:})==0)
            object=oldobj;
        elseif(size(object,1)>1 && ndims(object)==2)
            object=object';
        end
      catch
      end
    end
    if(count==-1)
        parse_char(']');
    end

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

function parse_char(c)
    global pos inStr len
    skip_whitespace;
    if pos > len || inStr(pos) ~= c
        error_pos(sprintf('Expected %c at position %%d', c));
    else
        pos = pos + 1;
        skip_whitespace;
    end

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

function c = next_char
    global pos inStr len
    skip_whitespace;
    if pos > len
        c = [];
    else
        c = inStr(pos);
    end

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

function skip_whitespace
    global pos inStr len
    while pos <= len && isspace(inStr(pos))
        pos = pos + 1;
    end

%%-------------------------------------------------------------------------
function str = parseStr(varargin)
    global pos inStr esc index_esc len_esc
 % len, ns = length(inStr), keyboard
    type=inStr(pos);
    if type ~= 'S' && type ~= 'C' && type ~= 'H'
        error_pos('String starting with S expected at position %d');
    else
        pos = pos + 1;
    end
    if(type == 'C')
        str=inStr(pos);
        pos=pos+1;
        return;
    end
    bytelen=double(parse_number());
    if(length(inStr)>=pos+bytelen-1)
        str=inStr(pos:pos+bytelen-1);
        pos=pos+bytelen;
    else
        error_pos('End of file while expecting end of inStr');
    end

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

function num = parse_number(varargin)
    global pos inStr len isoct
    id=strfind('iUIlLdD',inStr(pos));
    if(isempty(id))
        error_pos('expecting a number at position %d');
    end
    type={'int8','uint8','int16','int32','int64','single','double'};
    bytelen=[1,1,2,4,8,4,8];
    if(isoct)
        num=typecast(int8(inStr(pos+1:pos+bytelen(id))),type{id});
    else
        num=typecast(uint8(inStr(pos+1:pos+bytelen(id))),type{id});
    end
    pos = pos + bytelen(id)+1;

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

function val = parse_value(varargin)
    global pos inStr len
    true = 1; false = 0;

    switch(inStr(pos))
        case {'S','C','H'}
            val = parseStr(varargin{:});
            return;
        case '['
            val = parse_array(varargin{:});
            return;
        case '{'
            val = parse_object(varargin{:});
            if isstruct(val)
                if(~isempty(strmatch('x0x5F_ArrayType_',fieldnames(val), 'exact')))
                    val=jstruct2array(val);
                end
            elseif isempty(val)
                val = struct;
            end
            return;
        case {'i','U','I','l','L','d','D'}
            val = parse_number(varargin{:});
            return;
        case 'T'
            val = true;
            pos = pos + 1;
            return;
        case 'F'
            val = false;
            pos = pos + 1;
            return;
        case {'Z','N'}
            val = [];
            pos = pos + 1;
            return;
    end
    error_pos('Value expected at position %d');
%%-------------------------------------------------------------------------

function error_pos(msg)
    global pos inStr len
    poShow = max(min([pos-15 pos-1 pos pos+20],len),1);
    if poShow(3) == poShow(2)
        poShow(3:4) = poShow(2)+[0 -1];  % display nothing after
    end
    msg = [sprintf(msg, pos) ': ' ...
    inStr(poShow(1):poShow(2)) '<error>' inStr(poShow(3):poShow(4)) ];
    error( ['JSONparser:invalidFormat: ' msg] );

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

function str = valid_field(str)
global isoct
% From MATLAB doc: field names must begin with a letter, which may be
% followed by any combination of letters, digits, and underscores.
% Invalid characters will be converted to underscores, and the prefix
% "x0x[Hex code]_" will be added if the first character is not a letter.
    pos=regexp(str,'^[^A-Za-z]','once');
    if(~isempty(pos))
        if(~isoct)
            str=regexprep(str,'^([^A-Za-z])','x0x${sprintf(''%X'',unicode2native($1))}_','once');
        else
            str=sprintf('x0x%X_%s',char(str(1)),str(2:end));
        end
    end
    if(isempty(regexp(str,'[^0-9A-Za-z_]', 'once' ))) return;  end
    if(~isoct)
        str=regexprep(str,'([^0-9A-Za-z_])','_0x${sprintf(''%X'',unicode2native($1))}_');
    else
        pos=regexp(str,'[^0-9A-Za-z_]');
        if(isempty(pos)) return; end
        str0=str;
        pos0=[0 pos(:)' length(str)];
        str='';
        for i=1:length(pos)
            str=[str str0(pos0(i)+1:pos(i)-1) sprintf('_0x%X_',str0(pos(i)))];
        end
        if(pos(end)~=length(str))
            str=[str str0(pos0(end-1)+1:pos0(end))];
        end
    end
    %str(~isletter(str) & ~('0' <= str & str <= '9')) = '_';

%%-------------------------------------------------------------------------
function endpos = matching_quote(str,pos)
len=length(str);
while(pos<len)
    if(str(pos)=='"')
        if(~(pos>1 && str(pos-1)=='\'))
            endpos=pos;
            return;
        end
    end
    pos=pos+1;
end
error('unmatched quotation mark');
%%-------------------------------------------------------------------------
function [endpos e1l e1r maxlevel] = matching_bracket(str,pos)
global arraytoken
level=1;
maxlevel=level;
endpos=0;
bpos=arraytoken(arraytoken>=pos);
tokens=str(bpos);
len=length(tokens);
pos=1;
e1l=[];
e1r=[];
while(pos<=len)
    c=tokens(pos);
    if(c==']')
        level=level-1;
        if(isempty(e1r)) e1r=bpos(pos); end
        if(level==0)
            endpos=bpos(pos);
            return
        end
    end
    if(c=='[')
        if(isempty(e1l)) e1l=bpos(pos); end
        level=level+1;
        maxlevel=max(maxlevel,level);
    end
    if(c=='"')
        pos=matching_quote(tokens,pos+1);
    end
    pos=pos+1;
end
if(endpos==0)
    error('unmatched "]"');
end

author: plotly


text:
function feasible = checkfeasiblefast(p,x,tol)

feasible = 0;
if isempty(x)
    return
end
if ~(all(x < p.ub+tol) && all(x > p.lb-tol))
    return
end

if ~isempty(p.F_struc)
    vecres = p.F_struc*[1;x];

    if any(p.K.f)
        if any(-abs(vecres(1:p.K.f)) < - tol)
            return
        end
    end

    if any(p.K.l)
        if any(vecres(p.K.f+1:p.K.f+p.K.l) <-tol)
            return
        end
    end

    if any(p.K.q)
        top = startofSOCPCone(p.K);
        for i = 1:length(p.K.q)
            n = p.K.q(i);
            X = vecres(top:top+n-1);top = top+n;
            if X(1)-norm(full(X(2:end))) < -tol
                return
            end
        end
    end

    if any(p.K.e)
        top = startofEXPCone(p.K);
        for i = 1:p.K.e           ;
            X = vecres(top:top+2);top = top+3;
            if X(2)==0 && X(3) < -tol
                return
            elseif X(3) - X(2)*exp(X(1)/X(2)) < -tol
                return
            end
        end
    end

    if any(p.K.p)
        top = startofPOWCone(p.K);
        for i = 1:length(p.K.p)
            X = vecres(top:top+p.K.p(i)-1);
            a = X(end);
            top = top+p.K.p(i);
            if X(1)^a*X(2)^(1-a) - norm(X(3:end-1)) < -tol
                return
            end
        end
    end

    if any(p.K.s)
        top = startofSDPCone(p.K);
        for i = 1:length(p.K.s)
            n = p.K.s(i);
            X = reshape(vecres(top:top+n^2-1),n,n);top = top+n^2;
            X = (X+X')/2 + tol*eye(n);
            [~,fail] = chol(X);
            if fail
                return
            end
        end
    end
end
feasible = 1;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [varargout] = styleguide(varargin)
% STYLEGUIDE Open or search the MATLAB Programming Style Guide Wiki
%   STYLEGUIDE by itself opens the MATLAB Programming Style Guide Wiki in
%   the MATLAB web browser.
%
%   STYLEGUIDE SEARCHTERM1 SEARCHTERM2 SEARCHTERM3 ... searches the style
%   guide for any number of provided search terms and opens the search
%   results page in the MATLAB web browser.
%
%   STYLEGUIDE uses similar syntax to the WEB function with the exception
%   that it does not expect a URL. Rather, any number of search terms may
%   be used. All other inputs and outputs for WEB work for STYLEGUIDE, e.g.
%   [STAT, BROWSER, URL] = STYLEGUIDE(SEARCHTERM,'-NOTOOLBAR','-NEW').
%
%   To request permission to edit the wiki, use the contact form at
%   http://www.mathworks.com/matlabcentral/fileexchange/authors/contact/101715
%
%   Examples:
%      styleguide naming conventions
%         searches the style guide for "naming+conventions."
%
%      styleguide -browser
%         opens the style guide in the system web browser
%
%      [stat,h] = styleguide;
%         opens the style guide. Use close(h) to close the browser window.
%
%   See also WEB, HELP.

%   Author: Sky Sartorius
%           www.mathworks.com/matlabcentral/fileexchange/authors/101715

validWebArgs = {'-browser' '-display' '-helpbrowser' '-new' ...
    '-notoolbar' '-noaddressbox' '1'};

searchTerms = '';
isaWebArg = false(1,nargin);
for i = 1:nargin
    arg = strtrim(varargin{i});
    isWebArg = strcmpi({arg},validWebArgs);
    if any(isWebArg)
        isaWebArg(i) = true;
    else
        searchTerms = [searchTerms '+' arg]; %#ok<AGROW>
    end
end
if ~isempty(searchTerms)
    searchTerms(1) = '';
end

link = 'https://sites.google.com/site/matlabstyleguidelines/';

if ~isempty(searchTerms)
    searchStr = 'system/app/pages/search?scope=search-site&q=';
    link = [link searchStr searchTerms];
end

jnk = lower(varargin(isaWebArg));
[varargout{1:nargout}] = web(link,jnk{:});

% Revision history
%{
2013-03-13 first version
%}