text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes the angular displacements of a damped pendulum under
%           gravity
%
% PURPOSE: to illustrate that changing the pendulum bob's radius results in
%          significant dynamical differences in the pendulum system
%
% Author: N.A. Battista
% Date: 12/13/2019
% Institution: The College of New Jersey
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function damped_pendulum_diff_r()

%
% Physical Parameters
%
m = 100;    % Mass of Pendulum bob
g=9.80665;  % Gravitational Acceleration;
L=0.2;      % Length of Pendulum

%
% Radii of Bobs (here - all different)
%
rVec = 100*[0.0025 0.005 0.01 0.015 0.02 0.025];

%
% Damping Parameter (here - all the same)
%
bVec = 50*ones(1,6);

%
% Temporal Parameters
%
t=0;        % initial time
tFinal = 24; % final-time
dt = 1e-3;  % time-step
ct=1;       % counter
tVec(ct)=t; % time vector storage

%
% Moment of Inertia
%
IVec = m*(0.5*rVec.^2+L^2);

%
% Initial Angular Displacement
%
ang = -(pi/2-pi/5)*ones(1,6);
phi = zeros( size( ang ) );

%
% Perform the time-stepping, e.g., solve the ODEs
%
while t<=tFinal

    % increment the storage counter
    ct = ct+1;

    % increment time
    t = t+dt;

    % save current time
    tVec(ct) = t;

    % solve 1st order ODE for angular displacement
    ang(ct,:) = ang(ct-1,:) + dt * ( phi(ct-1,:) );

    % solve 1st order ODE for angular velocity
    phi(ct,:) = phi(ct-1,:) + dt * ( -m*g*L ./IVec .* sin( ang(ct-1,:) ) - bVec.*phi(ct-1,:)./IVec );

end

%
% Defining Colors for Plotting Data
%
color2 = [0.635 0.078 0.184];  % dark red
color5 = [0.9 0.525 0.098];    % orange
color6 = [0.929 0.840 0.1250]; % dark yellow
color7 = [0.7 0.7 0];          % 'ugly' green
color13 = [0 0.5 0.5];             % green

%
% Note discrepancy is a bit off due to small angle approximation:
% sin(theta) ~ theta gives below period of oscillation
%
T1 = 2*pi*sqrt( IVec(2) / (m*g*L) );
%
% Actually computed value of period is:
T1 = 1.93;

