text:
 function xs = wls_grpr(x, G, W, yi, D, xmin, xmax, niter)
%function xs = wls_grpr(x, G, W, yi, D, xmin, xmax, niter)
%
%	weighted least squares with constraint xmin <= x <= xmax
%	using gradient projection method (polyak:87 p. 207)
%		xnew = max(xmin, xold + D * \nabla J(xold))
%	cost function: J(x) = (y-Gx) W (y-Gx) / 2
%	in
%		x	[np,1]		initial estimate
%		G	[nd,np]		system matrix
%		W	[nd,nd]		weighting matrix
%		yi	[nd,1]		noisy measurements (e.g. sinogram)
%		D	[np,np]		preconditioning / step size matrix
%		xmin			minimum allowable x value
%		xmax			maximum allowable x value
%		niter			# of iterations
%	out
%		xs	[np,niter+1]	iterates
%
%	Copyright Dec. 2000,	Jeff Fessler, University of Michigan

if nargin < 2, ir_usage, end

if ~isvar('W') || isempty(W)
	W = 1;
end
if ~isvar('niter') || isempty(niter)
	niter = 1;
end
if ~isvar('D') || isempty(D)
	D = 1;
end
if ~isvar('xmin') || isempty(xmin), xmin = 0; end
if ~isvar('xmax') || isempty(xmax), xmax = inf; end

% loop over iterations
xs = zeros(length(x), niter);
xs(:,1) = x;
for iter = 2:niter
%	printf('WLS-PG iteration %d', iter-1)

	if ~rem(iter,2)
		lgrad = G' * (W .* (yi - G * x));

		x = x + D * lgrad;	% the update!
	else

		x = max(x,xmin);	% lower bound
		x = min(x,xmax);	% upper bound
	end

	xs(:,iter) = x;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function sVF = rotate(sVF, rot)
% rotate a function by a rotation
%
% Syntax
%   sVF = sVF.rotate(rot)
%
% Input
%  sVF - @S2VectorFieldHarmonic
%  rot - @rotation
%
% Output
%   sVF - @S2VectorFieldHarmonic
%

sVF.sF = rotate(sVF.sF, rot);

end

author: mtex-toolbox

repo: mtex

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


text:
function [Y1,U,err] = reduce_dim(Y)
%REDUCE_DIM Summary of this function goes here
%   Detailed explanation goes here
errors = [];
for ndims = 1:20
    [Y1,U,err] = pca_dr(Y,ndims);
    errors = [errors;err];
    if ndims > 3 && errors(end-1)-errors(end) < 1e-3*(errors(1)-errors(2))
        break;
    end
end

function [Y1,U,err] = pca_dr(Y,ndims)
[mappedY, mapping] = pca(Y, ndims);
U = mapping.M;
Y1 = Y*U;
recon_Y = Y1*U';
err = mean(mean(abs(recon_Y - Y)));

author: zhouyuanzxcv

repo: Hyperspectral


text:
% TIMEFDETAILS - details of the TIMEF function for time/frequency analysis
%                  of multiple epochs of single-channel event-related data.
%
% Global Description:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TIMEF performs normalized time/frequency averaging using either
% FFT-, wavelet-, or multitaper DFT estimates. The wavelets are N-cycle
% Hanning-windowed sinusoids. (Note: To substitute for HANNING windowing
% GAUSS or other windowing, replace the timef.m reference to HANNING).
%
% By default, the two image panels of the output plot show, respectively, the
% event-related spectral perturbation (ERSP) and inter-trial coherence (ITC)
% of the input data.
%
% The ERSP:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The ERSP (S. Makeig, Electroencephalogr Clin Neurophysiol 86:283-93, 1993)
% shows mean log event-loced deviations from epoch-mean (or baseline-mean) power
% at each frequency. If bootstrap statistics are computed (via the 'alpha'
% probability input parameter), (time, freq) points with non-significant
% differences from 0 are colored green in the image (but not in the 'ersp'
% output variable - use the output 'erspboot' variable to re-mask the 'ersp'
% output if desired). The baseline mean spectrum removed from each epoch
% is available as output parameter "powbase". Note that log(power) differences
% are equivalent to log(power ratios) between baseline and observed spectra,
% meaning the implicit ERSP model is one of (multiplicative) amplitude modulation
% of the EEG/MEG spectrum by e.g. subcortical and/or intra-cortical influences.
%
% In the default view, the thin bottom panel below the upper (ERSP) image shows
% the ERSP envelope (the most positive and most negative values at each output
% time point). The thin left panel shows the mean (or baseline) log spectrum
% (blue trace). When bootstrap statistics are computed (via the "alpha" argument),
% the left panel (green trace) also shows the bootstrap significant levels (+/-)
% at each frequency.
%
% The ITC:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (Inter-trial Coherence, cf. Tallon-Baudry et al., "Phase-locking factor")
% The lower panel shows the degree of tendency for spectral phase at each
% (time, freq) point to repeat across the data epochs. If bootstrap statistics
% are computed (as per the 'alpha' input parameter), non-significant points
% are colored green (but again not in the 'itc' output variable - use the
% itcboot output to re-mask if desired in later plotting).
%
% The lower thin panel shows the time-domain average (ERP) of the input data
% (blue) plus a zero-line (green). The average (ERP) is created principally by
% partial phase resetting of the EEG as measured by the ITC. (While phase resetting
% dominates, event-related spectral power changes (as measured by the ERSP)
% may also play a minor role). The thin left panel shows the frequency-mean ITC
% (blue trace) and, if bootstrap statistics are computed, the ITC significance
% limits at each frequency (green trace).
%
% ITC Math Derivation:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% By definition, linear coherence is
%       R=mean(Fxy)/sqrt(mean(abs(Fxx))*mean(abs(Fyy)));
% where Fxy is the cross-spectrum (FxFy*) and Fxx and Fyy the autospectra of
% processes x and y.  We define the phase coherence to be
%       R=mean(Fxy/(abs(Fxx))*abs(Fyy));  % mean of individually normed Fxy
% To derive the ITC, we consider y to be a stimulus-locked process
% such that, at each time and frequency, angle(Fyy) = 0 and abs(Fyy) = 1.
% Thus Pxy = Pxx, and the (complex) inter-trial phase coherence between x and
% the constant stimulus-locked process y is
%       ITC=mean(Pxx/abs(Pxx));
%
% USAGE:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   >> [ersp,itc,powbase,times,freqs,erspboot,itcboot]  ...
%                = timef(data,frames,tlimits,srate,cycles,...
%                              'key1',value1,'key2',value2, ... );
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOTE:
% * Left-click on subplots to view and zoom in separate windows (uses AXCOPY).
%
% Required Inputs:
%   data        = Single-channel data vector (1,frames*ntrials) (required)
%
%   frames      = Frames per trial                        {750}
%  Here, a frame is a data point (one channel at one time point) and the data
%  is assumed to be composed of concatenated epochs of the same length (frames).
%
%   tlimits     = [mintime maxtime] (ms) Epoch time limits {[-1000 2000]}
%  These should be the starting and ending times of the input epochs.
%
%   srate       = data sampling rate (Hz)                 {250}
%  This is sample rate per channel (i.e., data frames per second).
%
%   cycles      = >0 -> Number of cycles in each analysis wavelet
%                 =0 -> Use FFTs (with constant window length) {0}
%                 If [wavecycles factor] -> wavelet cycles increase with frequency
%                 beginning at wavecyles (0<factor<1; factor=1 -> no increase,
%                 standard wavelets; factor=0 -> fixed epoch length, as in FFT.
%                 OR multitaper decomposition (with 'mtaper').
%  Here, the user chooses either to use the FFT method (fixed window size
%  for all frequencies) or the wavelet DFT method (the data window length
%  depends inversely on the frequency, e.g. '3' means that the data
%  windows will each be three cycles wide (at each frequency). A higher
%  number here (e.g., '5') will narrow the frequency band and widen the
%  time window.
%
% Optional Inter-Irial Coherence Type:
%   'type'      = ['coher'|'phasecoher'] Compute either linear coherence
%                  ('coher') or phase coherence ('phasecoher') also known
%                  as the phase coupling factor           {'phasecoher'}.
%
% Optional Detrending:
%   'detret'    = ['on'|'off'], Detrend data in time.       {'off'}
%   'detrep'    = ['on'|'off'], Detrend data across trials (at each time point
%                 compute the linear trend across the set of ordered trials and
%                 remove it; not that this also subtract the ERP) {'off'}
%
% Optional FFT/DFT Parameters:
%   'winsize'   = If cycles==0: data subwindow length (fastest is 2^n < frames);
%                 If cycles >0: *longest* window length to use. This determines
%                 the lowest output frequency  {default: ~frames/8}
%  When cycles>0, winsize determines the lowest computed frequency. For example,
%  with srate=100 and cycles=3, a winsize of 100 means that 3 cycles must fit
%  within a 1-sec ( 100-sample) window. So, the lowest output frequency is 3 Hz.
%
%  When cycles=0, winsize is the length of data in each FFT window. This may be
%  extended with zeroes (to give more output frequencies) using padratio (below).
%
%  'timesout'  = Number of output times (int<frames-winframes) {200}
%  The number of FFTs or wavelet DFTs computed and plotted.
%
%   'padratio'  = FFT-length/winframes (2^k)                    {2}
%                  Multiplies the number of output frequencies by
%                  dividing their spacing. When cycles==0, frequency
%                  spacing is (low_freq/padratio).
%  This factor multiplies the number of output frequencies. In the FFT method
%  (cycles=0), this is done by zero-padding each analysis window. In the wavelet
%  DFT method (cycles>0), this gives the number of frequencies per Hz.
%
%   'maxfreq'   = Maximum frequency (Hz) to plot (& to output, if cycles>0)
%                  If cycles==0, all FFT frequencies are output. {50}
%   'baseline'  = Spectral baseline end-time (in ms). Use NaN for no baseline
%                 removal{0}
%   'powbase'   = Baseline spectrum to log-subtract. 'baseline' parameter is
%                 ignored if this parameter is used {def|NaN->from data}
%  This is useful only when you want to use a known baseline spectrum (e.g. from
%  another condition) instead of using the actual mean baseline spectrum of the data.
%  Otherwise, leave this out or specify as 'NaN' (not a number).
%
% Optional Multitaper Parameters:
%   'mtaper'    = If [N W], performs multitaper decomposition.
%                  (N is the time resolution and W the frequency resolution;
%                  maximum taper number is 2NW-1). Overwrites 'winsize' and
%                  'padratio'.
%                  If [N W K], forces the use of K Slepian tapers (if possible).
%                  Phase is calculated using standard methods.
%                  The use of mutitaper with wavelets (cycles>0) is not
%                  recommended (as multiwavelets are not implemented).
%                  Uses Matlab functions DPSS, PMTM.   {no multitaper}
%
% Optional Bootstrap Parameters:
%   'alpha'     = If non-0, compute two-tailed bootstrap significance prob.
%                  level. Show non-signif. output values as green.   {0}
%  This optional parameter lengthens the computation time but gives bootstrap
%  estimates of which ERSP and ITC values are significantly different from 0,
%  by setting to 0 (green) all non-significant values in the ERSP and ITC images.
%  Normal values for alpha are 0 ([], or none) -> no bootstrap computation, or
%  0.01 (which should allow about 1% of random images to appear "significant"
%
%   'naccu'     = Number of bootstrap replications to accumulate     {200}
%   'baseboot'  = Bootstrap baseline subtract (0 -> use 'baseline';
%                                                  1 -> use whole trial) {0}
% Optional Scalp Map Plotting Parameters:
%   'topovec'   = Scalp topography (map) to plot                     {none}
%   'elocs'     = Electrode location file for scalp map   {no default}
%                     File should be ascii in format of  >> topoplot example
%  This is an optional map-plotting feature. Given an input map vector
%  (one weight at each channel, and a electrode location file, TIMEF plots
%  a TOPOPLOT-style 2-d scalp map on the left side of the figure. See
%  >> topoplot example % for the format of the electrode location file.
%
% Other Optional Plotting Parameters:
%   'vert'      = [vector of ms times] -> plot vertical dashed lines  {0 only}
%  Use this to add extra vertical dashed lines at significant epoch times.
%  Time 0 is marked by default.
%
%   'plotersp'  = ['on'|'off'] Plot power spectral perturbations    {'on'}
%   'plotitc'   = ['on'|'off'] Plot inter trial coherence            {'on'}
%   'title'     = Optional figure title                              {none}
%
%   'pboot'     = Bootstrap power limits (e.g., from TIMEF)   {from data}
%   'rboot'     = Bootstrap ITC limits (e.g., from TIMEF)     {from data}
%  These are useful if you want to apply significance limits from another condition
%  to new data. {default|NaN, compute from data}
%
%   'linewidth' = Line width for 'marktimes' traces (thick=2, thin=1) {2}
%   'axesfont'  = Axes text font size                                {10}
%   'titlefont' = Title text font size                               {8}
%
% Outputs:
%        ersp   = Matrix (nfreqs,timesout) of log spectral diffs. from baseline (dB)
%        itc    = Matrix (nfreqs,timesout) of inter-trial phase coherence (range: [0 1])
%   Note that when cycles=0, nfreqs is total number of FFT frequencies, which
%   typically include frequencies higher than maxfreq. When cycles>0, *no* extra
%   (higher) frequencies are computed.
%
%      powbase  = Baseline power spectrum (removed to compute the ERSP)
%        times  = Vector of output times (subwindow centers) (in ms).
%        freqs  = Vector of frequency bin centers (in Hz).
%
%     erspboot  = Matrix (2,nfreqs) of [lower;upper] ERSP significance diffs.
%      itcboot  = Matrix (2,nfreqs) of [lower;upper] ITC thresholds (abs., not diffs)
%  Note that the itcboot lower bound is practically meaningless.
%
%  Plot description:
%    Assuming both 'plotersp' and 'plotitc' options are 'on' (= default). The upper panel
%    presents the data ERSP (Event-Related Spectral Perturbation) in dB, with mean baseline
%    spectral activity (in dB) subtracted. Use "'baseline', NaN" to prevent TIMEF from
%    removing the baseline. The lower panel presents the data ITC (Inter-Trial Coherence).
%    Click on any plot axes to pop up a new window (using 'AXCOPY')
%    -- Upper left marginal panel presents the mean spectrum during the baseline period
%       (blue), and when significance is set, the significance threshold at each frequency
%       (dotted green-black trace).
%    -- The marginal panel under the ERSP image shows the maximum (green) and minimum
%       (blue) ERSP values relative to baseline power at each frequency.
%    -- The lower left marginal panel shows mean ITC across the imaged time range (blue),
%       and when significance is set, the significance threshold (dotted green-black).
%    -- The marginal panel under the ITC image shows the ERP (which is produced by ITC
%       across the data spectral pass band).
%
% Author: Sigurd Enghoff, Arnaud Delorme & Scott Makeig
%          CNL / Salk Institute 1998- | SCCN/INC, UCSD 2002-
%
% See also: CROSSF - event-related cross-spectral coherence between two input
%                      time series.
%
% History:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TIMEF was coded by Sigurd Enghoff and Scott Makeig at The Salk
% Institute, La Jolla CA in August, 1998, using methods developed
% in Makeig, 1993. Arno Delorme added the multitaper option, recoded
% the function to use 'keyword','parameter' argument pairs, and added the
% 'type' argument with advice from Joern Anemueller at SCCN/Institute for
% Neural Computation, UCSD in early 2002.

% Copyright (C) 8/01/00 Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 01-25-02 reformated help & license -ad

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function y = conj( x )

%Disciplined convex/geometric programming information for CONJ:
%   CONJ(X) imposes no convexity restrictions on its arguments. However,
%   since CONJ(X)=X when X is real, it is only useful for complex
%   affine expressions.

y = cvx( x.size_, conj( x.basis_ ) );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function varargout = rodr(varargin)
% VL_RODR  Rodrigues' formula
%   R = VL_RODR(OM) where OM a 3-dimensional column vector computes the
%   Rodrigues' formula of OM, returning the rotation matrix R =
%   expm(vl_hat(OM)).
%
%   [R,DR] = VL_RODR(OM) computes also the derivative of the Rodrigues
%   formula. In matrix notation this is the expression
%
%           d(vec expm(vl_hat(OM)) )
%     dR = ----------------------.
%                  d om^T
%
%   [R,DR]=VL_RODR(OM) when OM is a 3xK matrix repeats the operation for
%   each column (or equivalently matrix with 3*K elements). In this
%   case R and DR are arrays with K slices, one per rotation.
%
%   See also: VL_IRODR(), VL_HELP().
[varargout{1:nargout}] = vl_rodr(varargin{:});

author: yihui-he

repo: panorama

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


text:
function p=lin2pcmu(x,s)
%LIN2PCMU Convert linear to Mu-law PCM P=(X,S)
%	pcmu = lin2pcmu(lin) where lin contains a vector
%	or matrix of signal values within a range determined by
%	the scale factor s (see table below).
%	Values outside this range will be clipped.
%	The input values will be converted to integer
%	Mu-law pcm vlues in the range 0 to 255.
%
%	Input values are multiplied by the scale factor s:
%
%		   s		Input Range
%
%		   1		+-8159
%		4004.189931	+-2.03761563 (default)
%		8159		+-1
%
%	The default input scaling factor 4004.189931 is equal to
%	sqrt((2207^2 + 5215^2)/2) and follows ITU standard G.711.
%	The sine wave with PCM-Mu values [158 139 139 158 30 11 11 30]
%	has a mean square value of unity corresponding to 0 dBm0.

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: lin2pcmu.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 2 s=4004.189931; end
y=x*s;
y=(abs(y+8031)-abs(y-8031))/2;
q=floor((y+8032)/8032);
[m,e]=log2(abs(y)+33);
p=175+128*q-8*(e+abs(e-6))-floor(32*m-16);

author: decouples


text:
function y = exp( x )

%   Disciplined convex programming information:
%       EXP(X) is convex and nondecreasing in X. When used in CVX
%       expressions, X must be real. Typically, X must also be affine
%       or convex; X can also be concave, but this produces a log-concave
%       result with very limited usefulness.

global cvx___
error(nargchk(1,1,nargin));
cvx_expert_check( 'log', x );

%
% Determine the expression types
%

persistent remap
if isempty( remap ),
    remap_1 = cvx_remap( 'real' );
    remap_2 = cvx_remap( 'convex', 'concave' ) & ~remap_1;
    remap = remap_1 + 2 * remap_2;
end
v = remap( cvx_classify( x ) );

%
% Process each type of expression one piece at a time
%

vu = sort( v(:) );
vu = vu([true;diff(vu)~=0]);
nv = length( vu );
sx = x.size_;
if nv ~= 1,
    y = cvx( sx, [] );
end
for k = 1 : nv,

    %
    % Select the category of expression to compute
    %

    vk = vu( k );
    if nv == 1,
        xt = x;
    else
        t = v == vk;
        xt = cvx_subsref( x, t );
    end

    %
    % Perform the computations
    %

    switch vk,
        case 0,
            % Invalid
            error( 'Disciplined convex programming error:\n    Illegal operation: exp( {%s} ).', cvx_class( xt ) );
        case 1,
            % Constant
            xt = cvx( exp( cvx_constant( xt ) ) );
        case 2,
            % Affine, convex, concave
            xt = sparsify( xt, 'exponential' );
            [ rx, cx, vx ] = find( xt.basis_ );
            tt = rx == 1;  rx( tt ) = [];
            cc = cx( tt ); cx( tt ) = [];
            vc = vx( tt ); % vx( tt ) = [];
            exps = cvx___.exponential( rx, 1 );
            tt = exps == 0;
            if any( tt ),
                n1 = unique( rx( tt ) );
                n2 = newvar( cvx___.problems( end ).self, '', length( n1 ) );
                [ n2, dummy ] = find( n2.basis_ ); %#ok
                cvx___.exponential( n1, 1 ) = n2( : );
                cvx___.logarithm( n2, 1 ) = n1( : );
                cvx___.vexity( n2 ) = 1;
                n2 = n2( cvx___.vexity( n1 ) < 0 );
                cvx___.vexity( n2 ) = NaN;
                cvx___.nan_used = true;
                cvx___.canslack( n2 ) = false;
                exps = cvx___.exponential( rx, 1 );
            end
            nb = size( xt.basis_, 2 );
            bx = sparse( exps, cx, 1, max( exps ), nb );
            if ~isempty( cc ),
                bx = bx * diag(exp(sparse(cc,1,vc,nb,1)));
            end
            xt = cvx( xt.size_, bx );
        otherwise,
            error( 'Shouldn''t be here.' );
    end

    %
    % Store the results
    %

    if nv == 1,
        y = xt;
    else
        y = cvx_subsasgn( y, t, xt );
    end

end

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

repo: REST


text:
function taun = tau ( n, taun )

%*****************************************************************************80
%
%% TAU returns the value of TAU(N), the number of distinct divisors of N.
%
%  Discussion:
%
%    TAU(N) is the number of divisors of N, including 1 and N.
%
%  First values:
%
%     N   TAU(N)
%
%     1    1
%     2    2
%     3    2
%     4    3
%     5    2
%     6    4
%     7    2
%     8    4
%     9    3
%    10    4
%    11    2
%    12    6
%    13    2
%    14    4
%    15    4
%    16    5
%    17    2
%    18    6
%    19    2
%    20    6
%
%  Formula:
%
%    If the prime factorization of N is
%
%      N = P1^E1 * P2^E2 * ... * PM^EM,
%
%    then
%
%      TAU(N) = ( E1 + 1 ) * ( E2 + 1 ) * ... * ( EM + 1 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the value to be analyzed.  N must be 1 or
%    greater.
%
%    Output, integer TAUN, the value of TAU(N).  But if N is 0 or
%    less, TAUN is returned as 0, a nonsense value.  If there is
%    not enough room for factoring, TAUN is returned as -1.
%
  if ( n <= 0 )
    taun = 0;
    return
  end

  if ( n == 1 )
    taun = 1;
    return
  end
%
%  Factor N.
%
  [ nfactor, factor, power, nleft ] = i4_factor ( n );

  if ( nleft ~= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TAU - Fatal error!\n' );
    fprintf ( 1, '  Not enough factorization space.\n' );
    taun = -1;
    error ( 'TAU - Fatal error!' );
  end

  taun = 1;
  for i = 1 : nfactor
    taun = taun * ( power(i) + 1 );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clear
clc

A0=imread('cam.jpg');
% chk=[1 1 1 3 1 1 2 1 1 2 3 2 2 2 1 1 1 2 3 1 1 4 1 3 2 1 1,...
%     1 1 1 1 1 2 1 2 2 3 1 1 2 1 2 3 2 1 4 1 1 1 1 3 2 2 2 2,...
%     1 1 1 1]; % for cam1
% chk=[1 1 1 3 1 1 2 1 1 2 3 2 1 2 2 1 1 2 3 3 1 2 1 3 2, ...
%     1 1 1 1 1 1 1 3 2 1 1 2 2 2 1 2 2 2 1 1 1 3 2 1 1 3 2, ...
%     1 2 3 1 1 1 1]; % for cam
% chk=[1 1 1 3 1 1 2 1 1 2 3 1 1 1 4 1 1 2 3 1 1 4 1 1, ...
%     1 1 4 1 1 1 1 1 1 1 1 4 1 3 1 2 3 2 1 1 1 1 3 2 2,...
%     2 2 1 3 2 1 1 1 1 1]; % for cam3
% chk=[1 1 1 1 3 1 2 3 1 2 1 1 1 2 3 3 2 1 1 2 1 3 1 3 2 1,...
%     1 1 1 1 1 1 2 2 2 1 1 1 3 2 1 1 3 2 1 3 1 2 3 1 1 2 1,...
%     2 1 3 1 1 1]; % for cambook3
A0=rgb2gray(A0);
% A=A0([90:330],[1:604]); % for cam1
A=A0([99+30:99+260],[71:608]); % for cam
% A=A0([182:389],[1:600]); % for cam3
% A=A0([193:349],[1:504]); % for cambook3
return
figure(1), imshow(A)%, return

nSeg=10;
wSeg=floor(size(A,1)/nSeg);

for i=1:nSeg

    At=A((i-1)*wSeg+1:i*wSeg,:);
    %figure(2), imshow(At)

    gsv(i,:)=255-mean(double(At),1);
    x=[1:length(gsv(i,:))]; g=gsv(i,:);
    [gw,xw]=movAvg(g,1);
    %figure(3), plot(x,gsv(i,:),xw,gw)

    dgw=diffScheme(gw);
    [pk1,loc1]=peakDetect(dgw,20,5);
    % sub-pixel edge detection
    for i1=1:length(loc1)
        [X1(i1),P1(i1)]=subPx(dgw(loc1(i1)-5:loc1(i1)+5),...
            xw(loc1(i1)-5:loc1(i1)+5),...
            dgw(loc1(i1)),100);
        %g=dgw(loc1(i)-5:loc1(i)+20); gMax=dgw(loc1(i)); save subPxData.mat g gMax
    end
    [pk2,loc2]=peakDetect(-dgw,20,5);
    for i1=1:length(loc2)
        [X2(i1),P2(i1)]=subPx(-dgw(loc2(i1)-5:loc2(i1)+5),...
            xw(loc2(i1)-5:loc2(i1)+5),...
            -dgw(loc2(i1)),100);
    end
    P2=-P2;

    if length(loc1)~=length(loc2)
        disp('oops! ... left edge & right edge don''t match!!');
    end
    %figure(4), plot(xw,dgw,'.k-',xw(loc1),dgw(loc1),'o',...
    %    xw(loc2),dgw(loc2),'o',X1,P1,'o',X2,P2,'o')

    % black gap
    bGap=X2-X1;
    % white gap
    wGap=-X2(1:end-1)+X1(2:end);

    % discretizing barcodes
    %figure(5), subplot(2,1,1), plot(sort(bGap),'ok')
    %subplot(2,1,2), plot(sort(wGap),'or'), pause(0.1)

    % working on white gaps
    [wGap2,wid]=sort(wGap);
    dwGap2=wGap2(2:end)-wGap2(1:end-1);
    [wGap21,dwid]=sort(dwGap2);
    bound(1:3)=sort(dwid(end:-1:end-2)+1);
    clear wGap21 dwid
    s1=mean(wGap2(1:bound(1)-1));
    s2=mean(wGap2(bound(1):bound(2)-1));
    s3=mean(wGap2(bound(2):bound(3)-1));
    s4=mean(wGap2(bound(3):end));
    if s4/s1 <3.5 % 4 is not present
        if s2/s1>1.5
            S1=wid(1:bound(1)-1); % s1 and s2 are exclusive
            if s3/s2<1.2 % s2 and s3 are same
                S2=wid(bound(1):bound(3)-1);
                S3=wid(bound(3):end);
            else % s2 and s3 are separate
                S2=wid(bound(1):bound(2)-1);
                S3=wid(bound(2):end);
            end
        else
            S1=wid(1:bound(2)-1); % s1 and s2 are inclusive
            S2=wid(bound(2)-1:bound(3));
            S3=wid(bound(3)-1:end);
        end
        S4=[];
    else % 4 is present
        S1=wid(1:bound(1)-1);
        S2=wid(bound(1):bound(2)-1);
        S3=wid(bound(2):bound(3)-1);
        S4=wid(bound(3):end);
    end
    %figure(6),plot(wGap(S1),'o'),hold on,plot(wGap(S2),'.'),plot(wGap(S3),'x'),hold off
    wGap3(S1)=1;
    wGap3(S2)=2;
    wGap3(S3)=3;
    if ~isempty('S4'), wGap3(S4)=4; end

    % working on black gaps
    [bGap2,bid]=sort(bGap);
    dbGap2=bGap2(2:end)-bGap2(1:end-1);
    [bGap21,dbid]=sort(dbGap2);
    bound(1:3)=sort(dbid(end:-1:end-2)+1);
    clear bGap21 dbid
    t1=mean(bGap2(1:bound(1)-1));
    t2=mean(bGap2(bound(1):bound(2)-1));
    t3=mean(bGap2(bound(2):bound(3)-1));
    t4=mean(bGap2(bound(3):end));
    if t4/t1 <3.5 % 4 is not present
        if t2/t1>1.5
            T1=bid(1:bound(1)-1); % t1 and t2 are exclusive
            if t3/t2<1.2 % t2 and t3 are same
                T2=bid(bound(1):bound(3)-1);
                T3=bid(bound(3):end);
            else % t2 and t3 are separate
                T2=bid(bound(1):bound(2)-1);
                T3=bid(bound(2):end);
            end
        else
            T1=bid(1:bound(2)-1); % t1 and t2 are inclusive
            T2=bid(bound(2)-1:bound(3));
            T3=bid(bound(3)-1:end);
        end
        T4=[];
    else % 4 is present
        T1=bid(1:bound(1)-1);
        T2=bid(bound(1):bound(2)-1);
        T3=bid(bound(2):bound(3)-1);
        T4=bid(bound(3):end);
    end
    %figure(7),plot(bGap(T1),'o'),hold on,plot(bGap(T2),'.'),plot(bGap(T3),'x'),hold off
    bGap3(T1)=1;
    bGap3(T2)=2;
    bGap3(T3)=3;
    if ~isempty('T4'), bGap3(T4)=4; end

    %wGap3=round(wGap/mean(wGap2(1:10)));
    %bGap3=round(bGap/mean(bGap2(1:10)));

    I(i).gap=zeros(length(bGap)+length(wGap),1);
    I(i).gap(1:2:end)=bGap3;
    I(i).gap(2:2:end)=wGap3;

    %disp(['segment: ',num2str(i),' | errors: ',num2str(sum((I(i).gap'-chk)~=0))])
    if i==1
        G=I(i).gap;
    else
        for j=1:length(I(i).gap)
            G(j)=G(j)+I(i).gap(j);
        end
    end

end

G=round(G/nSeg);
barcode=barcodeEAN13(G);
%reading_error=sum(sum(G'-chk))

author: Sable

repo: mcbench-benchmarks


text:
function create_table_of_algorithms(up)
%CREATE_TABLE_OF_ALGORITHMS creates a table of algorithms and their performances
%	            create_table_of_algorithms(up)
%
%	Inputs:
%		data            data files should be stored in the specified format
%       up              universal parameters structure
%
%	Outputs:
%       tables          csv tables of results
%

fprintf('\n--- Creating Results Tables ');

% save setup
save_name = up.paths.filenames.results_table;
savepath = [up.paths.tables_save_folder, save_name, '.mat'];
if ~up.analysis.redo_stats
    exist_log = check_exists(savepath, save_name);
    if exist_log
        return
    end
end

% Load BA Data
load_name = 'BA_results';
loadpath = [up.paths.data_save_folder, up.paths.filenames.global_BA, '.mat'];
load(loadpath, load_name);

% Load algorithm names
load_name = up.paths.filenames.alg_names;
loadpath = [up.paths.data_save_folder, up.paths.filenames.alg_names, '.mat'];
load(loadpath, load_name);

% Load study data
load_name = 'trad_stats';
loadpath = [up.paths.data_save_folder, up.paths.filenames.study_stats, '.mat'];
load(loadpath, load_name);

% identify groups
groups = fieldnames(BA_results);

%% Create a BA table and a Stats table for each group

for group_no = 1 : length(groups)
    group = groups{group_no};

    %% Extract relevant data
    rel_data = extract_rel_data(group, BA_results, trad_stats, up);

    %% Make Table of BA Results
    BA_results_table = create_table_results(alg_names, rel_data, up, group, 'BA');

    %% Save to file
    save_name = [up.paths.filenames.BA_results_table, '_', group];
    savepath = [up.paths.tables_save_folder, save_name];
    xlswrite(savepath, BA_results_table)

    %% Make Table of Stats Results
    stats_results_table = create_table_results(alg_names, rel_data, up, group, 'stats');

    %% Save to file
    save_name = [up.paths.filenames.stats_results_table, '_', group];
    savepath = [up.paths.tables_save_folder, save_name];
    xlswrite(savepath, stats_results_table)

end

end

function rel_data = extract_rel_data(group, BA_results, stats, up)

rel_data.BA = extractfield(BA_results, group); rel_data.BA = rel_data.BA{1,1};
rel_data.stats = extractfield(stats, group); rel_data.stats = rel_data.stats{1,1};

end

function table_data = create_table_results(alg_names, rel_data, up, group, type)

% Find alg numbers:
alg_names.no = cell(length(alg_names.names),1);
for s = 1 : length(alg_names.no)
    alg_names.no{s} = ['a', num2str(s)];
end

% Extract relevant variables for table:
tab_vars.alg_no = alg_names.no;
if isfield(alg_names.meths, 'xa')
    tab_vars.m_xa = alg_names.meths.xa;
    tab_vars.m_xb = alg_names.meths.xb;
    tab_vars.m_ef = alg_names.meths.ef;
    tab_vars.m_et = alg_names.meths.et;
    tab_vars.m_fm = alg_names.meths.fm;
    tab_vars.m_ft = alg_names.meths.ft;
    % modify the fme:
    rel_els = ~isnan(tab_vars.m_fm);
    tab_vars.m_xb(rel_els) = 99;
end
tab_vars.alg_name = alg_names.names;
tab_vars.signal = alg_names.sigs;

% Extract relevant statistics for table:
if strcmp(type, 'BA')
    tab_vars.bias = rel_data.BA.bias.val;
    tab_vars.two_sd = rel_data.BA.two_sd.val;
    tab_vars.prop_wins_est = rel_data.stats.prop_wins_est;
    tab_vars.prop_wins_ref_and_good_sqi = rel_data.stats.prop_wins_good_sqi_and_ref;
elseif strcmp(type, 'stats')
    tab_vars.percerr = rel_data.stats.percerr;
    tab_vars.mae = rel_data.stats.mae;
    tab_vars.sdae = rel_data.stats.sdae;
    tab_vars.rmse = rel_data.stats.rmse;
    tab_vars.cp2 = rel_data.stats.cp2;
    tab_vars.prop_wins_est = rel_data.stats.prop_wins_est;
    tab_vars.total_wins_all = rel_data.stats.total_wins_all;
end

%% Determine what goes into the table
headers = fieldnames(tab_vars);
table_data = cell(1+length(tab_vars.alg_name), length(headers));
table_data(1,:) = headers;
for col_no = 1 : length(headers)
    eval(['col_data = tab_vars.' headers{col_no} ';']);
    if ~iscell(col_data)
        table_data(2:end,col_no) = num2cell(col_data);
    else
        table_data(2:end,col_no) = col_data;
    end
end

end

function [subj_data, BA_data, stats_data, old_names] = load_data(group, up)

% Load subj data
load_name = ['data_', group];
loadpath = [up.paths.data_save_folder, up.paths.filenames.group_data, group, '.mat'];
load(loadpath, load_name)
eval(['subj_data = data_' group  ';']);

% Load B-A results
load_name = 'BA';
loadpath = [up.paths.data_save_folder, up.paths.filenames.group_BA, group, '.mat'];
load(loadpath, load_name)
BA_data = BA; clear BA

% Load Stats data
load_name = 'acc';
loadpath = [up.paths.data_save_folder, up.paths.filenames.group_stats, group, '.mat'];
load(loadpath, load_name)
stats_data = acc;
old_names = acc.rr_alg_names;

end

function create_table_stats_results_subj(alg_names, sig_names, old_names, subj_data, up, group)

%% Find prop with an abs error of < 1 bpm
eval(['subjs = unique(subj_data.' old_names{1} '.subj);']);
prop_acc = zeros(length(old_names), length(subjs));
for alg_no = 1 : length(old_names)
     eval(['rel_data = subj_data.' old_names{alg_no} ';']);

     rel_data.error = rel_data.est - rel_data.ref;
     rel_data.abs_error = abs(rel_data.error);
     rel_data.acc_log = rel_data.abs_error < 1;

     for subj_no = 1:length(subjs)
         subj = subjs(subj_no);
         rel_els = logical(rel_data.subj == subj);
         prop_acc(alg_no, subj_no) = sum(rel_data.acc_log(rel_els)) / length(rel_data.acc_log(rel_els));
     end

     % extract relevant parts of algs:
     for comp_str = {'FMe', 'Flt', 'RSa', 'PDt', 'Est', 'MFu', 'TFu'}
         if ~exist([lower(comp_str{1,1}), '_meth'], 'var')
             eval([lower(comp_str{1,1}), '_meth = nan(length(old_names),1);']);
         end
         rel_el = strfind(alg_names{alg_no}, comp_str{1,1});
         if rel_el
             temp = str2num(alg_names{alg_no}(rel_el+3));
             eval([lower(comp_str{1,1}), '_meth(alg_no) = temp;']);
         end
     end

     clear rel_data
end

%% Determine what goes into the table
headers = {'Signal', 'Algorithm Name', 'Old Name', 'FMe', 'Flt', 'RSa', 'PDt', 'Est', 'MFu', 'TFu'};
for subj_no = 1:length(subjs)
    headers{end+1} = ['Prop_acc ' num2str(subjs(subj_no))];
end

table_data = cell(1+length(alg_names), length(headers));
table_data(1,:) = headers;
table_data(2:end,1) = sig_names;
table_data(2:end,2) = alg_names;
table_data(2:end,3) = old_names;
table_data(2:end,4) = num2cell(fme_meth);
table_data(2:end,5) = num2cell(flt_meth);
table_data(2:end,6) = num2cell(rsa_meth);
table_data(2:end,7) = num2cell(pdt_meth);
table_data(2:end,8) = num2cell(est_meth);
table_data(2:end,9) = num2cell(mfu_meth);
table_data(2:end,10) = num2cell(tfu_meth);
for s = 1 : size(prop_acc,2)
    table_data(2:end,10+s) = num2cell(prop_acc(:,s));
end

%% Save to file
save_name = [up.paths.filenames.subj_results_table '_' group];
savepath = [up.paths.tables_save_folder, save_name];
xlswrite(savepath,table_data)

end

function create_table_BA_results(alg_names, sig_names, BA_data, up, group)

%% Determine what goes into the table
headers = {'Signal', 'Algorithm Name', 'Bias', 'Bias (LCI)', 'Bias (UCI)', 'LLOA', 'LLOA (LCI)', 'LLOA (UCI)', 'ULOA', 'ULOA (LCI)', 'ULOA (UCI)', 'Precision', 'Precision (LCI)', 'Precision (UCI)'};
table_data = cell(1+length(alg_names), length(headers));
table_data(1,:) = headers;
table_data(2:end,1) = sig_names;
table_data(2:end,2) = alg_names;
table_data(2:end,3) = num2cell(BA_data.bias.val);
table_data(2:end,4) = num2cell(BA_data.bias.lci);
table_data(2:end,5) = num2cell(BA_data.bias.uci);
table_data(2:end,6) = num2cell(BA_data.lloa.val);
table_data(2:end,7) = num2cell(BA_data.lloa.lci);
table_data(2:end,8) = num2cell(BA_data.lloa.uci);
table_data(2:end,9) = num2cell(BA_data.uloa.val);
table_data(2:end,10) = num2cell(BA_data.uloa.lci);
table_data(2:end,11) = num2cell(BA_data.uloa.uci);
table_data(2:end,12) = num2cell(BA_data.prec.val);
table_data(2:end,13) = num2cell(BA_data.prec.lci);
table_data(2:end,14) = num2cell(BA_data.prec.uci);

%% Save to file
save_name = [up.paths.filenames.BA_results_table '_' group];
savepath = [up.paths.tables_save_folder, save_name];
xlswrite(savepath,table_data)

end

function create_table_stats_results(alg_names, sig_names, stats_data, up, group)

%% Determine what goes into the table
headers = {'Signal', 'Algorithm Name', 'MAE', 'MAE (lq)', 'MAE (uq)', 'SDAE', 'SDAE (lq)', 'SDAE (uq)', 'PERC ERR', 'PERC ERR (lq)', 'PERC ERR (uq)', 'RMSE', 'RMSE (lq)', 'RMSE (uq)', 'MEAN ERROR', 'MEAN ERROR (lq)', 'MEAN ERROR (uq)', 'STD ERROR', 'STD ERROR (lq)', 'STD ERROR (uq)', 'PROP WINS', 'PROP WINS (lq)', 'PROP WINS (uq)'};
table_data = cell(1+length(alg_names), length(headers));
table_data(1,:) = headers;
table_data(2:end,1) = sig_names;
table_data(2:end,2) = alg_names;
table_data(2:end,3) = num2cell(stats_data.mae.med);
table_data(2:end,4) = num2cell(stats_data.mae.lq);
table_data(2:end,5) = num2cell(stats_data.mae.uq);
table_data(2:end,6) = num2cell(stats_data.sdae.med);
table_data(2:end,7) = num2cell(stats_data.sdae.lq);
table_data(2:end,8) = num2cell(stats_data.sdae.uq);
table_data(2:end,9) = num2cell(stats_data.percerr.med);
table_data(2:end,10) = num2cell(stats_data.percerr.lq);
table_data(2:end,11) = num2cell(stats_data.percerr.uq);
table_data(2:end,12) = num2cell(stats_data.rmse.med);
table_data(2:end,13) = num2cell(stats_data.rmse.lq);
table_data(2:end,14) = num2cell(stats_data.rmse.uq);
table_data(2:end,15) = num2cell(stats_data.mean_diff.med);
table_data(2:end,16) = num2cell(stats_data.mean_diff.lq);
table_data(2:end,17) = num2cell(stats_data.mean_diff.uq);
table_data(2:end,18) = num2cell(stats_data.std_diff.med);
table_data(2:end,19) = num2cell(stats_data.std_diff.lq);
table_data(2:end,20) = num2cell(stats_data.std_diff.uq);
table_data(2:end,21) = num2cell(stats_data.prop_wins.med);
table_data(2:end,22) = num2cell(stats_data.prop_wins.lq);
table_data(2:end,23) = num2cell(stats_data.prop_wins.uq);

%% Save to file
save_name = [up.paths.filenames.stats_results_table '_' group];
savepath = [up.paths.tables_save_folder, save_name];
xlswrite(savepath,table_data)

end

author: peterhcharlton

repo: RRest


text:
%PDF2EPS  Convert a pdf file to eps format using pdftops
%
% Examples:
%   pdf2eps source dest
%
% This function converts a pdf file to eps format.
%
% This function requires that you have pdftops, from the Xpdf suite of
% functions, installed on your system. This can be downloaded from:
% http://www.foolabs.com/xpdf
%
%IN:
%   source - filename of the source pdf file to convert. The filename is
%            assumed to already have the extension ".pdf".
%   dest - filename of the destination eps file. The filename is assumed to
%          already have the extension ".eps".

% Copyright (C) Oliver Woodford 2009-2010

% Thanks to Aldebaro Klautau for reporting a bug when saving to
% non-existant directories.

function pdf2eps(source, dest)
% Construct the options string for pdftops
options = ['-q -paper match -eps -level2 "' source '" "' dest '"'];
% Convert to eps using pdftops
[status message] = pdftops(options);
% Check for error
if status
    % Report error
    if isempty(message)
        error('Unable to generate eps. Check destination directory is writable.');
    else
        error(message);
    end
end
% Fix the DSC error created by pdftops
fid = fopen(dest, 'r+');
if fid == -1
    % Cannot open the file
    return
end
fgetl(fid); % Get the first line
str = fgetl(fid); % Get the second line
if strcmp(str(1:min(13, end)), '% Produced by')
    fseek(fid, -numel(str)-1, 'cof');
    fwrite(fid, '%'); % Turn ' ' into '%'
end
fclose(fid);
return

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function SO3F = sqrt(SO3F, varargin)
% square root of a SO3Fun
%
% Syntax
%   SO3F = sqrt(SO3F)
%
% Input
%  SO3F - @SO3FunHarmonic
%
% Output
%  SO3F - @SO3FunHarmonic
%

SO3F = SO3F.^(1/2);

end


text:
function [report] = TrackDetect(Track, UAV)
%TRACKDETECT 判断航迹是否符合条件(一个agent的)

%% 无人机航迹检测

% 威胁检测
dim = UAV.PointDim;                                         % 仿真环境维度
M = [UAV.Menace.radar; UAV.Menace.other]; % 威胁区

Threat = cell(UAV.num, 1);                                 % 威胁结果树
Angle = cell(UAV.num, 1);                                  % 转角检测结果树
MiniTraj = cell(UAV.num, 1);                              % 最小航迹片段检测结果树
ProbPoint = cell(UAV.num, 1);                           % 问题点

L = cell(UAV.num, 1);                                          % 航迹片段树（累加结构）
Time = cell(UAV.num, 1);                                    % 到达各个点时间树

L_mt = 0;                                                               % 所有无人机航迹之和
totalTime = zeros(UAV.num, 1);                         % 所用时间
totalL = zeros(UAV.num, 1);                                % 每个无人机飞行距离

for i = 1 : UAV.num
      PointNum = UAV.PointNum(i);
      Judge = zeros(1, PointNum+1);
      L_i = zeros(1, PointNum+1);
      Time_i = zeros(1, PointNum+1);
      Angle_i = zeros(1, PointNum+1);
      Traj_i = zeros(1, PointNum+1);
      ProbPoint_i = zeros(1, PointNum+1);
      % 进行检测
      l = 0;
      V = Track.V(i);
      for k = 1 : PointNum
            % 前向检测
            P2 = Track.P{i}(:, k)' ;   % 转置成 1*dim
            if k == 1
                P1 = UAV.S(i, :);
                ZZ = P2 - P1;
                phi0 = atan2(ZZ(2), ZZ(1));   % 偏角检测
                phi1 = phi0;
                d_phi = phi1 - phi0;
                if dim > 2                                % 倾角检测
                    theta0 = atan(ZZ(3) / sqrt(ZZ(1)^2 + ZZ(2)^2));
                    theta1 = theta0;
                    d_theta = theta1 - theta0;
                else
                    d_theta = 0;
                end
            else
                P1 = Track.P{i}(:, k-1)' ;  % 转置成 1*dim
                ZZ = P2 - P1;
                phi1 = atan2(ZZ(2), ZZ(1));
                d_phi = phi1 - phi0;
                phi0 = phi1;
                if dim > 2
                    theta1 = atan(ZZ(3) / sqrt(ZZ(1)^2 + ZZ(2)^2));
                    d_theta = theta1 - theta0;
                    theta0 = theta1;
                else
                    d_theta = 0;
                end
            end

            [across, ~] = CheckThreat(P1, P2, M); % 威胁检测
            Judge(k) = across;

            dl = norm(P1 - P2);
            l = l + dl; % 累加距离
            t = l / V;  % 时间点
            L_i(k) = l;
            Time_i(k) = t;

            if abs(d_phi) > UAV.limt.phi(i, 2)  ||  abs(d_theta) > UAV.limt.theta(i, 2)
                Angle_i(k) = true;
            else
                Angle_i(k) = false;
            end

            if dl < UAV.limt.L(i, 1)
                Traj_i(k) = true;
            else
                Traj_i(k) = false;
            end

            ProbPoint_i(k) = Angle_i(k) | Traj_i(k) | Judge(k); % 问题点

            % 最后一段检测
            if k == PointNum
                P1 = UAV.G(i, :);
                [across, ~] = CheckThreat(P2, P1, M);
                Judge(k+1) = across;

                dl = norm(P1 - P2);
                l = l + dl;
                t = l / V;
                L_i(k+1) = l;
                Time_i(k+1) = t;

                ZZ = P1-P2;
                phi1 = atan2(ZZ(2), ZZ(1));
                d_phi = phi1 - phi0;
                if dim>2
                    theta1 = atan(ZZ(3) / sqrt(ZZ(1)^2 + ZZ(2)^2));
                    d_theta = theta1 - theta0;
                else
                    d_theta = 0;
                end

                if abs(d_phi) > UAV.limt.phi(i, 2)  ||  abs(d_theta) > UAV.limt.theta(i, 2)
                    Angle_i(k+1) = true;
                else
                    Angle_i(k+1) = false;
                end

                if dl < UAV.limt.L(i, 1)
                    Traj_i(k+1) = true;
                else
                    Traj_i(k+1) = false;
                end

                ProbPoint_i(k+1) = Angle_i(k+1) | Traj_i(k+1) | Judge(k+1);

            end
      end
      Threat(i) = {Judge};                % 检测结果（长度比点的数目多一）
      Angle(i) = {Angle_i};               % 转角
      MiniTraj(i) = {Traj_i};              % 航迹片段
      ProbPoint(i) = {ProbPoint_i}; % 问题点

      L(i) = {L_i};                 % 路径长度（累加）
      Time(i) = {Time_i};    % 时间（累加）
      L_mt = L_mt + l;           % 总长度
      totalTime(i) = t;            % 时间
      totalL(i) = l;                   % 长度
end

% 多无人机碰撞检测
d_safe = UAV.ds; % 安全距离
CollideTimes = 0; % 碰撞次数
for i = 2 : UAV.num
    PointNum_i = UAV.PointNum(i);
    for k = 1 : PointNum_i
        P1 = Track.P{i}(:, k)';  % K时刻 i 无人机位置
        t_i = Time{i, 1}(k);                  % K时刻 i 无人机时间
        for j = 1 : i-1
            PointNum_j = UAV.PointNum(j);
            flag = false; %
            % 搜索同一时刻的点
            for kj = 1 : PointNum_j
                if kj ==1
                    t_j_l = 0;
                    P_l =  UAV.S(j, :);
                else
                    t_j_l = Time{j}(kj-1);
                    P_l =  Track.P{j}(:, kj-1)';
                end
                t_j_r = Time{j}(kj);
                P_r =  Track.P{j}(:, kj)';

                if t_i <= t_j_r  &&  t_i >= t_j_l
                    flag = true;
                    P2 =  P_l + (t_i - t_j_l) / (t_j_r - t_j_l) * (P_r - P_l);% K时刻 J 无人机位置
                    % K时刻 j 无人机位置
                end
            end
            % ———————

            if flag  % 查找到P2位置，进行碰撞检测
                collide = CheckCollide(P1, P2, d_safe);
            else     % 没有P2位置，不会碰
                collide = false;
            end
            if collide
                CollideTimes = CollideTimes + 1;
            end
        end
    end
end

% 生成检测报告
report.L_mt = L_mt;                       %总行程之和
report.Threat = Threat;                 %受威胁的航迹点位置
report.AngleProb = Angle;            % 转角不满足的点
report.TrajProb = MiniTraj;           % 不满足最小航迹间隔的点
report.ProbPoint = ProbPoint;     % 有问题的点
report.L = totalL;                            %飞行距离
report.time = totalTime;                %飞行时间
report.col_times = CollideTimes;  %碰撞次数
end

%% 穿越威胁区检测
function [across, across_num] = CheckThreat(P1, P2, M)
    % 威胁区（球或圆区域，不适合圆柱区域）
    O = M(:,1:end-1);  % 圆心
    R = M(:, end);        % 半径

    % 检测线段是否穿过某个障碍区
    total = 0;
    for i = 1 : size(O, 1)
         a = norm(P1 - P2);
         b = norm(P2 - O(i, :));
         c = norm(P1 - O(i, :));
         % P1点是否在圆内
         if c < R(i)
             isHit = true;
         % P2点是否在圆内
         elseif b < R(i)
             isHit = true;
         % P1 P2都不在圆内
         else
             dim = size(O, 2);  % P1： 1*dim 维
             if dim < 3
                 % 平面情况
                 A = P1(2) - P2(2);
                 B = P2(1) - P1(1);
                 C = P1(1)*P2(2) - P2(1)*P1(2);
                 x = O(i, 1);
                 y = O(i, 2);
                 d = abs(A*x + B*y + C) / sqrt(A^2 + B^2);
             else
                 % 空间情况
                 PP = P2 - P1;
                 PO = O(i, :) - P1;
                 d = norm(cross(PP, PO)) / a;
             end
             % 距离判据
             if d >= R(i)
                 isHit = false;  % 不相交
             % 角度判据(距离满足条件时两个角都为锐角时相交)
             elseif d > 0
                 cosP1 = a^2 + c^2 - b^2 / (2*c*a);
                 cosP2 = a^2 + b^2 - c^2 / (2*b*a);
                 if cosP1 > 0  &&  cosP2 > 0
                     isHit = true;
                 else
                     isHit = false;
                 end
             % 两点在外，距离为0情况（共线情况）
             else
                 if a > b  &&  a > c
                    isHit = true;
                 else
                    isHit = false;
                 end
             end
         end

         if isHit
            total = total + 1; %总共碰撞几次
         end
    end

    if total > 0
        across = true;
    else
        across = false;          % 是否穿越禁区
    end
    across_num = total;     % 穿过禁区的个数
end

%% 碰撞检测
function [collide] = CheckCollide(P1, P2, d_safe)
    if norm(P1 - P2) >= d_safe
        collide = false;
    else
        collide = true;
    end
end

repo: Grey-Wolf-Optimizer-for-Path-Planning


text:
function [b, dev, stat] = glm_multilevel(D, Yvarname, Xvarnames, wh_keep)
% Predict Y from X using GLM
%
% :Usage:
% ::
%
%    [b, dev, stat] = glm_multilevel(D, Yvarname, Xvarnames, wh_keep)
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2013 Tor Wager
%
%     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/>.
% ..
%
% :Inputs:
%
%   **D:**
%        a canlab_dataset object
%
%   **Yvarname:**
%        the name of a variable to predict. must be event level
%
%   **Xvarnames:**
%        the name(s) of predictor variables. if multiple, put in
%        cell array. must be event level
%
%   **wh_keep:**
%        a logical vector of 1/0 values
%
% :Outputs:
%   **b:**
%       a vector of coefficient estimates (same as for glmfit())
%
%   **dev:**
%       the deviance of the fit (same as for glmfit())
%
%   **stat:**
%       structure containing stats fields (see glmfit() documentation)
%

[Y, ~, levelY] = get_var(D, Yvarname, wh_keep);
Y = Y';

%% MUST IMPLEMENT GET_VAR FOR MULTIPLE VARS AT AN EVENT LEVEL.  RETURNS A CELL ARRAY FOR EACH PERSON.  DAT BECOMES WARNING STRING, DATCELL IS OF INTEREST.

[X, ~, levelX] = get_var(D, Xvarnames, wh_keep);

if levelY ~= 2 || levelX ~= 2, error('Vars must be event level'); end

n=size(Y,2);
X1 = cell(1,n);
X2 = mean(Y)'; % matrix of 2nd level preds

%{
if isstruct(varargin{1}) % the "alternative format" described above
    xstruct = varargin{1};
    fields = varargin{2};
    for i = 1:n % each subject
        clear myX
        for j = 1:length(fields) %all the fields
            myX(:, j) = xstruct.(fields{j}){i};
            X1{i} = myX;
        end
    getvif(X1{i})
    end

else
%}
    for i = 1:n % each subject
        clear myX
        for j = 1:length(varargin)
            myX(:, j) = varargin{j}{i};
            X1{i} = myX;
        end
    end
%end

stats = glmfit_multilevel(Y, X1, scale(X2, 1), 'weighted', 'noplots', ...
    'names', {'Intrcpt' names{2:end}}, 'beta_names', {'Avg within-ss relationship' ['Effect of Avg. ' names{1}]});

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% General state space model for gps imu integration formulated in n-frame
% either in psi model or phi model. Psi model only appears in some
% papers, Phi model is well documented in books like Jekeli 2000, Titterton
% and Weston, 2000.

% states: IMU position in e frame, IMU velocity in n frame,
% attitude w.r.t n frame, accelerometer and gyro biases, accelerometer and gyro scale factors

% the error states are delta(r x,y,z in earth centered
% "abused" n frame), delta(vx, vy, vz in n frame),
% delta(roll, pitch, yaw) of the imu, accelerometer and gyro bias drift, acc
% scale and gyro scale factor error
% accelerometer and gyro bias are usually modeled as random walk, it can
% also be random constant, This does not make much difference
% acc scale and gyro scale factor are usually modeled as random walk, it can
% also be random constant, This does not make much difference

% The state dynamics are driven by a 18 (accel and gyro white noise,
% bias white noise and scale factor white noise) or 6 (only accel and
% gyro white noise) dimensional white Gaussian
% noise and the observations are corrupted by white Gaussian noise
% test cases: (1) random walk bias and scale facotr error, random constant
% bias and scale factor error,
% (2) covariance update step > sample interval
% (3) different imy types, 3DM GX3-35 AND h764g
% (4) psi and phi model
% There is a very subtle point in implementing a class in matlab. Always
% call a class the same name, do not call it a filter in a function and a
% model at another. It causes strange behaviors.
classdef EKF_filter_nframe < handle
    properties (Hidden)
        type = 'ekf';
        tag  = 'EKF_IMU_GPS_NFRM';  % ID tag
        covDim=9+12; % the dimension of covariance matrix
    end
    % The following properties can be set only by class methods
    properties (SetAccess = private)
        invalidateIMUerrors; % reset IMU errors when they goes beyond expectation
        imuType; % the type of IMU, e.g., MEMS 3DX GM 3-35, HG1700, etc.
        imuErrorModel=3;
        % the IMU error model corresponding to random walk bias and
        % random walk scale factor error,
        % 1 without turn on bias estimates,first order GM bias and random walk scale factor errors
        % 2 with random constant turn on bias estimates, first order GM bias and random walk scale factor errors
        % 3 assumes bias and scale factors are random walks
        % 4 assumes random constant bias and random constant scale factor errors

        % transform from vehicle body frame to IMU sensor frame
        Cb2imu;
        dt; % sampling interval of IMU, unit sec
        % the translation from antenna to IMU sensor frame, i.e., the antenna's
        % position in the s-frame.
        Tant2imu;
        Cen; % rotation from e frame to n frame
        height; % h in geodetic coordinates
        qs2n; % q s(sensor) 2 navigation frame
        Vn; % velocity of the IMU sensor in n frame
        imuErrors; %ba, bg, sa, sg

        imuOrientSIP=7; % imu orientation covariance start index
        imuBiasDriftSIP=10; %imu bias error covariance start index
        imuScaleFactorSIP=16;
        p_k_k; % the covariance of the entire state vector
        mode=1; % 1 for phi model, 2 for psi model
        mechanization=2; % 1 for wander azimuth, 2 for local geodetic
        rvqs2e; % added only for compatibility
        camPose; % added only for compatibility
    end
    methods
        function filter = EKF_filter_nframe(options)
            filter.invalidateIMUerrors=options.InvalidateIMUerrors;
            filter.imuType=options.imutype;
            filter.imuErrorModel=options.imuErrorModel;
            filter.dt=options.dt;
            filter.Cb2imu=options.Cb2imu;
            filter.Tant2imu=filter.Cb2imu*(options.Tant2body-options.Timu2body);
            filter.Vn=options.Vn;
            % initialize states and covariance for imu gps integration,
            % For imu, rs in n, vs in n, q s2n, ba, bg, sa, sg. (s denotes imu)
            % the position of GPS antenna
            inillh=options.inillh_ant;
            % Cb2n initial value given by external reference
            initqbn=options.qb2n;   % H764G INS data, this b refers to the vehicle body frame
            filter.qs2n=quatmult_v001(initqbn,rotro2qr(filter.Cb2imu),2);

            Ce2n0=llh2dcm_v000(inillh(1:2),[0;1]);
            qs2e=quatmult_v001(rotro2qr(Ce2n0), filter.qs2n,1);
            xyz_imu=ecef2geo_v000(inillh,1)-quatrot_v000(qs2e,filter.Tant2imu,0);
            inillh_imu=ecef2geo_v000(xyz_imu, 0);
            filter.height=inillh_imu(3); %elipsoidal height
            filter.Cen=llh2dcm_v000(inillh_imu(1:2),[0;1]);
            filter.imuErrors=options.imuErrors;
            filter.p_k_k=zeros(filter.covDim);
            filter.p_k_k(1:3,1:3)=eye(3)*1^2; %position errors in meter
            filter.p_k_k(4:6,4:6)=eye(3)*0.1^2; %vel error in m/s
            filter.p_k_k(filter.imuOrientSIP+(0:2),filter.imuOrientSIP+(0:2))=...
                diag([options.initAttVar,options.initAttVar,options.initAttVar*2].^2);

            IMU_ERRDEF=imu_err_defs_v000(options.imutype);
            filter.p_k_k(filter.imuBiasDriftSIP+(0:2),filter.imuBiasDriftSIP+(0:2))=4*eye(3)*IMU_ERRDEF.acc_bias_var;% enlarge initial std by 2
            filter.p_k_k(filter.imuBiasDriftSIP+(3:5),filter.imuBiasDriftSIP+(3:5))=4*eye(3)*IMU_ERRDEF.gyro_bias_var;

            filter.p_k_k(filter.imuScaleFactorSIP+(0:2),filter.imuScaleFactorSIP+(0:2))=4*eye(3)*IMU_ERRDEF.acc_scale_var;
            filter.p_k_k(filter.imuScaleFactorSIP+(3:5),filter.imuScaleFactorSIP+(3:5))=4*eye(3)*IMU_ERRDEF.gyro_scale_var;
            filter.mode= options.mode;
        end
        %===============================================================================================
        %-- State transition function
        % propogate state with accelerometer and gyro input at time k-1 to state at k,
        % i.e., X(k|k-1), from state at k-1. U1 contains IMU measurement
        % delta v, and delta angle of gyro or acc and angular rate
        % and 7th row is the previous epoch(k-1) and 8th row is the current epoch(k)
        % imuaccum records the accumulated delta v and delta angle for two speed
        % covariance update
        function imuaccum= ffun_state(filter, imuaccum, U1)
            if(sum(abs(U1(1:3,1)))<1)
                dt1=U1(8,end)-U1(7,1);
                gyroinc=sum(U1(4:6,:),2);
                accinc=sum(U1(1:3,:),2);

                %reset the imu errors if estimates diverges
                IMU_ERRDEF=imu_err_defs_v000(filter.imuType);

                if(filter.invalidateIMUerrors)
                    spuracc=find(abs(filter.imuErrors(1:3))>2*IMU_ERRDEF.initacc_bias_err,1);
                    spurgyro=find(abs(filter.imuErrors(4:6))>2*IMU_ERRDEF.initgyro_bias_err,1);
                    if(~isempty(spuracc)||~isempty(spurgyro))
                        disp(['IMU errors diverge at ' num2str(U1(8,end)) '!']);
                        filter.imuErrors=zeros(12,1);
                    end
                end
                angleinc=gyroinc-filter.imuErrors(4:6)*dt1-diag(gyroinc)*filter.imuErrors(6+(4:6))/1000;
                velinc=accinc-filter.imuErrors(1:3)*dt1-diag(accinc)*filter.imuErrors(6+(1:3))/1000;
                %imu data accumulator for the covariance update
                imuaccum=imuaccum+[velinc;angleinc];
            else
                dt1=U1(8,end)-U1(7,1);
                gyroinc=mean(U1(4:6,:),2);
                accinc=mean(U1(1:3,:),2);

                %reset the imu errors if estimates diverges
                IMU_ERRDEF=imu_err_defs_v000(filter.imuType);

                if(filter.invalidateIMUerrors)
                    spuracc=find(abs(filter.imuErrors(1:3))>2*IMU_ERRDEF.initacc_bias_err,1);
                    spurgyro=find(abs(filter.imuErrors(4:6))>2*IMU_ERRDEF.initgyro_bias_err,1);
                    if(~isempty(spuracc)||~isempty(spurgyro))
                        disp(['IMU errors diverge at ' num2str(U1(8,end)) '!']);
                        filter.imuErrors=zeros(12,1);
                    end
                end
                gyroinc=gyroinc-filter.imuErrors(4:6)-diag(gyroinc)*filter.imuErrors(6+(4:6))/1000;
                accinc=accinc-filter.imuErrors(1:3)-diag(accinc)*filter.imuErrors(6+(1:3))/1000;
                angleinc=gyroinc*dt1;
                velinc=accinc*dt1;
                %imu data accumulator for the covariance update
                imuaccum=imuaccum+[accinc;gyroinc]*dt1;
            end
            %run strapdown with angle and velocity increments
            [filter.qs2n, filter.Vn, filter.Cen, filter.height]=strapdown_Cen_quat_v000(...
                filter.qs2n, filter.Vn, filter.Cen, filter.height, velinc, angleinc, dt1);
        end
        function ffun_covariance(filter, imuaccum, covupt_time, curimutime)
            %propagate the covariance corresponds to states, rs in n, vs in n, \psi or \phi,
            % ba, bg, sa, sg
            covdt=curimutime-covupt_time;
            [STM, Qd]=sys_metric_phipsi_v000(filter.Cen, filter.height, ...
                filter.Vn, filter.qs2n,  imuaccum(1:3)/covdt,imuaccum(4:6)/covdt, ...
                covdt,filter.mode,filter.imuType, filter.imuErrorModel); %use psi implementation with 2
            %              if(filter.imuErrorModel==4)
            %                     [STM, Qd]=sys_metric_phipsi_v002(filter.Cen, filter.height, ...
            %                         filter.Vn, filter.qs2n,  imuaccum(1:3)/covdt,imuaccum(4:6)/covdt, ...
            %                         covdt,filter.mode,filter.imuType); %use psi implementation with 2
            %                     assert(norm(STM1-STM)+norm(Qd1-Qd)<1e-5);
            %              end
            filter.p_k_k=STM*filter.p_k_k*STM'+Qd;  % the covariance of the navigation states and imu error terms
        end
        %==============================================================================================
        %update both the covariance and state
        function correctstates(filter, predict,measure, H,R, ~)
            p_km1_k=filter.p_k_k;
            inno= predict-measure;
            %Kalman
            K=p_km1_k*H'/(H*p_km1_k*H'+R);
            deltaX=K*inno;
            % update covariance
            filter.p_k_k=(eye(filter.covDim)-K*H)*p_km1_k*(eye(filter.covDim)-K*H)'+K*R*K';
            % compute updated state
            [filter.qs2n, filter.Vn, filter.Cen, filter.height]=correctnav_Cen_v000(...
                filter.qs2n, filter.Vn, filter.Cen, filter.height,...
                deltaX(1:filter.imuOrientSIP+2), filter.mode, filter.mechanization);
            filter.imuErrors = filter.imuErrors + deltaX(filter.imuBiasDriftSIP:end);
        end
        % output: position of the IMU in a NED frame anchored at some point or in geodetic coordinates,
        % velocity of the IMU in the body frame, and euler angles of
        % rotation from body frame to the moving NED frame, and covariances
        % of 9 navigation ERROR states in the filter
        function SaveToFile(filter, inillh_ant, preimutime, ffilres)
            vr_c=rotqr2eu('xyz',quatmult_v001(filter.qs2n, rotro2qr(filter.Cb2imu), 0))*180/pi; % body frame w.r.t the moving N frame
            xyz_imu= ecef2geo_v000([asin(-filter.Cen(3,3));asin(-filter.Cen(2,1));filter.height],1);
         %   xyz_ant= xyz_imu + quatrot_v000(quatmult_v001(rotro2qr(filter.Cen),filter.qs2n,1),filter.Tant2imu,0);
            if(~isempty(inillh_ant))
                vr_a=posdiff_v001(xyz_imu, inillh_ant); % position of IMU in the initial N frame
                fwrite(ffilres,[preimutime;vr_a;filter.Vn;vr_c;sqrt(diag(filter.p_k_k(1:9,...
                    1:9)))],'double');
            else
                fwrite(ffilres,[preimutime;ecef2geo_v000(xyz_imu,0);filter.Vn;vr_c;sqrt(diag(filter.p_k_k(1:9,...
                    1:9)))],'double');
            end
        end
        function mag=GetVelocityMag(filter)
            mag=norm(filter.Vn,2);
        end
        function  SetCamAndRvqs2e(filter, options)
            filter.rvqs2e=zeros(10,1);
            Ce2n=filter.Cen;
            heights2n=filter.height;
            Re=6378137/(sqrt(1.0-0.00669437999014*Ce2n(3,3)^2));
            filter.rvqs2e(1:3,1)=-[(Re+heights2n)*Ce2n(3,1);(Re+heights2n)*Ce2n(3,2);(Re*(1-0.00669437999014)+heights2n)*Ce2n(3,3)];
            filter.rvqs2e(4:6,1)=Ce2n'*filter.Vn;
            filter.rvqs2e(7:10,1)=quatmult_v001(rotro2qr(Ce2n),filter.qs2n,1);
            filter.camPose=zeros(7,1);
            filter.camPose(5:7)=options.Cimu2cam*options.Cb2imu*(options.Timu2body-options.Tcam2body);
            filter.camPose(1:4)=rotro2qr(options.Cimu2cam);
        end
    end
end

repo: ekfmonoslam


text:
%% The MTEX Documentation
%
% Documenting a project like MTEX is a challenging and ongoing task for the
% entire community. Therefore we are extremely happy about any spelling
% fixes, examples, theoretical explainations, special use cases, etc. As a
% bonus everybody who contributed to MTEX will automatically appear at the
% <https://github.com/mtex-toolbox/mtex/graphs/contributors contributors
% page>.
%
%% Making Changes Online
%
% The easiest way to suggest changes to the MTEX documentation is online
% visa Github. This is done as follows:
%
% # sign into <https://github.com GitHub>
% # open the help page you want to change, e.g., this one
% # click
% <https://github.com/mtex-toolbox/mtex/blob/develop/doc/GeneralConcepts/Contribute2Doc.m
% edit page> - this leads you to the file on <https://github.com GitHub>
% # click the little pencil in the top right corner which allows you to
% edit the script online
% # enter of copy and paste your changes into the site
% # click at the bottom *suggest changes*
% # click create *pull request*
%
%% Making Changes Within Matlab
%
% Since all documentation pages are essentially MTEX scripts they can be
% opened in Matlab by the command
%
%   edit hepPageName
%
% This allows you to run through all the commands and investiage the effect
% of different parameters or different input files or add further analysis
% steps.
%
% If you like the changes you made simply copy and paste the changes into
% the GitHub page as explained above or send it to one of the developers by
% email.
%
%% Publish command
%
% You may can a visual preview of your altered documentation page by using
% the matlab function
% <https://de.mathworks.com/help/matlab/matlab_prog/publishing-matlab-code.html
% publish>
%
%   publish filename
%
% which creates an folder, named 'html', in your current directory where
% you can find the generated html page.
%


text:
function lev = QLmaxlev(sizeX, filtername)
%-----------------------------------------------------------------------------
% QLmaxlev
% This function determines the maximum level possible for the QL-schemes.
%
% Syntax: lev = QLmaxlev(sizeX, filtername)
%
% QLmaxlev is a utility for the lifting scheme decomposition. It helps one to
% avoid silly values for the maximum level in the lifting scheme decomposition.
%
% lev is the integer outcome.
%
% sizeX is an integer row vector of dimension 2. Usually it will be the size of
% an image.
%
% filtername must be a string from the set:
% {Neville2, Neville4, Neville6, Neville8, MaxMin, MinMin, MaxMax}
%
% See also: QLiftDec2, QLiftRec2, wmaxlev
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: February 17, 2003.
% (c) 1999-2003 Stichting CWI, Amsterdam
%-----------------------------------------------------------------------------
%Firstly, check input data
if nargin ~= 2
  error(' QLmaxlev - number of arguments should be 2 ');
end
if ~ischar(filtername)
  error(' QLmaxlev - format of filtername should be character array ');
end
%
switch lower(filtername)
 case 'neville2'
   diam =  3;
 case 'neville4'
   diam =  5;
 case 'neville6'
   diam =  7;
 case 'neville8'
   diam = 11;
 case 'maxmin'
   diam =  3;
 case 'maxmax'
   diam =  3;
 case 'minmax'
   diam =  3;
 case 'minmin'
   diam =  3;
 otherwise
   error(' QLmaxlev - unknown filter ')
end
%
if isempty(sizeX)
  lev = [];
  return;
else
  lev=0;
  n=min(sizeX);
  z=2*diam-1;
  while z<=n
    lev=lev+2;
    z=2*z-1;
  end
end
%-----------------------------------------------------------------------------


text:
function gipl_write_volume(I,fname,scales)
% Function for writing Guys Image Processing Lab (Gipl) files.
%
% gipl_write_volume(volume, filename, voxelsize in mm)
%
% examples:
% I=uint8(rand(64,64,64)*256);
%
% 1: gipl_write_volume(I);
% 2: gipl_write_volume(I,'random.gipl',[2 2 2];

image_types=struct('double',65,'float',65,'single',64,'logical',1,'int8',7,'uint8',8,'int16',15,'uint16',16,'int32',32,'uint32',31);
bit_lengths=struct('double',64,'float',64,'single',32,'logical',1,'int8',8,'uint8',8,'int16',16,'uint16',16,'int32',32,'uint32',32);

% Filename
    if((nargin<2)||strcmp(fname,'')),
        [filename, pathname] = uiputfile('*.gipl', 'Write gipl-file');
        fname = [pathname filename];
    end
% Sizes
    sizes=size(I);
    while(length(sizes)<4), sizes=[sizes 1]; end
% Scales
    if(exist('scales','var')==0), scales=ones(1,length(sizes)); end;
    while(length(scales)<4), scales=[scales 0]; end
% Offset
    offset=256;
% Image Type
    image_type=getfield(image_types,class(I));
% File size
    fsize=offset+prod(sizes)*getfield(bit_lengths,class(I))/8;
% Patient
    patient='Generated by Matlab';
    while(length(patient)<80), patient=[patient ' ']; end
% Matrix
    matrix=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
% orientation
    orientation=0;
% voxel_min
    voxmin=min(I(:));
% voxel_max
    voxmax=max(I(:));
% origing
    origin=[0 0 0 0];
% RescaleIntercept
    RescaleIntercept=0;
% RescaleSlope
    RescaleSlope=1;
% interslicegap
    interslicegap=0;
% user_def2
    user_def2=0;
% par2
    par2=0;
% magic_number
    magic_number=4026526128;

trans_type{1}='binary'; trans_type{7}='char'; trans_type{8}='uchar'; trans_type{15}='short';
trans_type{16}='ushort'; trans_type{31}='uint'; trans_type{32}='int'; trans_type{64}='float';
trans_type{65}='double'; trans_type{144}='C_short'; trans_type{160}='C_int'; trans_type{192}='C_float';
trans_type{193}='C_double'; trans_type{200}='surface'; trans_type{201}='polygon';

trans_orien{0+1}='UNDEFINED'; trans_orien{1+1}='UNDEFINED_PROJECTION'; trans_orien{2+1}='AP_PROJECTION';
trans_orien{3+1}='LATERAL_PROJECTION'; trans_orien{4+1}='OBLIQUE_PROJECTION';  trans_orien{8+1}='UNDEFINED_TOMO';
trans_orien{9+1}='AXIAL'; trans_orien{10+1}='CORONAL'; trans_orien{11+1}='SAGITTAL'; trans_orien{12+1}='OBLIQUE_TOMO';

disp(['filename : ' num2str(fname)]);
disp(['filesize : ' num2str(fsize)]);
disp(['sizes : ' num2str(sizes)]);
disp(['scales : ' num2str(scales)]);
disp(['image_type : ' num2str(image_type) ' - ' trans_type{image_type}]);
disp(['patient : ' patient]);
disp(['matrix : ' num2str(matrix)]);
disp(['orientation : ' num2str(orientation) ' - ' trans_orien{orientation+1}]);
disp(['voxel min : ' num2str(voxmin)]);
disp(['voxel max : ' num2str(voxmax)]);
disp(['origing : ' num2str(origin)]);
disp(['RescaleIntercept : ' num2str(RescaleIntercept)]);
disp(['RescaleSlope : ' num2str(RescaleSlope)]);
disp(['interslicegap : ' num2str(interslicegap)]);
disp(['user_def2 : ' num2str(user_def2)]);
disp(['par2 : ' num2str(par2)]);
disp(['offset : ' num2str(offset)]);

fout=fopen(fname,'wb','ieee-be');
fwrite(fout,uint16(sizes),'ushort'); % 4
fwrite(fout,uint16(image_type),'ushort'); % 1
fwrite(fout,single(scales),'float'); % 4
fwrite(fout,patient,'char'); % 80
fwrite(fout,single(matrix),'float'); % 20
fwrite(fout,uint8(orientation),'uint8'); % 1
fwrite(fout,uint8(par2),'uint8'); % 1
fwrite(fout,double(voxmin),'double'); % 1
fwrite(fout,double(voxmax),'double'); % 1
fwrite(fout,double(origin),'double'); % 4
fwrite(fout,single(RescaleIntercept),'float'); % 1
fwrite(fout,single(RescaleSlope),'float'); % 1
fwrite(fout,single(interslicegap),'float'); % 1
fwrite(fout,single(user_def2),'float'); % 1
fwrite(fout,uint32(magic_number),'uint'); % 1
fwrite(fout,I, trans_type{image_type});
fclose('all');

author: tomdoel

repo: pulmonarytoolkit


text:
%   X: data matrix, each row is one observation, each column is one feature
%   d: reduced dimension
%   type: type of kernel, can be 'simple', 'poly', or 'gaussian'
%   para (input): parameter for computing the 'poly' kernel, for 'simple'
%       and 'gaussian' it will be ignored
%   Y: dimensionanlity-reduced data
%   eigVector: eigen-vector, will later be used for pre-image
%       reconstruction
%   para (output): automatically selected Gaussian kernel parameter

%   Copyright by Quan Wang, 2011/05/10
%   Please cite: Quan Wang. Kernel Principal Component Analysis and its
%   Applications in Face Recognition and Active Shape Models.
%   arXiv:1207.3538 [cs.CV], 2012.

function [Y, eigVector, para]=kPCA(X,d,type,para)

%% check input
if ( strcmp(type,'simple') || strcmp(type,'poly') || ...
        strcmp(type,'gaussian') ) == 0
    Y=[];
    eigVector=[];
    para=[];
    fprintf(['\nError: Kernel type ' type ' is not supported. \n']);
    return;
end

%% parameters
N=size(X,1);
if strcmp(type,'gaussian')
    DIST=zeros(N,N);
    for k=1:size(X,2)
        [a,b]=meshgrid(X(:,k));
        DIST=DIST+(a-b).^2;
    end
    DIST=DIST.^0.5;
    DIST=DIST+diag(ones(N,1)*inf);
    para=10*median(min(DIST));
end

%% kernel PCA
K0=kernel(X,type,para);
oneN=ones(N,N)/N;
K=K0-oneN*K0-K0*oneN+oneN*K0*oneN;

%% eigenvalue analysis
[V,D]=eig(K/N);
eigValue=diag(D);
% eigValue=eigValue(1:min(size(X)));
[eigValue,IX]=sort(eigValue,'descend');
eigVector=V(:,IX);

%% normailization
norm_eigVector=sqrt(sum(eigVector.^2));
eigVector=eigVector./repmat(norm_eigVector,size(eigVector,1),1);

%% dimensionality reduction
Y=K0*eigVector(:,1:d);


text:
function [UTC1,UTC2]=UT12UTC(UT11,UT12,deltaUTCUT1)
%%UT12UTC Convert from UT1, which is a nonuniform timescale based on the
%         rotation rate of the Earth, to Coordinated Universal Time (UTC),
%         which is the "standard" time that pretty much every country uses
%         (the time zone is that of the United Kingdom) including
%         leapseconds. Due to the use of leapseconds, UT1 is always within
%         0.9 seconds of UTC.
%
%INPUTS: Jul1,Jul2 Two parts of a Julian date given in UTC. 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.
%      deltaUTCUT1 An optional parameter specifying the offset between UTC
%                  and UT1 in seconds. If this parameter is omitted, then
%                  the value of the function getEOP will be used and
%                  iterations will be performed, because the correct
%                  tabulated value is given in UTC, not UT1.
%
%OUTPUTS: Jul1,Jul2 Two parts of a pseudo-Julian date in UTC with units of
%                   Julian days.
%
%This essentially just adds deltaUTCUT1 from the given UTC date, or
%looks up deltaUTCUT1 using the getEOP function and then adds
%deltaUTCUT1. The "looking up" step involves an interation, since the value
%of deltaUTCUT1 is tabulated by UTC, which is unknown. To know the peroper
%amount to subtract from the pseudo-Julian day, the function leapSecsOnDay
%is used to determine whether the day contains a leapsecond. As that
%function is also tabulated by UTC, an iteration is necessary even if
%deltaUTCUT1 is given. Unlike the function iauUtcut1 in the International
%Astronomical Union's Standards of Fundamental Astronomy library, this does
%not try to deal with deltaUTCUT1 being off by a second near a leapsecond
%as it was noticed that the "corrections" made in the library make the pair
%of UTC2UT1 and UT12UTC inconsistent by a second when it is near a
%leapsecond.
%
%Many temporal coordinate systems standards 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.
%
%April 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%The initial estimate of UTC is just UT1. This will be within 0.9 seconds
%of the correct value.
UTC1=UT11;
UTC2=UT12;

for curIter=1:2
    if(nargin<3)
        [~,~,deltaUTCUT1]=getEOP(UTC1,UTC2);
    end

    %Find the number of seconds in the day.
    [year,month,day]=UTC2Cal(UTC1,UTC2,true);
    numSecInDay=24*60*(60+leapSecsOnDay(year,month,day));

    %Add preserving precision.
    if(UTC1<UTC2)
        UTC2=UT11+deltaUTCUT1/numSecInDay;
        UTC1=UT12;
    else
        UTC2=UT12+deltaUTCUT1/numSecInDay;
        UTC1=UT11;
    end

    %Put the fractional part into UTC2 and the integer part into UTC1.
    frac1=UTC1-fix(UTC1);
    frac2=UTC2-fix(UTC2);
    sumVal=frac1+frac2;
    sumInt=fix(sumVal);
    sumFrac=sumVal-sumInt;
    UTC1=fix(UTC2)+fix(UTC1)+sumInt;
    UTC2=sumFrac;
end
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:
function [Afp,Bfp]=freeprecess(T,T1,T2,df)
%
%	Function simulates free precession and decay
%	over a time interval T, given relaxation times T1 and T2
%	and off-resonance df.  Times in ms, off-resonance in Hz.

phi = 2*pi*df*(T/1000);	% Off-resonance precession, radians.

% Relaxation exponentials
E1 = exp(-T/T1);
E2 = exp(-T/T2);

% Decay and phase due to off-resonance
Afp = [E2 0 0;0 E2 0;0 0 E1]*zrot(phi); % Mathieu, check order of matrix multiplication

% Regrowth
Bfp = [0 0 1-E1]';

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function test_suite = test_bayesoptimization3
    %   Run specific demo, save values and compare the results to the expected.
    %   Works for both xUnit Test Framework package by Steve Eddins and for
    %   the built-in Unit Testing Framework (as of Matlab version 2013b).
    %
    %   See also
    %     TEST_ALL, DEMO_BAYEESOPTIMIZATION
    %
    %  Copyright (c) 2015 Jarno Vanhatalo
    %  Copyright (c) 2017 Eero Siivola

    % This software is distributed under the GNU General Public
    % License (version 3 or later); please refer to the file
    % License.txt, included with the software, for details.

    % Check which package to use
    if exist('initTestSuite', 'file')
    initTestSuite_path = which('initTestSuite');
    else
    initTestSuite_path = '';
    end
    if ~isempty(initTestSuite_path) && ...
     exist([initTestSuite_path(1:end-15) 'runtests'], 'file')
    % xUnit package
    initTestSuite;
    else
    % Built-in package
    % Use all functions except the @setup
    tests = localfunctions;
    tests = tests(~cellfun(@(x)strcmp(func2str(x),'setup'),tests));
    test_suite = functiontests(tests);
    end

end

function testRunDemo(testCase)
  % Run the correspondin demo and save the values. Note this test has to
  % be run at lest once before the other test may succeed.

  % Suppress an expected warning
  Sw = warning('off','MATLAB:nearlySingularMatrix');
  rundemo(getName())
  warning(Sw); % Return the original warning visibility status
end

function testValues(testCase)
  % Test predictive mean and variance
  verifyVarsEqual(testCase, getName(), ...
    {'x'}, ...
    'RelTolElement', 0.01, 'RelTolRange', 0.02)
end

function testEI(testCase)
  % Test predictive mean and variance
  verifyVarsEqual(testCase, getName(), ...
    {'EI'}, ...
    'RelTolElement', 0.01, 'RelTolRange', 0.02)
end

% ------------------------
%     Helper functions
% ------------------------

function testCase = setup
  % Helper function to suply empty array into variable testCase as an
  % argument for each test function, if using xUnit package. Not to be
  % used with built-in test framework.
  testCase = [];
end

function name = getName
  % Helperfunction that returns the name of the demo, e.g. 'binomial1'.
  name = mfilename;
  name = name(6:end);
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function sdlfmKernDisplay(kern, spacing)

% SDLFMKERNDISPLAY Display parameters of the SDLFM kernel.
% FORMAT
% DESC displays the parameters of the switching dynamical latent force
% model kernel and the kernel type to the console.
% ARG kern : the kernel to display.
%
% FORMAT does the same as above, but indents the display according
% to the amount specified.
% ARG kern : the kernel to display.
% ARG spacing : how many spaces to indent the display of the kernel by.
%
% SEEALSO : sdlfmKernParamInit, modelDisplay, kernDisplay
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% MODIFICATIONS : Mauricio A. Alvarez, 2010.

% KERN

if nargin > 1
  spacing = repmat(32, 1, spacing);
else
  spacing = [];
end
spacing = char(spacing);
fprintf(spacing);
fprintf('SDLFM mass: %2.4f\n', kern.mass)
fprintf(spacing);
fprintf('SDLFM spring: %2.4f\n', kern.spring)
fprintf(spacing);
fprintf('SDLFM damper: %2.4f\n', kern.damper)
for i=1:kern.nlfPerInt
    for j=1:kern.nIntervals
        fprintf(spacing);
        fprintf('SDLFM inverse width %d interval %d: %2.4f (length scale %2.4f)\n', ...
           i,j, kern.inverseWidth(i,j), 1/sqrt(kern.inverseWidth(i,j)));
    end
end
for j=1:kern.nIntervals
    fprintf(spacing);
    fprintf('SDLFM switching point interval %d: %2.4f \n', ...
        j, kern.switchingTimes(j));
end
for i=1:kern.nlfPerInt
    for j=1:kern.nIntervals
        fprintf(spacing);
        fprintf('SDLFM sensitivity %d interval %d: %2.4f \n', ...
           i,j, kern.sensitivity(i,j));
    end
end
fprintf(spacing);
fprintf('System Characteristics:\n')
fprintf(spacing);
fprintf('SDLFM omega: %2.4f\n', kern.omega)
fprintf(spacing);
fprintf('SDLFM alpha: %2.4f\n', kern.alpha)
fprintf(spacing);
fprintf('SDLFM gamma: %2.4f\n', kern.gamma)
fprintf(spacing);
fprintf('SDLFM Damping Ratio: %2.4f\n', kern.zeta)
fprintf(spacing);
fprintf('SDLFM Undamped Natural Frequency: %2.4f\n', kern.omega_0)

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function out = subsref(x, s)
	%SUBSREF Subscripted reference for TTeMPS.
	%
	%   Examples:
	%   x([2,3,1]) computes and returns the element (2,3,1) of x.
	%   x( ind ) computes and returns all elements x(ind(i,:)) for i = 1:size(ind,1).
	%   x(:) returns the vectorization of full(x) (careful!)
    %   x{i} returns the i-th core of x, shorthand for x.U{i}
	%
	%   See also TTEMPS.

    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

	switch s(1).type

	case '.'

		prop = properties('TTeMPS');
		if( any(strcmp(s(1).subs, prop) ) ),
		    out = builtin('subsref', x, s);
		else
			ll = length(prop);
			proplist = repmat({', '}, 2*ll-1, 1);
			proplist(1:2:end) = prop;
			proplist = cat(2,proplist{:});
			error(['Object TTeMPS does not have field ' s(1).subs ...
			'. The following fields are available: ' proplist '.']);
		end

	case '()'
		% x(:)
		ind = s(1).subs{1};
		if(length(ind) == 1) && (ind == ':')
			out = full(x);
			out = out(:);
			% e.g. x( [1,2,3; 4 5 6; 7 8 9; 3 5 3]) for d = 3 tensor
		elseif(size(ind,2) == x.order)

			r = x.rank;
			%out = zeros(size(ind,1), 1);

			%C = cell(1,x.order);
			%for i=1:x.order
			%	C{i} = permute( x.U{i}, [1 3 2]);
			%end
			%for i = 1:size(ind,1)
			%	p = C{1}(:,:,ind(i,1));
			%	for j = 2:size(ind,2)
			%		p = p * C{j}(:,:,ind(i,j));
			%	end
			%	out(i) = p;
			%end
			n = x.size;

			C = cell(1,x.order);
			for i=1:x.order
				C{i} = permute( x.U{i}, [1 3 2]);
				C{i} = unfold( C{i}, 'right');
			end

            out = TTeMPS.subsref_mex( n, r, transpose(ind), C);

			%for i = 1:size(ind,1)
			%	p = C{1}(:, (ind(i,1)-1)*r(2)+1:ind(i,1)*r(2));
			%	for j = 2:size(ind,2)
			%		p = p * C{j}(: , (ind(i,j)-1)*r(j+1)+1:ind(i,j)*r(j+1));
			%	end
			%	out(i) = p;
			%end

		else
			error('Number of indices does not match order of TTeMPS tensor.');
		end

	case '{}'

		if(length(s(1).subs) ~= 1 || ~isnumeric(s(1).subs{1}) || ...
			s(1).subs{1} <= 0)
			error('{} only takes one positive integer.');
		end

		ii = s(1).subs{1};
		if(ii > x.order)
			error('Index exceeds number of dimensions');
		end
		out = builtin('subsref', x.U, s);

	end

author: NicolasBoumal

repo: manopt


text:
function jed = ymdf_to_jed_bahai ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_JED_BAHAI converts a Bahai YMDF date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Algorithm E,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, pages 323-324.
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, real JED, the corresponding Julian Ephemeris Date.
%

%
%  Convert the calendar date to a computational date.
%
  y_prime = y + 6560 - floor ( ( 39 - m ) / 20 );
  m_prime = mod ( m, 20 );
  d_prime = d - 1;
%
%  Convert the computational date to a JED.
%
  j1 = floor ( ( 1461 * y_prime ) / 4 );

  j2 = 19 * m_prime;

  g = floor ( ( y_prime + 184 ) / 100 );
  g = floor ( ( 3 * g ) / 4 ) - 50;
  jed = j1 + j2 + d_prime - 1412 - g - 0.5;
  jed = jed + f;

  return
end


text:
% Laplacian Smoothing Transform for face recognition

function [gammak,lambda]=lst_basis(M,N,k)
% By Suicheng Gu at Peking University, Jan. 5, 2009
% Email: gusuch@gmail.com

% Output: v: basis. vv: k minimal eigenvalues
% Input:
% N: width of input image, or columns of a matrix
% M: height of input image, or rows of a matrix
% k: Number of low frequency features

%************ An Example ************
% To extract 30 LST low frequency coefficients of an grey image:
%  Im=rand(50,40); [M,N]=size(Im);
%  v=lst_basis(M,N,30);
%  temp=zeros(1,M*N);temp(:)=Im(:);
%  Im_feature=temp*v;
%*********** End ******************

W = sparse([],[],[],M*N,M*N,4*M*N-2*M-2*N);
for x=0:N-2;
    for y=0:M-1;
        W(x*M+y+1,x*M+M+y+1)=1;
        W(x*M+M+y+1,x*M+y+1)=1;
    end;
end;
for x=0:N-1;
    for y=0:M-2;
        W(x*M+y+1,x*M+y+2)=1;
        W(x*M+y+2,x*M+y+1)=1;
    end;
end;
L=sparse(diag(sum(W)))-W;
options.disp = 0; options.isreal = 1; options.issym = 1;
[v,lambda]=eigs(L,k,0,options);
[vv,index] = sort(diag(lambda));
gammak = v(:,index);


text:
 function xs = eml_inc_em(x, Gb, yi, ci, ri, varargin)
%|function xs = eml_inc_em(x, Gb, yi, ci, ri, varargin)
%| E-ML-INC-EM algorithm for image reconstruction from Poisson emission data
%| (incremental EM algorithm)
%| model: Y_i ~ Poisson(c_i [G x]_i + r_i)
%| in
%|	x	[np,1]		initial estimate
%|	Gb	[nd,np]		Gblock object (see eml_osem_test.m)
%|	yi,ci,ri [nb,na]	see em_fbp.m (for model too)
%|
%| options
%|	'niter'	(#)		# of iterations
%|	'isave'	[]		list of iterations to archive
%|				(default: [] 'last)
%|	chat	(0 | 1)		verbosity
%|	pixmax	(value)		upper constraint for pixel values
%|	hds	(value)		'1' or '3' hidden data space
%|	os	(#)		how many "warmup" osem iterations
%|
%| out
%|	xs [np nsave]		estimates each (saved) iteration
%|
%| Copyright 2004-3-20, Jeff Fessler, University of Michigan

if nargin < 3, ir_usage, end

% options
arg.chat = 0;
arg.pixmax = inf;
arg.hds = 1;
arg.niter = 1;
arg.isave = 'last';
arg.os = 0;
arg.userfun = [];

arg = vararg_pair(arg, varargin);
arg.isave = iter_saver(arg.isave, arg.niter);

if ~isvar('ci') || isempty(ci)
	ci = ones(size(yi));
end
if ~isvar('ri') || isempty(ri)
	ri = zeros(size(yi));
end

Gb = block_ob(Gb, 'ensure'); % make it a block object (if not already)
nblock = block_ob(Gb, 'n');
starts = subset_start(nblock);

eml_check(yi, ci, ri, 'os', nblock);
[nb na] = size(yi);

ticker(mfilename, 1, arg.niter)

% precompute hidden data space factor
gam = eml_hds(Gb, ci, ri, arg.hds);
if arg.chat, printf('hds = %g', gam), end

% system (block) sensitivities
np = length(x);
asum = zeros(np, 1);
precon = zeros(np, nblock); % classic OSEM preconditioner
for iset=1:nblock
	ticker([mfilename ' : precon'], iset, nblock)
	istart = starts(iset);
	ia = istart:nblock:na;
	asum_m = Gb{istart}' * col(ci(:,ia));
	asum = asum + asum_m;
	asum_m(asum_m == 0) = Inf; % avoid divide by 0
	precon(:, iset) = 1 ./ asum_m;
end, clear asum_m
% asum = Gb' * ci(:);

x = max(x,0);
x = min(x,arg.pixmax);
xs = zeros(numel(x), length(arg.isave));
if any(arg.isave == 0)
	xs(:,find(arg.isave == 0)) = x;
end

%
% precompute sufficient statistics,
% possibly using osem warmup iterations
%
ff = zeros(np,nblock);
if arg.os == 0
	for iblock=1:nblock
		ff(:,iblock) = compute_fm(x, Gb, yi, ci, ri, gam, iblock);
	end

else
	for iter = 1:arg.os
		for iset = 1:nblock
			ticker([mfilename ': osem'], [iter iset], [arg.os nblock])
			iblock = starts(iset);
			ia = iblock:nblock:na;

			tmp = compute_fm(x, Gb, yi, ci, ri, gam, iblock);

			pre = precon(:,min(iset,ncol(precon)));	% 1 or iset
			x = max(pre .* tmp - gam, 0);

			if iter == arg.os
				ff(:,iblock) = tmp;
			end
		end

		if any(arg.isave == iter)
			xs(:, arg.isave == iter) = x;
		end
	end
end
fsum = sum(ff, 2);

%
% loop over incremental EM iterations
%

for iter = (2+arg.os):arg.niter

	%
	% loop over subsets
	%
	for iset=1:nblock
		ticker(mfilename, [iter iset], [arg.niter nblock])
		iblock = starts(iset);

		x = fsum ./ asum - gam;
		x = max(x,0);
		x = min(x,arg.pixmax);

		fsum = fsum - ff(:,iblock);
		ff(:,iblock) = compute_fm(x, Gb, yi, ci, ri, gam, iblock);
		fsum = fsum + ff(:,iblock);
	end

	if arg.chat, printf('Range %g %g', min(x), max(x)), end
	if any(arg.isave == iter)
		xs(:, arg.isave == iter) = x;
	end
end

%
% compute_fm()
% compute sufficient statistics
%
function fm = compute_fm(x, Gb, yi, ci, ri, gam, iblock)
eterm = eml_eterm(x, Gb, yi, ci, ri, iblock);
fm = (x + gam) .* eterm;

author: JeffFessler

repo: mirt


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = INVERSEKINEMATIC_KUKA_KR_30_L16_2(robot, T)
%   Solves the inverse kinematic problem for the KUKA KR5 ARC robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=INVERSEKINEMATIC_KUKA_KR_30_L16_2 returns 8 possible solutions, thus,
%   Q is a 6x8 matrix where each column stores 6 feasible joint values.
%
%
%   Example code:
%
%   robot=load_robot('kuka', 'KR_30_L16_2');
%   q = [0 0 0 0 0 0];
%   T = directkinematic(robot, q);
%   %Call the inversekinematic for this robot
%   qinv = inversekinematic(robot, T);
%   check that all of them are feasible solutions!
%   and every Ti equals T
%   for i=1:8,
%        Ti = directkinematic(robot, qinv(:,i))
%   end
%	See also DIRECTKINEMATIC.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function q = inversekinematic_kuka_kr20_3(robot, T)

%initialize q,
%eight possible solutions are generally feasible
q=zeros(6,8);

% %Evaluate the parameters
% theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
L6=abs(d(6));

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Compute the position of the wrist, being W the Z component of the end effector's system
W = T(1:3,3);

% Pm: wrist position
Pm = [Px Py Pz]' - L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution
q1=atan2(Pm(2), Pm(1));

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);
%the other possible solution is q1 + pi
q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);
%solver for q3 for both cases
q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];

%leave only the real part of the solutions
% q=real(q);
% q = [pi/4 -pi/6 pi/4 0.1 0.1 0.1];
% T = directkinematic(robot,q)
% T =
%    -0.3864   -0.6500    0.6544    0.4980
%    -0.1062    0.7361    0.6685    0.5003
%    -0.9162    0.1888   -0.3535    1.0008
%          0         0         0    1.0000
%q =
%    0.7854   -0.5236    0.7854    0.1000    0.1000    0.1000

%Note that in this robot, the joint q3 has a non-simmetrical range. In this
%case, the joint ranges from 60 deg to -219 deg, thus, the typical normalizing
%step is avoided in this angle (the next line is commented). When solving
%for the orientation, the solutions are normalized to the [-pi, pi] range
%only for the theta4, theta5 and theta6 joints.

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));

% solve for the last three joints
% for any of the possible combinations (theta1, theta2, theta3)
for i=1:2:size(q,2),
    qtemp = solve_spherical_wrist2(robot, q(:,i), T, 1,'geometric'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    qtemp = solve_spherical_wrist2(robot, q(:,i), T, -1, 'geometric'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry
L2=abs(a(2));
L3=abs(d(4));
A2 = (a(3));

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%The inverse kinematic problem can be solved as in the IRB 140 (for example)

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = atan2(-p1(2), p1(1));
gamma = real(acos((L2^2+r^2-L4^2)/(2*r*L2)));

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta3
q2(1) = pi/2 - beta - gamma; %elbow up
q2(2) = pi/2 - beta + gamma; %elbow down

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry
L2=abs(a(2));
L3=abs(d(4));

A2 = (a(3));

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%the angle phi is fixed
phi=acos((A2^2+L4^2-L3^2)/(2*A2*L4));

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = real(acos((L2^2 + L4^2 - r^2)/(2*L2*L4)));

%return two possible solutions
%elbow up and elbow down solutions
%the order here is important
q3(1) = pi - phi - beta;
q3(2) = pi - phi + beta;

repo: ARTE