%
% Make Figure of Angular Displacement vs. Non-dimensional Time (# of oscillations of 'r' case)
%
lw=6;
fs=18;
plot(tVec/T1,ang(:,1),'k-','LineWidth',lw); hold on;
plot(tVec/T1,ang(:,2),'b-','LineWidth',lw,'Color',color2); hold on;
plot(tVec/T1,ang(:,3),'g-','LineWidth',lw,'Color',color5); hold on;
plot(tVec/T1,ang(:,4),'r-','LineWidth',lw,'Color',color6); hold on;
plot(tVec/T1,ang(:,5),'c-','LineWidth',lw,'Color',color7); hold on;
plot(tVec/T1,ang(:,6),'c-','LineWidth',lw,'Color',color13); hold on;
plot(tVec/T1,ang(:,5),'c-','LineWidth',lw,'Color',color7); hold on;
plot(tVec/T1,ang(:,4),'r-','LineWidth',lw,'Color',color6); hold on;
plot(tVec/T1,ang(:,3),'g-','LineWidth',lw,'Color',color5); hold on;
plot(tVec/T1,ang(:,2),'b-','LineWidth',lw,'Color',color2); hold on;
plot(tVec/T1,ang(:,1),'k-','LineWidth',lw); hold on;
xlabel('Non-Dimensional Time (# periods of r-case)');
ylabel('Angular Displacement (Radians)');
axis([0 12 -1 1.5]);
leg=legend('$\frac{1}{2}$r','r','2r','3r','4r','5r','Orientation','horizontal','interpreter','latex');
set(gca,'FontSize',fs);
set(leg,'FontSize',fs+2);
leg.NumColumns=3;

author: nickabattista


text:
function trans = createRotationOz(varargin)
%CREATEROTATIONOZ Create the 4x4 matrix of a 3D rotation around z-axis.
%
%   TRANS = createRotationOz(THETA);
%   Returns the transform matrix corresponding to a rotation by the angle
%   THETA (in radians) around the Oz axis. A rotation by an angle of PI/2
%   would transform the vector [1 0 0] into the vector [0 1 0].
%
%   The returned matrix has the form:
%   [cos(THETA) -sin(THETA)  0  0]
%   [sin(THETA)  cos(THETA)  0  0]
%   [    0           0       1  0]
%   [    0           0       0  1]
%
%   TRANS = createRotationOz(ORIGIN, THETA);
%   TRANS = createRotationOz(X0, Y0, Z0, THETA);
%   Also specifies origin of rotation. The result is similar as performing
%   translation(-X0, -Y0, -Z0), rotation, and translation(X0, Y0, Z0).
%
%
%   See also
%   transforms3d, transformPoint3d, createRotationOx, createRotationOy
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-04-06
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

% default values
dx = 0;
dy = 0;
dz = 0;
theta = 0;

% get input values
if length(varargin) == 1
    % only one argument -> rotation angle
    theta = varargin{1};

elseif length(varargin) == 2
    % origin point (as array) and angle
    var = varargin{1};
    dx = var(1);
    dy = var(2);
    dz = var(3);
    theta = varargin{2};

elseif length(varargin) == 4
    % origin (x and y) and angle
    dx = varargin{1};
    dy = varargin{2};
    dz = varargin{3};
    theta = varargin{4};
end

% compute coefs
cot = cos(theta);
sit = sin(theta);

% create transformation
trans = [...
    cot -sit 0 0;...
    sit  cot 0 0;...
    0 0 1 0;...
    0 0 0 1];

% add the translation part
t = [1 0 0 dx;0 1 0 dy;0 0 1 dz;0 0 0 1];
trans = t * trans / t;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%ZNCC  Normalized cross correlation
%
% M = ZNCC(I1, I2) is the zero-mean normalized cross-correlation between the
% two equally sized image patches I1 and I2.  The result M is a scalar in
% the interval -1 to 1 that indicates similarity.  A value of 1 indicates
% identical pixel patterns.
%
% Notes::
% - The ZNCC similarity measure is invariant to affine changes in image
%   intensity (brightness offset and scale).
%
% See also NCC, SAD, SSD, ISIMILARITY.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB 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.
%
% MVTB 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 MVTB.  If not, see <http://www.gnu.org/licenses/>.

function m = zncc(w1, w2)

	w1 = w1 - mean(w1(:));
	w2 = w2 - mean(w2(:));

	denom = sqrt( sum(sum(w1.^2))*sum(sum(w2.^2)) );

	if denom < 1e-10,
		m = 0;
	else
		m = sum(sum((w1.*w2))) / denom;
	end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [theta] = trainLinearReg(X, y, lambda)
%TRAINLINEARREG Trains linear regression given a dataset (X, y) and a
%regularization parameter lambda
%   [theta] = TRAINLINEARREG (X, y, lambda) trains linear regression using
%   the dataset (X, y) and regularization parameter lambda. Returns the
%   trained parameters theta.
%

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

% Create "short hand" for the cost function to be minimized
costFunction = @(t) linearRegCostFunction(X, y, t, lambda);

% Now, costFunction is a function that takes in only one argument
options = optimset('MaxIter', 200, 'GradObj', 'on');

% Minimize using fmincg
% XXX(SaveTheRbtz@): Disable warnings here?
theta = fmincg(costFunction, initial_theta, options);

end

author: SaveTheRbtz

repo: ml-class


text:
% pop_readneurone() - A function for reading NeurOne data.
%
% Usage: >> [EEGOUT,command] = pop_readneurone() % Launches GUI for user input
%       OR
%        >> [EEGOUT,command] = pop_readneurone(dataPath, sessionPhaseNumber, chans)
%
% Please note that if this function is called without using the menu item
% in EEGLAB, it will not be automatically stored in its memory.
% ========================================================================
% Inputs:
%        dataPath           - A string containing the full path for .ses
%                             -file to be read. This input argument is
%                             required.
%        sessionPhaseNumber - A number indicating the session
%                             phase to be observed. This input argument is
%                             optional. If the session phase number is left
%                             empty, default value 1 will be used.
%        chans              - Another optional input argument to indicate
%                             the channels to be imported. Acceptable
%                             values are strings where the channels are
%                             in a vector form without brackets e.g.
%                             '1,3,5' or '1:4' (when the GUI is executed,
%                             the quotation marks should not be used).
%                             Note that only acceptable values
%                             are those between 1 and the total number of
%                             channels. It should also be noted that the
%                             number of the channel may not be the same as
%                             its input number since when an input
%                             number is missing, the channel numbering will
%                             continue with next available input number.
%                             Here is an example of the channel numbering
%                             (note the missing input numbers 3 and 4):
%
%              Input number               Channel number
%             ------------               --------------
%                   1          ----->           1
%                   2          ----->           2
%                   5          ----->           3
%                   6          ----->           4
%
%
% Outputs:
%        EEGOUT             - A struct with the same fields as a standard
%                             EEG data structure under EEGLAB
%                             (see eeg_checkset.m).
%        command            - A string containing all the input parameters
%                             and the command used to execute this
%                             function.
% ========================================================================
% GUI help:
%
% When this function is used without any input arguments, it will open
% several windows to ask for required information as described above. The
% first pop-up window asks the user to specify the .ses -file to be loaded.
% If 'Cancel' -button is pressed at this point, the importing process will
% terminate. The next window is a GUI asking for information about the
% channels to be loaded and the number of the desired session phase to be
% observed. When 'Ok' -button is pressed, the data is sent to readneurone.m
% for further processing. Again, if 'Cancel' -button is pressed, the whole
% importing process will be terminated.
%
% IMPORTANT NOTE: For a computer with high memory capacity it is much more
%                 faster to import data by leaving the channel field empty
%                 and removing unnecessary channels afterwards. However,
%                 in case of limited memory it is recommended to first
%                 specify the channels to be imported.
%
% ========================================================================
% NOTE:
% This file is part of the NeurOne data import plugin for EEGLAB.
% ========================================================================
%
% Current version: 1.0.3.4 (2016-06-17)
% Author: Mega Electronics
%
% see also: eegplugin_neurone.m or readneurone.m

function [EEG,command] = pop_readneurone(dataPath,sessionPhaseNumber,chans)
%% Initialize empty structure for the data

EEG = {};              % Create empty space for data
command = '';

%% Check the number of input arguments
nargchk(0,1,nargin);

%% Use the GUI

if nargin==0

    if isfield(evalin('base','EEG'),'history')
        history = evalin('base','EEG.history');
        if  ~(isempty(strfind(history,'dataPath')))
            dataPathIndex = strfind(history,'dataPath=''');
            idx = strfind(history,'''');
            idx = idx(idx>dataPathIndex);
            dataPath = history(idx(1)+1:idx(2)-1);
            dataPath = fileparts(dataPath);
            [dataPath folder] = fileparts(dataPath);
        else
            dataPath = cd;
        end
    else
        dataPath = cd;
    end

    [fname dataPath] = uigetfile({'*.ses', ...
        'NeurOne session file (*.ses)'}, ...
        'Load NeurOne session file',dataPath);

    % If cancel is pressed, the import process will terminate
    if fname==0
        return
    end

    % If ok, manipulate the filename
    fullpath = [dataPath fname];
    a = dir(dataPath);
    for k = 3:numel(a)
        if a(k).isdir==1
            if ~(isempty(regexpi(fullpath,a(k).name)))
                dataPath = [dataPath a(k).name filesep];
            end
        end
    end

    fprintf('\n----------IMPORTING DATA----------\n')
    fprintf('File to Import: %s\n',fullpath)

    guifig=guireadneurone;

    try
       SETTINGS=guidata(guifig);
    catch
       return;
    end

    switch SETTINGS.loadStatus
    % When cancel button is pressed, the GUI terminates
    case 0
        disp('NeurOne import cancelled')
        delete(guifig);
        return
    % When ok button is pressed, proceed to read the data
    case 1
        delete(guifig);

        SETTINGS.dataPath=dataPath;

        if  isempty(SETTINGS.chans)
            fprintf('Channels to be read: all\n');
        else
            fprintf('Channels to be read: %s\n',SETTINGS.chans);
        end

        EEG = readneurone(SETTINGS.dataPath, ...
            SETTINGS.sessionPhaseNumber,SETTINGS.chans);
        EEG = eeg_checkset(EEG); % EEGLAB function to check consistency
        EEG = eeg_checkset(EEG,'eventconsistency'); % check event structure
        EEG = eeg_checkset(EEG,'makeur');
    end

% ====================================================================
% Read data using direct command

else
    if nargin>=1
         fullpath=dataPath;
         fprintf('\n----------IMPORTING DATA----------\n')
         fprintf('File to Import: %s\n',fullpath)

      % Manipulate the given path
      [dataPath fname] = fileparts(fullpath);
      a = dir(dataPath);
      for k = 3:numel(a)
          if a(k).isdir==1
             if ~(isempty(regexpi(fullpath,a(k).name)))
                dataPath = [dataPath filesep a(k).name filesep];
             end
          end
      end
         SETTINGS=struct();
         SETTINGS.dataPath=dataPath;

         if ~(exist('sessionPhaseNumber'))
            SETTINGS.sessionPhaseNumber = 1;
         else
            SETTINGS.sessionPhaseNumber=sessionPhaseNumber;
         end

         if ~(exist('chans'))
            SETTINGS.chans = '';
         else
            SETTINGS.chans=chans;
         end

    end

    EEG = readneurone(SETTINGS.dataPath,SETTINGS.sessionPhaseNumber, ...
        SETTINGS.chans);
    EEG = eeg_checkset(EEG);                    % check EEG data structure
    EEG = eeg_checkset(EEG,'eventconsistency'); % check event structure
    EEG = eeg_checkset(EEG,'makeur');

end

comline = sprintf('dataPath=''%s'';',SETTINGS.dataPath);
comline = sprintf('%s sessionPhaseNumber=''%d'';',...
    comline,SETTINGS.sessionPhaseNumber);
comline = sprintf('%s chans=''%s'';',...
    comline,SETTINGS.chans);
command = sprintf('%s [EEG,COM]=pop_readneurone(dataPath,sessionPhaseNumber,chans);',comline);

return

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [u, flag] = minL1Tikhonov(f, gamma, A, opts)
%minL1Tikhonov Tikhonov regularization with L^1 data term
%
% Solves the optimization problem
%  \gamma \| u \|_2^2 + \| A u - f\|_1 -> min
% by a semi-smooth Newton method
%
% u = minL1Tikhonov(f, gamma, A)
%
%
% Reference:
%  Clason, Jin, Kunisch
%  A semismooth Newton method for L^1 data fitting with automatic choice of regularization parameters and noise calibration
%  SIAM Journal on Scientific Computing, 2010
%

% written by M. Storath
% $Date: 2012/10/29 01:19:08 $	$Revision: 0.1 $

%%-------------------------------------------------------------------------
% init
m = numel(f);
maxit = [];
init = [];
flag = 0;
if exist('opts', 'var')
    if isfield(opts, 'mat')
        AAT = opts.mat;
    end
    if isfield(opts, 'maxit')
        maxit = opts.maxit;
    end
    if isfield(opts, 'init')
        init = opts.init;
    end
    if isfield(opts, 'L')
        L = opts.L;
    end
end
if ~exist('AAT', 'var')
    AAT = A * A';
end
if ~isa(A, 'linop') && ~exist('L', 'var')
    L = spconvmatrix([-1, 2, -1], m);
end

% diagonal idx
diagidx = eye(m) == 1;

% standard parameters (see Clason, Jin, Kunisch 2010)
iMax = 20; %(iMax higher than standard)
c = 1e9;
beta = 1;
q = 0.2;

% init variables
p = zeros(m,1);
pos = p;
neg = p;
alpha = 2 * gamma;

%deactivate warnings temporaily
warning off;

%%-------------------------------------------------------------------------
% main program
while beta > 1e-16
    pOld = p;

    % precompute matrices
    if ~isa(A, 'linop')
        M = (1/alpha) * AAT + beta * L;
        diagM = M(diagidx);
    end

    % newton step
    for k = 1:iMax
        % compute active sets
        posNew = p > +1;
        negNew = p < -1;
        comb = posNew | negNew;

        % solve linear equation

        b = f(:) + c * (posNew - negNew);
        if isa(A, 'linop')
            M = @(x) reshape(AAT * (x / alpha), m, 1) + beta * conv(x(:), [-1; 2; -1], 'same') + c * comb .* x;
            [p, flag] = cgs(M, b, [], maxit, [], [], init);
        else
            M(diagidx) = diagM + c * comb;
            p = M \ b;
        end

        % if active sets did not change, break
        if ~(any(posNew - pos) || any(negNew - neg))
            break;
        end

        % set new active sets
        pos = posNew;
        neg = negNew;
    end

    if (k == iMax) && any(abs(p)>2)
        p = pOld;
        break;
    end

    % decrease beta
    beta = beta * q;
end

% restore primal variable
u = A' * (p / alpha);

% reactivate warnings
warning on;

author: mstorath

repo: Pottslab

save_path: github-repos/MATLAB/mstorath-Pottslab


text:
function varargout = solver_sDantzig( A, b, delta, mu, x0, z0, opts, varargin )
% SOLVER_SDANTZIG Dantzig selector problem. Uses smoothing.
%[ x, out, opts ] = solver_sDantzig( A, b, delta, mu, x0, z0, opts )
%    Solves the smoothed Dantzig
%        minimize norm(x,1) + (1/2)*mu*norm(x-x0).^2
%        s.t.     norm(D.*(A'*(A*x-b)),Inf) <= delta
%    by constructing and solving the composite dual
%        maximize - g_sm(z) - delta*norm(z,1)
%    where
%        gsm(z) = sup_x <z,D.*A'*(Ax-b)>-norm(x,1)-(1/2)*mu*norm(x-x0)
%    A must be a linear operator, b must be a vector, and delta and mu
%    must be positive scalars. Initial points x0 and z0 are optional.
%    The standard calling sequence assumes that D=I. To supply a scaling,
%    pass the cell array { A, D } instead of A. D must either be a scalar,
%    a vector of weights, or a linear operator.

% Supply default values
error(nargchk(4,8,nargin));
if nargin < 5, x0   = []; end
if nargin < 6, z0   = []; end
if nargin < 7, opts = []; end

% -- legacy options from original software --
if isfield(opts,'lambda0')
    z0 = opts.lambda0;
    opts = rmfield(opts,'lambda0');
end
if isfield(opts,'xPlug')
    x0 = opts.xPlug;
    opts = rmfield(opts,'xPlug');
end
if isfield(opts,'solver')
    svr     = opts.solver;
    opts    = rmfield(opts,'solver');
    if isfield(opts,'alg') && ~isempty(opts.alg)
        disp('Warning: conflictiong options for the algorithm');
    else
        % if specified as "solver_AT", truncate:
        s = strfind( svr, '_' );
        if ~isempty(s), svr = svr(s+1:end); end
        opts.alg = svr;
    end
end

% Extract the linear operators
D = [];
if isa( A, 'cell' ),
    if length(A) > 1, D = A{2}; end
    A = A{1};
end
if isempty(D),
    D = @(x)x;
elseif isa( D, 'double' ),
    D = @(x)D.*x;
end
if isa( A, 'double' ),
    A = linop_matrix(A);
end

% Call TFOCS
objectiveF = prox_l1;
affineF    = { @(y,mode)linear_DS( D, A, y, mode ), -D(A(b,2)) };
dualproxF  = prox_l1( delta );
[varargout{1:max(nargout,1)}] = ...
    tfocs_SCD( objectiveF, affineF, dualproxF, mu, x0, z0, opts, varargin{:} );

% Implements x -> D*A'*A*x and its adjoint if A is a linop
function y = linear_DS( D, A, y, mode )
switch mode,
case 0,
    y = A([],0);
    if iscell( y ),
        y = { y{1}, y{1} };
    else
        y = { [y(2),1], [y(2),1] };
    end
case 1, y = D(A(A(y,1),2));
case 2, y = A(A(D(y),1),2);
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function wfdbdemo()
% WFDB App Toolbox Demo
%
% Written by Ikaro Silva, 2013
% Last modified: January 10, 2014
%

[~,config]=wfdbloadlib;
echo on;
display('Reading samples ECG signal from MIT-BIH Arrhythmia Database');
N=10000;
[ecg,Fs,tm]=rdsamp('mitdb/100',1,N);

display('Reading and plotting annotations (human labels) of QRS complexes performend on the signals');
%by cardiologists.
[ann,type,subtype,chan,num]=rdann('mitdb/100','atr',1,N);

%Plot 2D version of signal and labels
figure;
plot(tm(1:N),ecg(1:N));hold on;grid on;
plot(tm(ann(ann<N)+1),ecg(ann(ann<N)+1),'ro');

%Stack the ECG signals based on the labeled QRSs
display('Ploting 3D version of signal and labels');
[RR,tms]=ann2rr('mitdb/100','atr',N);
delay=round(0.1/tm(2));
M=length(RR);
offset=0.3;
stack=zeros(M,min(RR))+NaN;
qrs=zeros(M,2)+NaN;
for m=1:M
    stack(m,1:min(RR)+1)=ecg(tms(m)-delay:tms(m)+min(RR)-delay);
    qrs(m,:)=[delay+1 ecg(tms(m))];
end

%Generate plot inspired by
%Clifford GD, Azuaje F, McSharry PE, editors.
%Advanced Methods and Tools for ECG Analysis.
%1st ed., Norwood, MA, USA: Artech House; 2006. (Engineering in Medicine and Biology; 1).
figure;
[X,Y] = meshgrid(1:min(RR)+1,1:M);
surf(Y,X,stack);hold on;grid on;
shading interp;
plot3(1:M,qrs(:,1),qrs(:,2)+offset,'go-','MarkerFaceColor','g');
view(120, 30);
axis off;

%Generate plot inspired by
% Samenie et al
% "Filtering Noisy ECG Signals Using Extended Kalman Filter Based on a
% Modified Dynamic ECG Model"
% Computers in Cardiology, 2005
figure;
stack=stack';
stack=stack(:)+[0:length(stack(:))-1]'.*0.0005;
theta=linspace(0,M*2*pi,length(stack));
x=sin(theta);y=cos(theta);
plot3(x,y,stack,'b');
grid on;
axis off;

%Display information about databases availabe in PhysioNet
fprintf('**Querying PhysioNet for available databases...\n');
db_list=physionetdb;
db_size=length(db_list);
fprintf(['\tYou currently have access to ' num2str(db_size) ' databases for download in PhysioNet (type ''help physionetdb'' for more info)!\n']);

display('Demoing finished !!');
display('For more information about the toolbox, type ''wfdb'' at the command prompt.');
echo off;

author: ikarosilva

repo: wfdb-app-toolbox


text:
function showresult3(node,elem,u,expr,varargin)
%% SHOWRESULT3 display the mesh and the solution in 3-D
%
%  showresult3(node,elem,u,viewangle) displays the mesh and the solution in
%  one figure. The left one is the mesh, the middle
%  one is the contour of the solution, and the right one is the graph of
%  the function. The last viewangle is used to adjust the view angle of the
%  graph of the function.
%
%  Example:
%
% See also showrate, showmesh, showsolution
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('expr','var'), expr = []; end
set(gcf,'Units','normal');
set(gcf,'Position',[0.25,0.25,0.6,0.4]);
% show mesh
if size(elem,1) < 1e6
    subplot(1,2,1);
    showboundary3(node,elem,expr,varargin{:});
else
    subplot(1,2,1);
    title('The mesh is too dense to display')
end
% show solution
subplot(1,2,2);
showsolution3(node,elem,u,expr,varargin{:});
colorbar;
pause(0.05)

repo: ifem


text:
function [ind]=ij2ind(a,i,j)
[m,n]=size(a);
ind=m*i-j+1;

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function asa266_test ( )

%*****************************************************************************80
%
%% ASA266_TEST tests the ASA266 library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 August 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA266_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Tests for the ASA266 package.\n' );

  asa266_test01 ( );
  asa266_test02 ( );
  asa266_test03 ( );
  asa266_test04 ( );
  asa266_test05 ( );
  asa266_test06 ( );
  asa266_test07 ( );
  asa266_test08 ( );
  asa266_test09 ( );
  asa266_test10 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA266_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function v = axis(q,varargin)
% rotational axis of the quaternion
%
% Syntax
%   v = axis(q)
%
% Input
%  q - @quaternion
%
% Output
%  v - @vector3d

v = vector3d(q.b,q.c,q.d);
v(q.a<0) = -v(q.a<0);
v(isnull(norm(v))) = xvector;
v = v./norm(v);

author: mtex-toolbox

repo: mtex

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


text:
report_this_filefun(mfilename('fullpath'));

% compute the source zones for landers

dy = 0.15;
tmin1 = 3;
newt0 = newt2;
tpre =  120;

l=storedcat.Date <= min(maepi.Date) + days(60);
a=storedcat.subset(l);
update(mainmap()); newt2 = a;
timeplot

d = [];
prol = 0.90;

l = a.Magnitude > 7;
maepi = a.subset(l);

mati = maepi(1,3);
M = 7.4 - 4;
da = []; anz = [];
B2 = [];
t0 = (max(a.Date) - mati)*365;
tlen = t0;
cd /home2/stefan/ZMAP/aspar
fid2 = fopen('sourczones.txt','w');

dx = 0.33;

for x = 29:dx:31
    l = a.Longitude >= x & a.Longitude < x+dx;
    b1 = prctile2(a(l,2),20);
    b2 = prctile2(a(l,2),80);

    newt2 = a.subset(l);
    mcperc_ca3;
    if isnan(Mc95) == 0
        magco = Mc95;
    elseif isnan(Mc90) == 0
        magco = Mc90;
    else
        [bv magco stan av me mer me2,  pr] =  bvalca3(b,1,1);
    end
    magco
    l = newt2.Magnitude >= magco;
    newt2 = newt2(l,:);

    figure_w_normalized_uicontrolunits(map)
    try
        delete(pl)
    catch ME
        disp(ME.message);
    end

    hold on
    pl = plot(newt2.Longitude,newt2.Latitude,'xk');drawnow
    timeplot

    calcp
    %da = [da ; x+dx/2 P];
    B2 = [B2 ; b1 x  b2 x];

    anz = [];
    for m2 = 4.75:0.5:7.25
        M2 = maepi(1,6) - m2;
        t0 = tlen;
        pla = 0; pla2 = 0;
        dt = 0.5;
        for t = t0:dt:t0+tpre
            pla = pla + (10^(A + bv*(M2)) * (t + c)^(-p))  *dt;
            pla2 = pla2 + (10^(A + bv*(M2-0.5)) * (t + c)^(-p))  *dt;
        end

        anz = [anz ;  m2+0.25  (pla-pla2)/tpre];
    end
    % anz(5,2) = anz(5,2)/10;

    % write info to file
    s = ['0    1.     -1          zn03']; s = s';
    fprintf(fid2,'%s\n',s);
    s = ['2 1 1']; s = s';
    fprintf(fid2,'%s\n',s);
    s = [num2str(x) ' ' num2str(b1,5) ' ' num2str(x) ' ' num2str(b2,5)]; s = s';
    fprintf(fid2,'%s\n',s);
    s = [num2str(x+dx) ' ' num2str(b1,5) ' ' num2str(x+dx) ' ' num2str(b2,5)]; s = s';
    fprintf(fid2,'%s\n',s);
    fprintf(fid2,'%7.6f    ',anz(:,2));
    fprintf(fid2,'\n');
    fprintf(fid2,'\n');
    fprintf(fid2,'%3.2f  ',anz(:,1));
    fprintf(fid2,'\n');
end

%end

fclose(fid2)

do = [' ! cat head_turk.txt | sed -e"s/sub1/' num2str(prol) ' 1 ' num2str(tpre) '/" > head2.txt ' ]; eval(do)
do = [' ! cat head2.txt sourczones.txt tail.txt > /home2/stefan/srisk/myrisk.inp']; eval(do)

cd /home2/stefan/srisk/

do = [ '! /nfs/alaska/home2/stefan/srisk/seis4b.exe  myrisk.inp myrisk.out f2 f3' ]; eval(do)
do = [' !cat myrisk.out | grep -e "LAT   "  -e "' num2str(tpre) ' YE"  > tmp2 ']; eval(do)
do = ['condata2']; err = [' ']; eval(do,err);
save inpudata.xyz da -ascii

eq = [a.Longitude a.Latitude a.Magnitude];
save eqs2.dat eq -ascii

cd /home2/stefan/srisk/

do = [' ! /home2/stefan/srisk/myriskturk 0.05 0.02 ' num2str(tpre) ]; eval(do)

load lat
load lon
load hpga

[X,Y] = meshgrid(min(lon)-0.1:0.08:max(lon)+0.1,min(lat)-0.1:0.08:max(lat)+0.1);

Z = griddata(lon,lat,hpga,X,Y,'linear');
figure

pcolor(X,Y,Z); colorbar; shading flat
hold on
shading interp
overlay

author: CelsoReyes


text:
function pkdo = umSYMBPKDO2d(i,j)

  r = sym('r');
  s = sym('s');

  a = 2*(1+r)/(1-s) - 1;
  b = s;

  pkdo = umSYMBJACOBI1d(a,i,0,0);
  pkdo = pkdo.*umSYMBJACOBI1d(b,j,2*i+1,0).*((0.5*(1-b)).^i);

  % normalize (sets L2 norm of each basis function to be 1)
  pkdo = pkdo*sqrt((i+0.5)*(i+j+1));

  % simplify
  pkdo = simple(pkdo);

repo: Aero-matlab


text:
function [Rt]=RateSimCIR(theta,kappa,sigma,lambda,dt,ratestart,months,tau)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CIR term structure simulation, inputs below
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clear all
% theta=0.10;
% kappa=0.05;
% sigma=0.075;
% lambda=-0.4;
% dt=1/12;
% months=120;
% ratestart=0.10;
% tau=[3/12,6/12,2,5];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Short Rate Dynamics
rt(1)=ratestart;
for i=1:months*4
    rt(i+1)=rt(i)+kappa*(theta-rt(i))*dt/4+sqrt(rt(i))*sigma*sqrt(dt/4)*randn(1);
end
% Term Structure Dynamics
for i=1:months
    rttemp=rt(i*4-3);
    rttemp1(i)=rt(i*4-3);
    for j=1:numel(tau)
        AffineG=sqrt((kappa+lambda)^2+2*sigma^2);
        AffineB=2*(exp(AffineG*tau(j))-1)/((AffineG+kappa+lambda)...
            *(exp(AffineG*tau(j))-1)+2*AffineG);
        AffineA=2*kappa*theta/(sigma^2)*log(2*AffineG*...
            exp((AffineG+kappa+lambda)*tau(j)/2)/((AffineG+kappa+lambda)*...
            (exp(AffineG*tau(j))-1)+2*AffineG));
        A(j)=-AffineA/tau(j);
        B(j)=AffineB/tau(j);
        Rt(i,j)=A(j)+B(j)*rttemp;
    end
end
Rt;

%%%%%%%%%%%%% MAKE THE GRAPH %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rt;
% figure(2)
% surf(tau,1:120,Rt)
% xlabel('Time to Maturity (Fixed)')
% ylabel('Months Passed')
% title('A Given Simulation of the Term Structure (CIR)')

author: Sable

repo: mcbench-benchmarks


text:
function fx = p11_fun ( n, x )

%*****************************************************************************80
%
%% P11_FUN evaluates the integrand for problem 11.
%
%  Discussion:
%
%    S&S gives exact value as pi =  3.1415926535897932385...
%    S&S gives Laguerre(16) as      2.6652685196...
%    S&S gives EXP_TRANSFORM(16) as 2.3629036166...
%
%  Integral:
%
%    Integral ( 0 <= x < +oo ) 1/((1+x)*sqrt(x)) dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the function values.
%
  for i = 1 : n
    if ( x(i) == 0.0 )
      fx(i) = 0.0;
    else
      fx(i) = 1.0 / ( ( 1.0 + x(i) ) * sqrt ( x(i) ) );
    end
  end

  return
end


text:
function fileNameOut = lrsRun(modelName, param)
% Runs lrs in the command line either using a shell script or launched
% directly from within Matlab.
%
% USAGE:
%
%     fileNameOut = runLrs(modelName, positivity, inequality, shellScript)
%
% INPUTS:
%    modelName:      name of the file to give as input to lrs
%
%    param:         parameter structure with the following fields:
%    *.positivity:    if equals to 1, then positive orthant base
%    *.inequality:    if equals to 1, then represent as two inequalities rather than a single equality
%    *.shellScript:   if equals to 1, then lrs is run through a bash script
%    *.redund         if equals to 0, then remove redundant linear equalities
%
% OUTPTS
% fileNameOut       full path to output file

% Ronan Fleming 2021

if ~exist('param','var')
    param = struct();
end
if ~isfield(param,'positivity')
    param.positivity  = 0;
end
if ~isfield(param,'inequality')
    param.inequality  = 0;
end
if ~isfield(param,'shellScript')
    param.shellScript  = 0;
end
if ~isfield(param,'facetEnumeration')
    %assume vertex enumeration, unless specified that it is facet enumeration
    param.facetEnumeration  = 1;
end
if ~isfield(param,'redund')
    param.redund  = 1;
end

if contains(modelName,'.ine')
    param.facetEnumeration  = 0;
end

if contains(modelName,'.ine') || contains(modelName,'.ext')
    modelName = strrep(modelName,'.ine','');
    modelName = strrep(modelName,'.ext','');
else
    if param.inequality == 0
%         if param.positivity == 0
%             modelName = [modelName '_pos_eq'];
%         else
%             modelName = [modelName '_neg_eq'];
%         end
    else
        if param.positivity == 1
            modelName = [modelName '_pos_ineq'];
        else
            modelName = [modelName '_neg_ineq'];
        end
    end
    if ~param.redund
        modelName = [modelName '_noR'];
    end
end

if param.facetEnumeration
    [status, result] = lrsSpecifyRows([modelName '.ext']);
else
    [status, result] = lrsSpecifyRows([modelName '.ine']);
end

if param.redund==0
    %remove all redundant halfspaces
    %http://cgm.cs.mcgill.ca/%7Eavis/C/lrslib/USERGUIDE.html#redopt
    redundCmd = ['redund ' int2str(0) ' ' int2str(0)];
end

if isunix
    [status, result] = system('which lrs');
    if ~isempty(result)
        if param.shellScript
            % call lrs through a bash script and wait until extreme pathways have been calculated
            systemCallText = ['sh ' pwd filesep modelName '.sh'];
            [status, result] = system(systemCallText);
            if status == 0
                error(['Failure to run Bash script ' pwd filesep modelName '.sh']);
            end
        else
            if param.facetEnumeration
                fileNameOut = [modelName '.ine'];
                % call lrs to compute halfspaces
                systemCallText = ['lrs ' pwd filesep modelName '.ext > ' pwd filesep fileNameOut];
                localCallText = ['lrs ' modelName '.ext > ' fileNameOut];
                disp(localCallText)
            else
                fileNameOut = [modelName '.ext'];
                % call lrs to compute vertices
                systemCallText = ['lrs ' pwd filesep modelName '.ine > ' pwd filesep fileNameOut];
                localCallText = ['lrs ' modelName '.ine > ' fileNameOut];
                disp(localCallText)
            end
            [status, result] = system(systemCallText);
            if status == 1
                error(['lsr failed on file ', pwd filesep modelName '.ine']);
            end
            if param.facetEnumeration
                if param.redund == 0
                    fileNameIn = fileNameOut;
                    fileNameOut = [modelName '_noR.ine'];
                    uselrsRedundOption =0;

                    [status, result] = lrsSpecifyRows(fileNameIn);

                    if uselrsRedundOption
                        [status, result] =  system(['sed -i ''s/end/end\nredund 0 0/g'' ' fileNameIn]);
                        systemCallText = ['lrs ' pwd filesep fileNameIn ' > ' pwd filesep fileNameOut];
                        localCallText = ['lrs ' fileNameIn ' > ' fileNameOut];
                    else
                        systemCallText = ['redund ' pwd filesep fileNameIn ' > ' pwd filesep fileNameOut];
                        localCallText = ['redund ' fileNameIn ' > ' fileNameOut];
                    end
                    disp(localCallText);
                    [status, result] = system(systemCallText);
                    if status == 1
                        disp(result)
                        error(['lsr failed on file ', pwd filesep modelName '.ine']);
                    end
                end
            end
        end
    else
        fprintf('lrs not installed or not in path\n');
    end
else
    fprintf('non unix machines not yet supported\n');
end
end

author: opencobra

repo: cobratoolbox


text:
function varargout = erfc(varargin)
%ERFC (overloaded)

switch class(varargin{1})

    case 'double'
        error('Overloaded SDPVAR/ERFC CALLED WITH DOUBLE. Report error')

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = struct('convexity','none','monotonicity','decreasing','definiteness','positive','model','callback');
        operator.bounds = @bounds;
        operator.range = [0 2];
        operator.derivative =@(x)-exp(-x.^2)*2/sqrt(pi);
        operator.inverse = @(x)(erfcinv(x));

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error('SDPVAR/ERF called with CHAR argument?');
end

function [L,U] = bounds(xL,xU)
L = erfc(xU);
U = erfc(xL);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function pass = test_eig( pref )
% Test for EIG command of Chebfun2.

if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 1e4*pref.cheb2Prefs.chebfun2eps;

% Decomposition on a square domain.
f = cheb.gallery2('challenge');
[V, D] = eig(f);

% Is it correct?
pass(1) = norm(f * V - V * D) < tol;

%%
% The following should give an error message as the domain is not square.
f = chebfun2(@(x,y,z) x+y, [-1 1 -2 2]);
try
    d = eig(f);
    pass(2) = false;
catch ME
    pass(2) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN2:eig:domainerr');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function out = flattenMaskOverlay(img, mask, varargin)
% FLATTENMASKOVERLAY  Create a flattened image with overlayed mask (with transparency)
%  FLATTENMASKOVERLAY(IMG, MASK) applies a MASK (defined as logical) on top
%     of the original IMG and displays the result in the current axes. The
%     mask (locations of 1's) will be displayed as a red color with 50%
%     opacity. The resulting image will be a flat image (no use of
%     AlphaData), and therefore will have better interactive performance
%     (zooming, panning, data cursors). IMG must be MxN (gray scale) or
%     MxNx3 (RGB). Supported data types are double, uint8, uint16. MASK
%     must be an MxN logical array.
%
%  FLATTENMASKOVERLAY(IMG, MASK, ALPHALEVEL, OVERLAYCOLOR) specifies the
%     opacity level and the overlay color of the mask. ALPHALEVEL must be a
%     scalar between 0 to 1. Pass in [] for default (0.5). OVERLAYCOLOR can
%     be a string indicating the color (see list below) or a 1x3 double
%     array (0 to 1) representing the RGB values. Default is red.
%            b     blue
%            g     green
%            r     red
%            c     cyan
%            m     magenta
%            y     yellow
%            k     black
%            w     white
%
%  IM = FLATTENMASKOVERLAY(IMG, MASK, ...) returns the image data without
%     displaying the image.
%
%  Example:
%     load mandrill
%     rgb = ind2rgb(X, map);
%
%     subplot(2, 2, 1);
%     image(rgb); axis image off; title('Original image');
%
%     subplot(2, 2, 2);
%     % try to segment out the nose
%     mask = rgb(:, :, 1) > 0.8 & rgb(:, :, 2) < 0.3 & rgb(:, :, 3) < 0.3;
%     image(mask); colormap(gray(2)); axis image off; title('Mask image');
%
%     subplot(2, 2, 3);
%     flattenMaskOverlay(rgb, mask, [], 'g');
%     title('Green mask with 50% opacity');
%
%     subplot(2, 2, 4);
%     flattenMaskOverlay(rgb, mask, 1, 'g');
%     title('Green mask with 100% opacity');

% Jiro Doke
% v1.0 - Nov 2011
% v1.1 - Nov 2011. Bug in default alphaLevel and typo in example.
%
% Copyright 2011 The MathWorks, Inc.

%%

% Default alpha level
DEFAULT_ALPHA = 0.5;

% Default overlay color
DEFAULT_COLOR = [1 0 0];

% Image size
sz = size(img);

%% Input Validation
p = inputParser();
p.addRequired('img', @isvalidImageData);
p.addRequired('mask', @(x) isvalidMaskData(x, sz));
p.addOptional('alphaLevel', DEFAULT_ALPHA, @(x) ...
    assert(isempty(x) || ...
    (isa(x, 'double') && numel(x) == 1 && x >= 0 && x <= 1), ...
    sprintf('%s:InvalidAlphaLevel', mfilename), ...
        '''alphaLevel'' must be a scalar between [0, 1]'));
p.addOptional('overlayColor', DEFAULT_COLOR, @(x) ...
    validateattributes(x, {'char', 'double'}, {'2d'}));

p.parse(img, mask, varargin{:});

img = p.Results.img;
mask = p.Results.mask;
alphaLevel = p.Results.alphaLevel;
overlayColor = p.Results.overlayColor;

if isempty(alphaLevel)
    alphaLevel = DEFAULT_ALPHA;
end

if ischar(overlayColor)
    switch lower(overlayColor)
        case {'y','yellow'}
            overlayColor = [1 1 0];
        case {'m','magenta'}
            overlayColor = [1 0 1];
        case {'c','cyan'}
            overlayColor = [0 1 1];
        case {'r','red'}
            overlayColor = [1 0 0];
        case {'g','green'}
            overlayColor = [0 1 0];
        case {'b','blue'}
            overlayColor = [0 0 1];
        case {'w','white'}
            overlayColor = [1 1 1];
        case {'k','black'}
            overlayColor = [0 0 0];
        case {''}
            overlayColor = DEFAULT_COLOR;
        otherwise
            error(sprintf('%s:UnknownOverlayColor', mfilename), ...
                'Unknown overlayColor "%s".', overlayColor);
    end
else
    if isempty(overlayColor)
        overlayColor = DEFAULT_COLOR;
    end
    assert(isequal(size(overlayColor), [1 3]) && ...
        nnz(overlayColor < 0 | overlayColor > 1) == 0, ...
        sprintf('%s:InvalidOverlayColor', mfilename), ...
        '''overlayColor'' must be a 1x3 array between [0, 1]');
end

%% Create overlay
rMask = mask * overlayColor(1);
gMask = mask * overlayColor(2);
bMask = mask * overlayColor(3);

[img_r, img_g, img_b] = getDoubleRGB(img);

% Red plane
tmp = img_r + (rMask - img_r) * alphaLevel;
img_r(mask) = tmp(mask);

% Blue plane
tmp = img_g + (gMask - img_g) * alphaLevel;
img_g(mask) = tmp(mask);

% Green plane
tmp = img_b + (bMask - img_b) * alphaLevel;
img_b(mask) = tmp(mask);

% Convert back to the original data type
if isa(img, 'double');
    finalImg = cat(3, img_r, img_g, img_b);
else
    finalImg = cast(cat(3, img_r, img_g, img_b) * ...
        double(intmax(class(img))), class(img));
end

if nargout == 0
    image(finalImg); axis image off
else
    out = finalImg;
end

end

%% Validate IMG
function tf = isvalidImageData(img)

tf = true;
validateattributes(img, {'double', 'uint8', 'uint16'}, {'nonempty'}, ...
    mfilename, 'img');

if isa(img, 'double')
    % Check that the values are between 0 and 1
    assert(nnz(img < 0 | img > 1) == 0, ...
        sprintf('%s:InvalidDoubleImageRange', mfilename), ...
        'Values in a double ''img'' must be between [0, 1].');

end

% IMG should be a 2D or an MxNx3 array
sz = size(img);
assert((numel(sz) == 2) || (numel(sz) == 3 && sz(3) == 3), ...
    sprintf('%s:InvalidImageSize', mfilename), ...
    'Size of ''img'' must be MxN or MxNx3.');

end

%% Validate MASK
function tf = isvalidMaskData(mask, sz)

tf = true;
thisSZ = size(mask);
assert(isequal(thisSZ, sz(1:2)), ...
    sprintf('%s:InvalidMaskSize', mfilename), ...
    'Size of ''mask'' must be %dx%d.', sz(1), sz(2));

end

%% Convert IMG to MxNx3 double array
function [r, g, b] = getDoubleRGB(img)

if ndims(img) == 2
    img = repmat(img, [1 1 3]);
end

if isa(img, 'double')
    newImg = img;
else
    newImg = double(img) / double(intmax(class(img)));
end
r = newImg(:, :, 1);
g = newImg(:, :, 2);
b = newImg(:, :, 3);

end


text:
% -------------------------------------------------------------------------- %
% hfssCoaxialCable(fid, Names, Axis, Center, Radii, Height, Units)
%
% Description :
% -------------
% This function creates the VBScript necessary to draw a Coaxial Cable in
% HFSS. The "Coaxial Cable" can have as many cylinders as you specify. This
% function only creates the geometric structure and does not set any material
% properties.
%
% Parameters :
% ------------
% fid - file identifier of the HFSS VBScript File.
% Name - a cell of strings that contains the names of each cylinder that is a
%        part of the co-axial cable (see example).
% Axis - axis of the Coaxial Cable (choose between 'X', 'Y' or 'Z').
% Center - Coordinates of the Center of the Coaxial Cable ([x, y, z]).
% Radii - an array of Radii of the Cylinders present in the Coaxial cable.
%         each radius corresponds to the respective name specfied in the
%         'Name' parameter.
% Height - length of the coaxial cable.
% Units - can be either 'in' or 'mm' or 'meter' or anything defined in HFSS.
%
% Example:
% ---------
% ...
% hfssCoaxialCable(fid, {'Cyl1_In', 'Cyl1_Er', 'Cyl1_Out'}, 'Z', [0, 0, 0],
%                  0.02, 0.03, 0.04], 10, 'in')
% ...
% -------------------------------------------------------------------------- %

% ----------------------------------------------------------------------------
% This file is part of HFSS-MATLAB-API.
%
% HFSS-MATLAB-API 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.
%
% HFSS-MATLAB-API 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
% Foobar; if not, write to the Free Software Foundation, Inc., 59 Temple
% Place, Suite 330, Boston, MA  02111-1307  USA
%
% Copyright 2004, Vijay Ramasami (rvc@ku.edu)
% ----------------------------------------------------------------------------
function hfssCoaxialCable(fid, Names, Axis, Center, Radii, Height, Units)

% Sort the Radii first in Ascending order.
[Radii, iR] = sort(Radii);
Names = cellstr(char(Names{iR}));

% Get the # of Cylinders.
nCylinders = length(Radii);

% First create the N-1 hollow cylinders.
for iR = nCylinders:-1:2
	hfssHollowCylinder(fid, Names{iR}, Axis, Center, Radii(iR-1), ...
	                   Radii(iR), Height, Units);
end

% Finally create the inner cylinder.
hfssCylinder(fid, Names{1}, Axis, Center, Radii(1), Height, Units);

author: yuip

repo: hfss-api

save_path: github-repos/MATLAB/yuip-hfss-api


text:
function x = i4vec_increment ( n, x )

%*****************************************************************************80
%
%% I4VEC_INCREMENT increments an I4VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input/output, integer X(N), the array to be incremented.
%
  x(1:n) = x(1:n) + 1;

  return
end


text:
function bbm = bbmean(x,B,w)
% BBMEAN  Bayesian bootstrap mean
%
%    Description
%    bbm = bbmean(x,B,w)
%    x NxM matrix of data
%    B number of bootstrap replicates
%    w Nx1 vector of weights
%
% Copyright (c) Aki Vehtari, 1998-2004

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

error('No mex-file for this architecture. See Matlab help and convert.m in ./linuxCsource or ./winCsource for help.')

author: gpstuff-dev

repo: gpstuff

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


text:
%% 1D Heat Transfer FDM
% by Manuel Diaz
clc; clear; close all;

%% Parameters
alpha = 0.25;       % difusion speed
iter = 400;         % iteration time steps

%% Domain Grid
L = 30;    % Length
n = 30;    % Nodes
dx = L/n; dt = 0.5;
x  = 1:dx:L;
t  = zeros(1,n);    % temperature of the nodes in our Domain

%% Wall Temperatura at x = 0
twall = 1;

%% Initial Condition
t_0 = zeros(1,n);   % @time = 0
t_0(1) = twall;     %Dirichlet BC
t_0(n) = t_0(n-1);   %Neumann BC

%% Main Loop
t_next = zeros(1,n); %Next time step
t = t_0;			 %Load I.C.
for k = 1:iter;
    for i = 2:n-1;
        t_next(i) = t(i) + (dt/dx^2)*alpha*(t(i+1)-2*t(i)+t(i-1));
    end
    % BC
    t_next(1) = twall;
	t_next(n) = t_next(n-1);
	% Update info
	t = t_next;
end

%% Make pretty figures
plot(x,t,'.'); xlabel 'x cell'; ylabel 'Temperature'; title '1D Heat Equation using FDM';


text:
% Test spectrogram

function tPlotSpectrogram

addpath('..');

[x, Fs] = wavread('addf8.wav');

Options = {'BW', 'NB', 'FLim', [500 1500]};
PlotSpectrogram(x, Fs, Options{:});
colormap(SpecColorMap);
colorbar;

% Use alternate colour map
figure;

Options = {'NSlice', 450, 'preF', 0.97,  'BW', 'WB'};
TLen = (length(x)-1)/Fs;
TLim = [0.3*TLen, 0.4*TLen];
PlotSpectrogram(x, TLim, Fs, Options{:});
colorbar;

% Sine wave test (expect sine wave peak at -6 dBov)
% Changing fc to 0, gives 0 dBov at dc
Fs = 16000;
fc = Fs/8;
NSamp = 5000;
t = (0:NSamp-1)/Fs;
Amax = 32767;
x = Amax * cos(2*pi*t*fc);

PmaxdBSPL = 92;             % Max level in dB SPL
PoffsdB = PmaxdBSPL + 20*log10(2);
g = 10^(PoffsdB/20);
AmaxN = Amax/g;

figure;
Options = {'Amax', AmaxN, 'FLim', [0 4000]};
PlotSpectrogram(x, Fs, Options{:});
colormap(SpecColorMap);
colorbar;

return


text:
function msfun_realtime_elapsed(block)
% Help for Writing Level-2 M-File S-Functions:
%   web([docroot '/toolbox/simulink/sfg/f7-67622.html']
%   http://www.mathworks.com/access/helpdesk/help/toolbox/simulink/sfg/f7-67622.html

% Copyright 2009, The MathWorks, Inc.

% instance variables
myRealTimeBaseline = 0;

setup(block);

%% ---------------------------------------------------
    function setup(block)
        % Register the number of ports.
        block.NumInputPorts  = 0;
        block.NumOutputPorts = 1;

        block.SetPreCompOutPortInfoToDynamic;
        block.OutputPort(1).Dimensions  = 1;
        block.OutputPort(1).SamplingMode = 'sample';

        % Set up the states
        block.NumContStates = 0;
        block.NumDworks = 0;

        % Register the parameters.
        block.NumDialogPrms     = 0; % scale factor

        % Block is fixed in minor time step, i.e., it is only executed on major
        % time steps. With a fixed-step solver, the block runs at the fastest
        % discrete rate.
        block.SampleTimes = [0 1];

        block.SetAccelRunOnTLC(false); % run block in interpreted mode even w/ Acceleration

        % methods called at run-time
        block.RegBlockMethod('Start', @Start);
        block.RegBlockMethod('Outputs', @Output);
        block.RegBlockMethod('SimStatusChange', @SimStatusChange);
    end

%%
    function Start(block)  %#ok<INUSD>
        myRealTimeBaseline = tic;
    end

%%
    function Output(block)
         block.OutputPort(1).Data = toc(myRealTimeBaseline);
    end

%%
    function SimStatusChange(block, status)         %#ok<INUSL>
        if status == 1,  % resume
            myRealTimeBaseline = tic;
        end
    end

end

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function [alpha, beta, mu, delta] = nigpar(m, v, s, k)
%NIGPAR Parameters for the Normal-Inverse-Gaussian distribution.
%   [ALPHA, BETA, MU, DELTA] = NIGPAR(M, V, S, K) returns the scale
%   paramter ALPHA, BETA which determines the skewness, the location
%   parameter MU and the scale parameter DELTA of the Normal-Inverse-Gaussian
%   distribution with mean M, variance V, skewness S and kurtosis K.
%
%   See also NIGPDF, NIGCDF, NIGINV, NIGRND, NIGSTATS.
%
%   References:
%   [1] Prause, K. (1999). The Generalized Hyperbolic Model

% -------------------------------------------------------------------------
%
% Allianz, Group Risk Controlling
% Risk Methodology
% Koeniginstr. 28
% D-80802 Muenchen
% Germany
% Internet:    www.allianz.de
% email:       ralf.werner@allianz.de
%
% Implementation Date:  2006 - 05 - 01
% Author:               Dr. Ralf Werner
%
% -------------------------------------------------------------------------

Tol = 1.0e-007;

%% Default values
if nargin < 1
    m = 0;
end
if nargin < 2
    v = 1;
end
if nargin < 4
    if nargin < 3
        s = 0;
        k = 3 + Tol;
    else
        k = 3 + 5/3*s^2 + Tol;
    end
end

%% Constraints for the parameters

if v <= 0
    error('The variance V must be positive.');
end

if (k - 5/3*s^2 - 3 <= 0)
    error('K must be greater than 3 + 5/3 S^2.');
end

alpha = sqrt((3*k - 4*s^2 - 9) / (v*(k-5/3*s^2 - 3)^2));
beta = s/(sqrt(v)*(k - 5/3*s^2 - 3));
mu = m - 3*s*sqrt(v)/(3*k - 4*s^2 - 9);
delta = 3^(3/2)*sqrt(v*(k - 5/3*s^2 - 3))/(3*k - 4*s^2 - 9);


text:
function y = wmean(x, w)
%WMEAN  Weighted average or mean value.
%
%  Description
%    WMEAN(X,W) is the weighted mean value of the elements in X
%    (along first dimension) given weights W.
%
%  See also wprctile
%
% Copyright (c) 2000-2013 Aki Vehtari

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

y=sum(bsxfun(@times,x,w),1);


text:
function smin = choose_smin(kernel,sn,prob)
%% Choose minimal spike spike to enforce sparsity constraint in spike inference
%  The function chooses a regularization weight for sparse deconvolution
%  with a given kernel and noise level. The weight of the regularizer
%  is chosen such that noise alone cannot lead to a non-zero solution is
%  at least prob.

% Inputs:
% kernel:       deconvolution kernel
%               if length(kernel) == 1 or length(kernel) == 2, then kernel
%               is treated as a set of discrete time constants g. Otherwise,
%               it is treated as the actual vector.
% sn:           noise level
% prob:         probability of zero solution (deafult: 0.99)

% Output:
% smin:          minimal spike size

% Author: Pengcheng Zhou, Carnegie Mellon University, 2016
% modified from choose_lambda.m Eftychios A. Pnevmatikakis, 2016, Simons Foundation
% based on ideas and discussions with J. Tubiana and G. Debregeas,
% Laboratorie Jean Parrin, UPMC, France

% Reference for this approach:
% Selesnick, I. (2012). Sparse deconvolution (an MM algorithm)

%%

if nargin < 3 || isempty(prob)
    prob = 0.99999;
end

if nargin < 2 || isempty(sn)
    sn = 1;
    warning('Noise value not provided. Using sn = 1...')
end

if length(kernel) <= 2
    kernel = filter(1,[1,-kernel(:)'],[1,zeros(1,999)]);
end

smin = sn/norm(kernel)*norminv(prob);

author: zhoupc


text:
function spm_eeg_review(D,flag,inv)
% General review (display) of SPM meeg object
% FORMAT spm_eeg_review(D,flags,inv)
%
% INPUT:
% D      - meeg object
% flag   - switch to any of the displays (optional)
% inv    - which source reconstruction to display (when called from
% spm_eeg_inv_imag_api.m)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Jean Daunizeau
% $Id: spm_eeg_review.m 7221 2017-11-16 14:25:37Z vladimir $

if nargin == 0
    [D, sts] = spm_select(1, 'mat$', 'Select M/EEG mat file');
    if ~sts, return; end
    D = spm_eeg_load(D);
end

%-- Initialize SPM figure
D.PSD.handles.hfig = spm_figure('GetWin','Graphics');
spm_clf(D.PSD.handles.hfig);
% Get default SPM graphics options --> revert back to defaults
D.PSD.SPMdefaults.col = get(D.PSD.handles.hfig,'colormap');
D.PSD.SPMdefaults.renderer = get(D.PSD.handles.hfig,'renderer');

%-- Create default userdata structure
try
    D.PSD.source.VIZU.current = inv;
end
[D] = PSD_initUD(D);
if ~strcmp(transformtype(D),'time') %
    D.PSD.type = 'epoched';
    D.PSD.trials.current = 1;
    D.PSD.VIZU.type = 2;
end

%-- Create figure uitabs
labels = {'EEG', 'MEG', 'MPLANAR', 'MCOMB', 'OTHER','info','source'};
callbacks = {'spm_eeg_review_callbacks(''visu'',''main'',''eeg'')',...
    'spm_eeg_review_callbacks(''visu'',''main'',''meg'')',...
    'spm_eeg_review_callbacks(''visu'',''main'',''megplanar'')',...
    'spm_eeg_review_callbacks(''visu'',''main'',''megcomb'')',...
    'spm_eeg_review_callbacks(''visu'',''main'',''other'')',...
    'spm_eeg_review_callbacks(''visu'',''main'',''info'')',...
    'spm_eeg_review_callbacks(''visu'',''main'',''source'')'};
try
    [h] = spm_uitab(D.PSD.handles.hfig,labels,callbacks,[],flag);
catch
    [h] = spm_uitab(D.PSD.handles.hfig,labels,callbacks,[],6);
end
D.PSD.handles.tabs = h;

% Add prepare and SAVE buttons
object.type = 'buttons';
object.list = 1;
D = spm_eeg_review_uis(D,object);
set(D.PSD.handles.BUTTONS.pop1,...
    'deletefcn',@back2defaults)

%-- Attach userdata to SPM graphics window
Dtmp = rmfield(D,'PSD');
D.PSD.D0 = Dtmp;
%CP, some how, removing and adding a field in the object doesn't work directly...
set(D.PSD.handles.hfig,...
    'units','normalized',...
    'color',[1 1 1],...
    'userdata',D);

try
    if ismac
        set(D.PSD.handles.hfig,'renderer','zbuffer');
    else
        set(D.PSD.handles.hfig,'renderer','OpenGL');
    end
catch
    set(D.PSD.handles.hfig,'renderer','OpenGL');
end

try
    switch flag
        case 1
            spm_eeg_review_callbacks('visu','main','eeg')
        case 2
            spm_eeg_review_callbacks('visu','main','meg')
        case 3
            spm_eeg_review_callbacks('visu','main','megplanar')
        case 4
            spm_eeg_review_callbacks('visu','main','megcomb')
        case 5
            spm_eeg_review_callbacks('visu','main','other')
        case 6
            spm_eeg_review_callbacks('visu','main','info')
        case 7
            spm_eeg_review_callbacks('visu','main','source')
    end
catch
    % Initilize display on 'info'
    spm_eeg_review_callbacks('visu','main','info')
end

%% Revert graphical properties of SPM figure back to normal
function back2defaults(e1,e2)
hf = spm_figure('FindWin','Graphics');
D = get(hf,'userdata');
try
    set(D.PSD.handles.hfig,'colormap',D.PSD.SPMdefaults.col);
    set(D.PSD.handles.hfig,'renderer',D.PSD.SPMdefaults.renderer);
end

%% initialization of the userdata structure
function [D] = PSD_initUD(D)
% This function initializes the userdata structure.

%-- Check spm_uitable capability (JAVA compatibility) --%
D.PSD.VIZU.uitable = spm_uitable;

%-- Initialize time window basic info --%
D.PSD.VIZU.xlim = [1,min([5e2,nsamples(D)])];
D.PSD.VIZU.info = 4; % show history
D.PSD.VIZU.fromTab = [];

%-- Initialize trials info --%
switch D.type

    %------ before epoching -----%

    case 'continuous'
        D.PSD.type = 'continuous';
        if ntrials(D) && ~isempty(events(D,1))
            Events = events(D);
            Nevents = length(Events);
            for i =1:Nevents
                if isempty(Events(i).duration)
                    Events(i).duration = 0;
                end
                if isempty(Events(i).value)
                    Events(i).value = '0';
                end
                if isempty(Events(i).type)
                    Events(i).type = '0';
                end
                if ~ischar(Events(i).value)
                    Events(i).value = num2str(Events(i).value);
                end
                if ~ischar(Events(i).type)
                    Events(i).type = num2str(Events(i).type);
                end
            end
            %D = events(D,1,Events);
        end
        D.PSD.VIZU.type = 1;

        %------ after epoching -----%

    case 'single'
        D.PSD.type = 'epoched';
        nTrials = D.ntrials;
        D.PSD.trials.TrLabels = cell(nTrials,1);
        bTrials = badtrials(D);
        for i = 1:nTrials
            if any(i==bTrials)
                str = ' (bad)';
            else
                str = ' (not bad)';
            end
            D.PSD.trials.TrLabels{i} = [...
                'Trial ',num2str(i),': ',char(conditions(D,i)) ,str];
        end
        D.PSD.trials.current = 1;
        D.PSD.VIZU.type = 1;
    case {'evoked','grandmean'}
        D.PSD.type = 'epoched';
        nTrials = D.ntrials;
        D.PSD.trials.TrLabels = cell(nTrials,1);
        for i = 1:nTrials
            D.PSD.trials.TrLabels{i} = [...
                'Trial ',num2str(i),' (average of ',...
                num2str(repl(D,i)),' events): ',...
                char(conditions(D,i))];
            D = events(D,i,[]);
        end
        D.PSD.trials.current = 1;
        D.PSD.VIZU.type = 1;
end

%-- Initialize channel info --%
nc = D.nchannels;
D.PSD.EEG.I  = indchantype(D,'EEG');
D.PSD.MEG.I  = sort(indchantype(D,'MEG'));
D.PSD.MEGPLANAR.I  = indchantype(D,'MEGPLANAR');
D.PSD.MEGCOMB.I  = indchantype(D,'MEGCOMB');
D.PSD.other.I = setdiff(1:nc, ...
    [D.PSD.EEG.I(:);D.PSD.MEG.I(:);D.PSD.MEGPLANAR.I(:);D.PSD.MEGCOMB.I(:)]);
%-- Get basic display variables (data range, offset,...)
if ~isempty(D.PSD.EEG.I)
    set(D.PSD.handles.hfig,'userdata',D);
    figure(D.PSD.handles.hfig)
    [out] = spm_eeg_review_callbacks('get','VIZU',D.PSD.EEG.I);
    D.PSD.EEG.VIZU = out;
else
    D.PSD.EEG.VIZU = [];
end
if ~isempty(D.PSD.MEG.I)
    set(D.PSD.handles.hfig,'userdata',D);
    figure(D.PSD.handles.hfig)
    [out] = spm_eeg_review_callbacks('get','VIZU',D.PSD.MEG.I);
    D.PSD.MEG.VIZU = out;
else
    D.PSD.MEG.VIZU = [];
end
if ~isempty(D.PSD.MEGPLANAR.I)
    set(D.PSD.handles.hfig,'userdata',D);
    figure(D.PSD.handles.hfig)
    [out] = spm_eeg_review_callbacks('get','VIZU',D.PSD.MEGPLANAR.I);
    D.PSD.MEGPLANAR.VIZU = out;
else
    D.PSD.MEGPLANAR.VIZU = [];
end
if ~isempty(D.PSD.MEGCOMB.I)
    set(D.PSD.handles.hfig,'userdata',D);
    figure(D.PSD.handles.hfig)
    [out] = spm_eeg_review_callbacks('get','VIZU',D.PSD.MEGCOMB.I);
    D.PSD.MEGCOMB.VIZU = out;
else
    D.PSD.MEGCOMB.VIZU = [];
end
if ~isempty(D.PSD.other.I)
    set(D.PSD.handles.hfig,'userdata',D);
    figure(D.PSD.handles.hfig)
    [out] = spm_eeg_review_callbacks('get','VIZU',D.PSD.other.I);
    D.PSD.other.VIZU = out;
else
    D.PSD.other.VIZU = [];
end

%-- Initialize inverse field info --%
if isfield(D,'inv') && ~isempty(D.inv)
    isInv = zeros(length(D.inv),1);
    for i=1:length(D.inv)
        if isfield(D.inv{i},'inverse') && ...
                isfield(D.inv{i}, 'method') && ...
                strcmp(D.inv{i}.method,'Imaging')
            isInv(i) = 1;
        end
    end
    isInv = find(isInv);
    Ninv = length(isInv);
    if Ninv>=1
        labels = cell(Ninv,1);
        callbacks = cell(Ninv,1);
        F  = zeros(Ninv,1);
        ID = zeros(Ninv,1);
        pst = [];
        for i=1:Ninv
            if ~isfield(D.inv{isInv(i)},'comment')
                D.inv{isInv(i)}.comment{1} = num2str(i);
            end
            if ~isfield(D.inv{isInv(i)},'date')
                D.inv{isInv(i)}.date(1,:) = '?';
                D.inv{isInv(i)}.date(2,:) = ' ';
            end
            if isfield(D.inv{isInv(i)}.inverse,'R2') ...
                    && isnan(D.inv{isInv(i)}.inverse.R2)
                D.inv{isInv(i)}.inverse.R2 = [];
            end
            if isfield(D.inv{isInv(i)}.inverse, 'ID')
                ID(i) = D.inv{isInv(i)}.inverse.ID;
            else
                ID(i) = nan;
            end
            labels{i} = [D.inv{isInv(i)}.comment{1}];
            callbacks{i} = ['spm_eeg_review_callbacks(''visu'',''inv'',',num2str(i),')'];
            try
                F(i) = D.inv{isInv(i)}.inverse.F;
                pst = [pst;D.inv{isInv(i)}.inverse.pst(:)];
            catch
                continue
            end
        end
        if isempty(pst)
            Ninv = 0;
        else
            pst = unique(pst);
        end
    end
else
    Ninv = 0;
end
if Ninv >= 1
    try
        if D.PSD.source.VIZU.current > Ninv
            D.PSD.source.VIZU.current = 1;
        end
    catch
        D.PSD.source.VIZU.current = 1;
    end
    D.PSD.source.VIZU.isInv = isInv;
    D.PSD.source.VIZU.pst = pst;
    D.PSD.source.VIZU.F = F;
    D.PSD.source.VIZU.ID = ID;
    D.PSD.source.VIZU.labels = labels;
    D.PSD.source.VIZU.callbacks = callbacks;
    D.PSD.source.VIZU.timeCourses = 1;
else
    D.PSD.source.VIZU.current = 0;
    D.PSD.source.VIZU.isInv = [];
    D.PSD.source.VIZU.pst = [];
    D.PSD.source.VIZU.F = [];
    D.PSD.source.VIZU.labels = [];
    D.PSD.source.VIZU.callbacks = [];
    D.PSD.source.VIZU.timeCourses = [];
end

author: spm


text:
% Note: no im2double() used !
img = imread('finger.png'); % Read image

img = imresize(img,0.25);    % Downscale image

% Get the valid region, this is a binary mask which indicates the region of
% the finger. For quick testing it is possible to use something like:
% fvr = ones(size(img));
% The lee_region() function can be found here:
% http://www.mathworks.com/matlabcentral/fileexchange/35752-finger-region-localisation
fvr = lee_region(im2double(img),4,20);    % Get finger region

%% Extract veins using wide line detector
r = 5; t=1; g=41; % Parameters
veins = huang_wide_line(img,fvr,r,t,g);

%% Visualise
rgb = zeros([size(img) 3]);
rgb(:,:,1) = im2double(img);
rgb(:,:,2) = im2double(img) + 0.2*double(veins);
rgb(:,:,3) = im2double(img);

figure;
subplot(2,1,1)
  imshow(img,[])
  title('Original image')
subplot(2,1,2)
  imshow(rgb)
  title('Wide line detector method')


text:
function t=talACPCPointsToInplane
% t=talACPCPointsToInplane
%
% Computes the coordinates of the talairach anchor points (ac,pc, sac
% etc..) in the coordinate frame of the inplanes. This is useful when
% converting all the data to analyze format for going into SPM
% Assumes that the talairach points exist in the usual anatomy directory and that mrSESSION.alignment exists.
% arw 110105 : Wrote it.
%

% Get voxel sizes to make sure that the transformation preserves volume
mrGlobals;

ipVoxSize = mrSESSION.inplanes.voxelSize;
volVoxSize = readVolAnatHeader(vANATOMYPATH);

% Transform ROI coordinates
xform = inv(mrSESSION.alignment);
talPoints = loadTalairachXform(mrSESSION.subject);
talPoints.refPoints

rp=talPoints.refPoints;

acpc=[talPoints.refPoints.acXYZ;talPoints.refPoints.pcXYZ];
disp(acpc)
acpc2=acpc(:,[2 1 3]);

[d,t]= xformROIcoords(acpc2',xform,volVoxSize,ipVoxSize)

disp('Note: Y axis is anterior->posterior, X axis is left to right on screen (right to left in brain)');
disp('AC Coords (y,x,slice) :');
disp(round(t(1:3,1)));
disp('PC Coords (y,x,slice) :');
disp(round(t(1:3,2)));

t=t(1:3,:);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ nxe, cnt ] = linop_normest( op, cmode, tol, maxiter )

%LINOP_NORMEST Estimates the operator norm.
%    EST = LINOP_NORMEST( OP ) estimates the induced norm of the operator:
%           || OP || = max_{||x||<=1} ||OP(X)||
%    using a simple power method similar to the MATLAB NORMEST functon.
%
%    When called with a single argument, LINOP_TEST begins with a real
%    initial vector. To test complex operators, use the two-argument
%    version LINOP_NORMEST( OP, cmode ), where:
%        cmode = 'R2R': real input, real output
%        cmode = 'R2C': real input, complex output
%        cmode = 'C2R': imag input, imag output
%        cmode = 'C2C': complex input, complex output
%
%    LINOP_NORMEST( OP, CMODE, TOL, MAXITER ) stops the iteration when the
%    relative change in the estimate is less than TOL, or after MAXITER
%    iterations, whichever comes first.
%
%    [ EST, CNT ] = LINOP_NORMEST( ... returns the estimate and the number
%    of iterations taken, respectively.

if isnumeric( op ),
    op = linop_matrix( op, 'C2C' );
elseif ~isa( op, 'function_handle' ),
    error( 'Argument must be a function handle or matrix.' );
end
x_real = true;
if nargin >= 2 && ~isempty( cmode ),
    switch upper( cmode ),
        case { 'R2R', 'R2C' }, x_real = true;
        case { 'C2R', 'C2C' }, x_real = false;
        otherwise, error( 'Invalid cmode: %s', cmode );
    end
end
if nargin < 3 || isempty( tol ),
    tol = 1e-8;
end
if nargin < 4 || isempty( maxiter ),
    maxiter = 50;
end
sz = op([],0);
if iscell( sz ),
    sz = sz{1};
elseif ~isempty(sz)
    sz = [sz(2),1];
else
    % if the input is the identity, it may not have a size associated with it,
    % so current behavior is to try an arbitrary size:
    sz = [50,1];
end
cnt = 0;
nxe = 0;
while true,
    if nxe == 0,
        if x_real,
            xx = randn(sz);
        else
            xx = randn(sz) + 1j*randn(sz);
        end
        nxe = sqrt( tfocs_normsq( xx ) );
    end
    yy = op( xx / max( nxe, realmin ), 1 );
    nye = sqrt( tfocs_normsq( yy ) );
    xx = op( yy / max( nye, realmin ), 2 );
    nxe0 = nxe;
    nxe = sqrt( tfocs_normsq( xx ) );
    if abs( nxe - nxe0 ) < tol * max( nxe0, nxe ),
        break;
    end
    cnt = cnt + 1;
    if cnt >= maxiter,
        break;
    end
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.


text:
% The COBRAToolbox: testFBA.m
%
% Purpose:
%     - tests the basic functionality of FBA
%       Tests four basic solution: Optimal minimum 1-norm solution, Optimal
%       solution on fructose, Optimal anaerobic solution, Optimal ethanol
%       secretion rate solution returns 1 if all tests were completed succesfully, 0 if not
%
% Authors:
%     - Original file: Joseph Kang 04/27/09
%     - CI integration: Laurent Heirendt January 2017
%

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testFBA'));
cd(fileDir);

% set the tolerance
tol = 1e-8;

% define the solver packages to be used to run this test
solverPkgs = {'ibm_cplex', 'mosek', 'gurobi', 'tomlab_cplex', 'glpk'};

% load the model
load('testFBAData.mat');
model = readCbModel('testFBAData.mat','modelName','model');
for k = 1:length(solverPkgs)

    % change the COBRA solver (LP)
    solverOK = changeCobraSolver(solverPkgs{k}, 'LP', 0);

    if solverOK == 1
        fprintf('   Testing flux balance analysis using %s ... ', solverPkgs{k});

        % check the optimal solution - BiomassEcoli
        fprintf('\n>> Optimal minimum 1-norm solution\n');
        model = changeObjective(model, {'BiomassEcoli'}, 1);
        solution = optimizeCbModel(model);

        % testing if f values are within range
        assert(abs(solution.f - solutionStd.f) < tol);

        % testing if c*x == f
        assert(abs(model.c' * solution.x - solution.f) < tol);

        % print the flux vector
        printFluxVector(model, solution.x, true, true);

        % check the optimal solution - fructose
        fprintf('\n>> Optimal solution on fructose\n');
        model2 = changeRxnBounds(model, {'EX_glc(e)', 'EX_fru(e)'}, [0 -9], 'l');
        solution2 = optimizeCbModel(model2);

        % testing if f values are within range
        assert(abs(solution2.f - solution2Std.f) < tol);

        % testing if c*x == f
        assert(abs(model2.c' * solution2.x - solution2.f) < tol);

        % print the flux vector
        printFluxVector(model2, solution.x, true, true);

        % check the optimal anaerobic solution
        fprintf('\n>> Optimal anaerobic solution\n');
        model3 = changeRxnBounds(model, 'EX_o2(e)', 0, 'l');
        solution3 = optimizeCbModel(model3);

        % testing if f values are within range
        assert(abs(solution3.f - solution3Std.f) < tol);

        % testing if c*x == f
        assert(abs(model3.c' * solution3.x - solution3.f) < tol);

        % check the optimal ethanol secretion rate solution
        fprintf('\n>> Optimal ethanol secretion rate solution \n');
        model4 = changeObjective(model, 'EX_etoh(e)', 1);
        solution4 = optimizeCbModel(model4);

        % testing if f values are within range
        assert(abs(solution4.f - solution4Std.f) < tol);

        % testing if c*x == f
        assert(abs(model4.c' * solution4.x - solution4.f) < tol);

        % output a success message
        fprintf('Done.\n');
    end
end

% change the directory
cd(currentDir)


text:
%KSOBEL Sobel edge detector
%
% K = KSOBEL() is the Sobel x-derivative kernel:
%         1/8  |1  0  -1|
%              |2  0  -2|
%              |1  0  -1|
%
% Notes::
% - This kernel is an effective vertical-edge detector
% - The y-derivative (horizontal-edge) kernel is K'
%
% See also ISOBEL.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB 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.
%
% MVTB 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 MVTB.  If not, see <http://www.gnu.org/licenses/>.

function k = ksobel

    k = [   1 0 -1
            2 0 -2
            1 0 -1] / 8;

author: petercorke

repo: machinevision-toolbox-matlab


text:
function rx = rxShowStats(rx,flag);
%
% rx = rxShowStats([rx],flag);
%
% mrRx:
% Show or Hide the text controls
% for the statistics of comparing
% the volume and reference images/vols.
%
% fig: figure to add/remove fields from.
%
% flag: if 1, add fields; if 0, remove.
%
%
% 03/05 ras.
if ieNotDefined('rx')
    cfig = findobj('Tag','rxControlFig');
    rx = get(cfig,'UserData');
end

if flag==1
    % move the axes over to make room
    set(rx.ui.compareAxes,'Position',[0.05 .25 .9 .7]);

    % show/create fields
    if isfield(rx.ui,'compareStats') & ishandle(rx.ui.compareStats.corrcoefVal)

        % show
        fields = fieldnames(rx.ui.compareStats);
        for i = 1:length(fields)
            set(rx.ui.compareStats.(fields{i}),'Visible','on');
        end
    else
        %%%%% create %%%%%

        % create text for corr coefficient (R)
        h = uicontrol('Style','text',...
                      'Units','Normalized',...
                      'Position',[0 .16 .25 .05],...
                      'String','Correlation [R]:',...
                      'FontWeight','bold',...
                      'FontAngle','italic',...
                      'HorizontalAlignment','left',...
                      'BackgroundColor',get(gcf,'Color'),...
                      'ForegroundColor',[0 .4 0]);
        rx.ui.compareStats.corrcoefLabel = h;
        h = uicontrol('Style','text',...
                      'Units','Normalized',...
                      'Position',[.25 .16 .2 .05],...
                      'String','(N/A)',...
                      'HorizontalAlignment','left',...
                      'BackgroundColor',get(gcf,'Color'),...
                      'ForegroundColor',[0 0 0]);
        rx.ui.compareStats.corrcoefVal = h;

        % create text for root mean square error
        h = uicontrol('Style','text',...
                      'Units','Normalized',...
                      'Position',[0 .1 .25 .05],...
                      'String','R.M.S. Error:',...
                      'FontWeight','bold',...
                      'FontAngle','italic',...
                      'HorizontalAlignment','left',...
                      'BackgroundColor',get(gcf,'Color'),...
                      'ForegroundColor',[.5 0 0]);
        rx.ui.compareStats.rmseLabel = h;
        h = uicontrol('Style','text',...
                      'Units','Normalized',...
                      'Position',[.25 .1 .2 .05],...
                      'String','(N/A)',...
                      'HorizontalAlignment','left',...
                      'BackgroundColor',get(gcf,'Color'),...
                      'ForegroundColor',[0 0 0]);
        rx.ui.compareStats.rmseVal = h;
    end
else
    % hide 'em
    fields = fieldnames(rx.ui.compareStats);
    for i = 1:length(fields)
        set(rx.ui.compareStats.(fields{i}),'Visible','off');
    end

    % restore axes to full size
	set(rx.ui.compareAxes,'Position',[0.05 .1 .9 .85]);
end

if ishandle(rx.ui.controlFig)
    set(rx.ui.controlFig,'UserData',rx);
end

rxRefresh(rx);

return


text:
function estimateClinicalRF_Aerts(pathRF,nSplit,nBoot,seed)

startpath = pwd;

% INITIALIZATIONS
cd(pathRF), load('training')
testSplit = 1/3; % Proportion of test cases in the splits
testCost = 0.4:0.1:2.5; % Emphasis factor on positive instances during random forest training
comb = {[1,2,3],[1,2,4],[1,2,5],[1,2,3,4]};

% COMPUTATIONS
[indClinic,bestCost] = estimateClinicalRF(pathRF,'DeathSign','CT',comb,testCost,nSplit,testSplit,nBoot,seed,0);
training.clinical.bestAdd.DeathSign.CT = indClinic;
training.cost.DeathSign.CT = bestCost;

cd(pathRF), save('training','training')
cd(startpath)
end

author: mvallieres

repo: radiomics


text:
%

dbName = '/home/admin/databases/PLUTONS/wf/plutons';
ds = datasource('antelope',dbName);

%t1 = datenum('2010/2/27 6:35');
%t2 = datenum('2010/2/27 6:55');
%filterObj = filterobject('b',[0.033 0.2],3);
%stachan = {
% 'UTLO'   'BHZ'
% 'UTKH'   'BHZ'
% 'UTSA'   'BHZ'
% 'UTQU'   'BHZ'
% 'UTLV'   'BHZ'
% 'UTSM'   'BHZ'
% 'UTLL'   'SHZ'
% 'UTSS'   'SHZ'
% 'UTLA2'  'SHZ'
% 'UTSW'   'SHZ'
% 'UTMK'   'SHZ'
% 'UTCM'   'SHZ'
% 'UTCA'   'SHZ'
% 'UTZN'   'SHZ'
% 'UTTM'   'SHZ'}

%t1 = datenum('2009/9/22 22:47:10');
%t2 = datenum('2009/9/22 22:47:30');
t1 = datenum('2010/02/27 22:47:10');
t2 = datenum('2010/02/27 22:47:30');

filterObj = filterobject('b',[0.5 10],2);
stachan = {
'UTSM'   'BHZ'
'UTQU'   'BHZ'
'UTTM'   'SHZ'
'UTLA2'  'SHZ'
'UTCA'   'SHZ'
'UTLV'   'BHZ'
'UTCM'   'SHZ'
'UTSW'   'SHZ'
'UTSA'   'BHZ'
'UTZN'   'SHZ'
'UTSS'   'SHZ'
'UTMK'   'SHZ'
'UTLO'   'BHZ'
'UTLL'   'SHZ'
'UTKH'   'BHZ'
}

for n = 1:size(stachan,1)
    scnl = scnlobject(stachan(n,1),stachan(n,2),'','');
    try
        w(n) = waveform(ds,scnl,t1,t2);
    catch
        w(n) = waveform;
    end
end
f = find(get(w,'DATA_LENGTH')~=0);
w = w(f);

% APPLY RESPONSES TO ORIGINAL WAVEFORMS (SHORT PERIOD)

wFilt = filtfilt(filterObj,w);
wNew = response_apply(wFilt,filterObj,'antelope',dbName);

% for n =1:numel(w)
%     n
%     wFilt(n) = filtfilt(filterObj,w(n));
%     wNew(n) = response_apply(wFilt(n),filterObj,dbName);
% end
save

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% check one channel
n = 3;
figure
plot(wFilt(n),'r');
hold on;
plot(wNew(n),'b');
response_plot(get(wNew(n),'RESPONSE'),[0.01 50])

load
% compare all channels
cNew = correlation(wNew);
%cNew = xcorr(cNew);
%cNew = adjusttrig(cNew);
plot(cNew,'wig',0.6)
title('Filtered and corrected for instrument response','FontSize',16);
%print(gcf,'-dpsc2','demo_corrected.ps');

cFilt = correlation(wFilt);
plot(cFilt,'wig',.6)
title('Filtered only (30.5 - 10 Hz)','FontSize',16);
%print(gcf,'-dpsc2','demo_filtered.ps');

c = correlation(w);
plot(c,'wig',.6)
title('Raw','FontSize',16);
%print(gcf,'-dpsc2','demo_raw.ps');

%
% % compare all channels
% cNew = correlation(wNew);
% cNew = crop(cNew,[0 600]);
% cNew = xcorr(cNew,[100 300]);
% cNew = adjusttrig(cNew);
% plot(cNew,'wig',0.6)
% title('Filtered and corrected for instrument response','FontSize',16);
% print(gcf,'-dpsc2','demo_corrected.ps');
%
% cFilt = correlation(wFilt);
% cNew = crop(cFilt,[0 600]);
% plot(cFilt,'wig',.6)
% title('Filtered only (30-5 s)','FontSize',16);
% print(gcf,'-dpsc2','demo_filtered.ps');
%
% c = correlation(w);
% cNew = crop(c,[0 600]);
% plot(c,'wig',.6)
% title('Raw','FontSize',16);
% print(gcf,'-dpsc2','demo_raw.ps');
%
%

author: geoscience-community-codes

repo: GISMO


text:
%Simple function to test the cross method
d=10;
%elem_fun=@(x) sum(x); %Just sum of everything
p=0:d-1; p = 2.^p;
a=-5;b=5;
n=2^d;
h=(b-a)/(n-1);
%mv=@(x) x.^3;
%elem_fun=@(x) 1.0./(dot((x-1),p)+1e-3); %Just sum of everything
%elem_fun=@(x) mv(1e-12+dot(x-1,p)*h);

%Compare functions of TT-tensors
x=tt_x(d,2); x=tt_tensor(x);
e=tt_ones(d,2);e=tt_tensor(e);
x=a*e+h*x; x=round(x,1e-13);
rs=x;
%fun=@(x) 1.0./sqrt(x);
%fun=@(x) 1.0./x;
fun=@(x) exp(-(x).^4) + 1;
elem_fun = @(ind) fun(rs(ind));
%elem_fun=@(ind) rs(ind);
%elem_fun=@(x) sqrt(x(1))+x(2);

eps=1e-6;
%y=tt_rc2(2*d,2,elem_fun,1e-12);
y=tt_rc(d,2,elem_fun,1e-6,'nswp',40,'change_dir_on',false);

z=funcrs2(rs,fun,1e-12,rs,20);

author: oseledets

repo: TT-Toolbox


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

nsteps = 10;
Vfrac_final = 0.1;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

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

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function Bcol = perform_segmentation_colorization(B)

% perform_segmentation_colorization - return a color image
%
%   Bcol = perform_segmentation_colorization(B);
%
%   Copyright (c) 2007 Gabriel Peyre

C = [[1; 0; 0] [0;1;0] [0;0;1] [1;1;0] [1;0;1] [0;1;1] [.3;.7;1]]';
ncol = size(C,1);
Bcol = C(cat(3,B,B+ncol,B+2*ncol));

author: gpeyre

repo: matlab-toolboxes


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open 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.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

function T = fixT(T)
assert(all(size(T) == [4 4]));
T(1:3, 1:3) = fixR(T(1:3, 1:3));
assert(all(T(4, :) == [0 0 0 1]));
end

author: uzh-rpg


text:
function protocol = FSL2Protocol(bvalfile, bvecfile)
%
% function protocol = FSL2Protocol(bvalfile, bvecfile)
%
% Note: for NODDI, the exact sequence timing is not important.
%  this function reverse-engineerings one possible sequence timing
%  given the b-values.
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

protocol.pulseseq = 'PGSE';
protocol.schemetype = 'multishellfixedG';
protocol.teststrategy = 'fixed';

% load bval
bval = load(bvalfile);
bval = bval';

% set total number of measurements
protocol.totalmeas = length(bval);

% set the b=0 indices
protocol.b0_Indices = find(bval==0);
protocol.numZeros = length(protocol.b0_Indices);

% find the unique non-zero b-values
B = unique(bval(bval>0));

% set the number of shells
protocol.M = length(B);
for i=1:length(B)
    protocol.N(i) = length(find(bval==B(i)));
end

% maximum b-value in the s/mm^2 unit
maxB = max(B);

% set maximum G = 40 mT/m
Gmax = 0.04;

% set smalldel and delta and G
GAMMA = 2.675987E8;
tmp = nthroot(3*maxB*10^6/(2*GAMMA^2*Gmax^2),3);
for i=1:length(B)
    protocol.udelta(i) = tmp;
    protocol.usmalldel(i) = tmp;
    protocol.uG(i) = sqrt(B(i)/maxB)*Gmax;
end

protocol.delta = zeros(size(bval))';
protocol.smalldel = zeros(size(bval))';
protocol.G = zeros(size(bval))';

for i=1:length(B)
    tmp = find(bval==B(i));
    for j=1:length(tmp)
        protocol.delta(tmp(j)) = protocol.udelta(i);
        protocol.smalldel(tmp(j)) = protocol.usmalldel(i);
        protocol.G(tmp(j)) = protocol.uG(i);
    end
end

% load bvec
bvec = load(bvecfile);
protocol.grad_dirs = bvec';

% make the gradient directions for b=0's [1 0 0]
for i=1:length(protocol.b0_Indices)
    protocol.grad_dirs(protocol.b0_Indices(i),:) = [1 0 0];
end

% make sure the gradient directions are unit vectors
for i=1:protocol.totalmeas
    protocol.grad_dirs(i,:) = protocol.grad_dirs(i,:)/norm(protocol.grad_dirs(i,:));
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function geo_params = drawparticals(geo_param, pf_param)
geo_param = repmat(geo_param, [1,pf_param.p_num]);
geo_params = geo_param + randn(6,pf_param.p_num).*repmat(pf_param.affsig(:),[1,pf_param.p_num]);
end

repo: FCNT


text:
function x = calcAxisDistribution(cs,varargin)
% compute the axis distribution of an uniform ODF or MDF
%
% Syntax
%   value = calcAxisDistribution(cs,a)
%   adf = calcAxisDistribution(cs)
%
% Input
%  cs - @crystalSymmetry
%  h  - @vector3d
%
% Output
%  value - values of the axis distribution function at axes a
%  adf - axes distribution function @S2Fun
%
% See also
% SO3Fun/calcAxisDistribution

[oR,dcs,nSym] = fundamentalRegion(cs,varargin{:});
varargin = delete_option(varargin,'complete');
if isa(varargin{1},'symmetry'), varargin(1) = []; end

if ~isempty(varargin) && isa(varargin{1},'vector3d')

  x = getValue(varargin{1});

else

  f = @(h) getValue(h);
  x = S2FunHarmonicSym.quadrature(f,dcs,'bandwidth',256,varargin{:});

end

function value = getValue(h)
  h = project2FundamentalRegion(h,dcs);
  omega = oR.maxAngle(h);
  value = nSym * (omega - sin(omega)) ./ pi;
end

end

author: mtex-toolbox

repo: mtex


text:
function [ vWm, vWc, scalingFactor ] = CalcSigmaPointsWeights( paramAlpha, paramBeta, paramKappa, stateOrder )
% ----------------------------------------------------------------------------------------------- %
% [ vWm, vWc, scalingFactor ] = CalcSigmaPointsWeights( paramAlpha, paramBeta, paramKappa, stateOrder )
%   Calculates the weights for the mean, covariance and scaling factor of
%   the Unscented Transform.
% Input:
%   - paramAlpha    -   Parameter Alpha.
%                       Influences how far the Sigma Points are form the
%                       mean.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: (0, 1].
%   - paramBeta     -   Parameter Beta.
%                       For Gaussian PDF the optimal value is 2.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: NA.
%   - paramKappa    -   Parameter Kappa.
%                       Influences how far the Sigma Points are form the
%                       mean.
%                       Type: 'Single' / 'Double'.
%                       Range: [0, inf).
%   - stateOrder    -   State Vector  Dimension.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: {1, 2, ...}.
% Output:
%   - vWm           -   Mean Weights Vector.
%                       Weights used for calculation of the mean of the
%                       transformed data.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - vWc           -   Covariance Weights Vector.
%                       Weights used for calculation of the covariance of
%                       the transformed data.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - scalingFctr   -   Scaling Factor.
%                       Parameter of the Generalized Unscented Transform
%                       which controls
%                       Equals to 'size(hF(mX(:, 1), 1)'.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: {1, 2, ...}.
% References
%   1.  Unscented Transform (Wikipedia) - https://en.wikipedia.org/wiki/Unscented_transform.
%   2.  Lecture 5: Unscented Kalman Filter and General Gaussian Filtering (Simo Sarkka).
%   3.  Robot Mapping - Unscented Kalman Filter (Cyrill Stachniss).
% Remarks:
%   1.  This method of parameters is often called Generalized / Scaled
%       Unscented Kalman Filter.
% TODO:
%   1.  U.
% Release Notes:
%   -   1.0.000     24/08/2018  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

numSigmaPoints  = (2 * stateOrder) + 1;

vWm = zeros(numSigmaPoints, 1);
vWc = zeros(numSigmaPoints, 1);

paramLambda = ((paramAlpha * paramAlpha) * (stateOrder + paramKappa)) - stateOrder;

vWm(1)      = paramLambda / (stateOrder + paramLambda);
vWm(2:end)  = 1 / (2 * (stateOrder + paramLambda));
vWc(1)      = vWm(1) + (1 - (paramAlpha * paramAlpha) + paramBeta);
vWc(2:end)  = vWm(2:end);

scalingFactor = sqrt(stateOrder + paramLambda);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function [L, U, P, Q] = gecp(A)
%GECP calculate Gauss elimination with complete pivoting
%
% (G)aussian (E)limination (C)omplete (P)ivoting
% Input : A nxn matrix
% Output
% L = Lower triangular matrix with ones as diagonals
% U = Upper triangular matrix
% P and Q permutations matrices so that P*A*Q = L*U
%
% See also LU
%
% written by : Cheilakos Nick
[n, n] = size(A);
p = 1:n;
q = 1:n;
for k = 1:n-1
    [maxc, rowindices] = max( abs(A(k:n, k:n)) );
    [maxm, colindex] = max(maxc);
    row = rowindices(colindex)+k-1; col = colindex+k-1;
    A( [k, row], : ) = A( [row, k], : );
    A( :, [k, col] ) = A( :, [col, k] );
    p( [k, row] ) = p( [row, k] ); q( [k, col] ) = q( [col, k] );
    if A(k,k) == 0
      break
    end
    A(k+1:n,k) = A(k+1:n,k)/A(k,k);
    i = k+1:n;
    A(i,i) = A(i,i) - A(i,k) * A(k,i);
end
L = tril(A,-1) + eye(n);
U = triu(A);
   P = eye(n);
   P = P(p,:);
   Q = eye(n);
   Q = Q(:,q);

author: Sable

repo: mcbench-benchmarks


text:
function [ x, rank, seed ] = mono_between_random ( m, n1, n2, seed )

%*****************************************************************************80
%
%% MONO_BETWEEN_RANDOM: random monomial with total degree between N1 and N2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N1, N2, the minimum and maximum degrees.
%    0 <= N1 <= N2.
%
%    Input, integer SEED, the random number seed.
%
%    Output, integer X(M), the random monomial.
%
%    Output, integer RANK, the rank of the monomial.
%
%    Output, integer SEED, the updated random number seed.
%
  n1_copy = max ( n1, 0 );
  rank_min = mono_upto_enum ( m, n1_copy - 1 ) + 1;
  rank_max = mono_upto_enum ( m, n2 );
  [ rank, seed ] = i4_uniform_ab ( rank_min, rank_max, seed );
  x = mono_unrank_grlex ( m, rank );

  return
end

author: johannesgerer

repo: jburkardt-m