text:
function varargout = panel_guidelines( varargin )
% PANEL_GUIDELINES: Load a scenario in the guidelines tab.
%
% USAGE:  [bstPanel] = panel_guidelines('CreatePanel', ScenarioName)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2017

eval(macro_method);
end

%% ===== CREATE PANEL =====
function [bstPanelNew, panelName] = CreatePanel(ScenarioName) %#ok<DEFNU>
    % Java initializations
    import java.awt.*;
    import javax.swing.*;
    panelName = 'Guidelines';

    % Load all the panels from the target scenario
    switch lower(ScenarioName)
        case 'epileptogenicity'
            ctrl = scenario_epilepto('CreatePanels');
        otherwise
            error(['Unknow scenario: ' ScenarioName]);
    end
    % Create main panel
    ctrl.jPanelContainer = gui_component('Panel');
    ctrl.ScenarioName = ScenarioName;

    % Create control panel
    jPanelControl = gui_river([1,1], [1,3,1,3], ' ');
    ctrl.jPanelContainer.add(jPanelControl, BorderLayout.EAST);
    % Add buttons
    buttonFormat = {Insets(0,0,0,0), Dimension(java_scaled('value', 47), java_scaled('value', 22))};
    ctrl.jButtonReset = gui_component('button', jPanelControl, 'br', 'Reset', buttonFormat, [], @(h,ev)ResetPanel());
    gui_component('label', jPanelControl, 'br', ' ');
    ctrl.jButtonPrev = gui_component('button', jPanelControl, 'br', '<<', buttonFormat, [], @(h,ev)SwitchPanel('prev'));
    ctrl.jButtonNext = gui_component('button', jPanelControl, 'br', '>>', buttonFormat, [], @(h,ev)SwitchPanel('next'));
    ctrl.jLabelStep  = gui_component('label', jPanelControl, 'br hfill', sprintf('0 / %d', length(ctrl.jPanels)));
    ctrl.jLabelStep.setHorizontalAlignment(JTextField.CENTER);
    gui_component('label', jPanelControl, 'br', ' ');
    ctrl.jButtonSkip = gui_component('button', jPanelControl, 'br', 'Skip', buttonFormat, [], @(h,ev)SwitchPanel('skip'));

    % Create the BstPanel object that is returned by the function
    bstPanelNew = BstPanel(panelName, ...
                           ctrl.jPanelContainer, ...
                           ctrl);
end

%% =================================================================================
%  === EXTERNAL CALLBACKS  =========================================================
%  =================================================================================
%% ===== GET PANEL CONTENTS =====
% GET Panel contents in a structure
function s = GetPanelContents() %#ok<DEFNU>
    s = [];
end

%% ===== CLOSE PANEL =====
function ClosePanel(varargin) %#ok<DEFNU>
    panelName = 'Guidelines';
    % Hide panel
    gui_hide(panelName);
    % Release mutex
    bst_mutex('release', panelName);
end

%% ===== GET CURRENT PANEL =====
function iPanel = GetCurrentPanel()
    % Get panel controls handles
    ctrl = bst_get('PanelControls', 'Guidelines');
    if isempty(ctrl)
        iPanel = 0;
        return;
    end
    % Read panel index
    strIndex = char(ctrl.jLabelStep.getText());
    iPanel = str2num(strIndex(1:2));
end

%% ===== SWITCH PANEL =====
function SwitchPanel(command)
    import java.awt.*;
    % Get panel controls handles
    ctrl = bst_get('PanelControls', 'Guidelines');
    if isempty(ctrl)
        return;
    end
    % Get current panel
    iPanel = GetCurrentPanel();

    % === VALIDATE CURRENT PANEL ===
    % If there is a validation callback for this panel
    if (iPanel >= 1) && ~isempty(ctrl.fcnValidate{iPanel}) && strcmpi(command, 'next')
        [isValidated, errMsg] = ctrl.fcnValidate{iPanel}();
        if ~isempty(errMsg)
            bst_error(errMsg, sprintf('Step #%d', iPanel), 0);
            return;
        end
    end
    % If this is the last panel and moving forward: stop
    if (iPanel == length(ctrl.jPanels)) && strcmpi(command, 'next')
        return;
    end

    % === MOVE TO NEXT PANEL ===
    % Remove existing panel
    if (iPanel >= 1) && (iPanel <= length(ctrl.jPanels))
        ctrl.jPanelContainer.remove(ctrl.jPanels(iPanel));
    end
    % Switch according to command
    switch (command)
        case 'next',  iPanel = iPanel + 1;
        case 'skip',  iPanel = iPanel + 1;
        case 'prev',  iPanel = iPanel - 1;
    end
    % If invalid panel: stop
    if (iPanel < 1) || (iPanel > length(ctrl.jPanels))
        return;
    end
    % Add new panel
    ctrl.jPanelContainer.add(ctrl.jPanels(iPanel), BorderLayout.CENTER);
    % Update step number
    ctrl.jLabelStep.setText(sprintf('%d / %d', iPanel, length(ctrl.jPanels)));
    % Repaint
    ctrl.jPanelContainer.invalidate();
    ctrl.jPanelContainer.repaint();
    % First panel: Disable previous button
    if (iPanel == 1)
        ctrl.jButtonPrev.setEnabled(0);
    else
        ctrl.jButtonPrev.setEnabled(1);
    end
    % Show/Hide Skip button
    ctrl.jButtonSkip.setVisible(ctrl.isSkip(iPanel));
%     % Last panel: Disable next button
%     if (iPanel == length(ctrl.jPanels))
%         ctrl.jButtonNext.setEnabled(0);
%     else
%         ctrl.jButtonNext.setEnabled(1);
%     end

    % === UPDATE NEW PANEL ===
    % If there is an update callback for this panel
    if (iPanel >= 1) && ~isempty(ctrl.fcnUpdate{iPanel})
        ctrl.fcnUpdate{iPanel}();
    end
end

%% ===== RESET PANEL =====
function ResetPanel()
    % Get panel controls handles
    ctrl = bst_get('PanelControls', 'Guidelines');
    if isempty(ctrl)
        return;
    end
    % Get current panel
    iPanel = GetCurrentPanel();
    if (iPanel >= 1) && ~isempty(ctrl.fcnReset{iPanel})
        ctrl.fcnReset{iPanel}();
    end
end

%% ===== OPTIONS: PICK FILE CALLBACK =====
function [OutputFiles, FileFormat] = PickFile(jControl, DefaultDir, SelectionMode, FilesOrDir, Filters, DefaultFormat) %#ok<DEFNU>
    % Parse inputs
    if (nargin < 6) || isempty(DefaultFormat)
        DefaultFormat = [];
    end
    % Get default import directory and formats
    LastUsedDirs = bst_get('LastUsedDirs');
    DefaultFormats = bst_get('DefaultFormats');
    % Default dir type
    DefaultFile = LastUsedDirs.(DefaultDir);
    % Default filter
    if ~isempty(DefaultFormat) && isfield(DefaultFormats, DefaultFormat)
        defaultFilter = DefaultFormats.(DefaultFormat);
    else
        defaultFilter = [];
        DefaultFormat = [];
    end
    % Pick a file
    [OutputFiles, FileFormat, FileFilter] = java_getfile('open', 'Select file', DefaultFile, SelectionMode, FilesOrDir, Filters, defaultFilter);
    % If nothing selected
    if isempty(OutputFiles)
        return
    end
    % If only one file selected
    if ~iscell(OutputFiles)
        OutputFiles = {OutputFiles};
    end
    % Get the files
    OutputFiles = file_expand_selection(FileFilter, OutputFiles);
    if isempty(OutputFiles)
        error(['No ' FileFormat ' file in the selected directories.']);
    end
    % Save default import directory
    if ischar(OutputFiles)
        newDir = OutputFiles;
    elseif iscell(OutputFiles)
        newDir = OutputFiles{1};
    end
    % Get parent folder if needed
    if ~isdir(newDir)
        newDir = bst_fileparts(newDir);
    end
    LastUsedDirs.(DefaultDir) = newDir;
    bst_set('LastUsedDirs', LastUsedDirs);
    % Save default import format
    if ~isempty(DefaultFormat)
        DefaultFormats.(DefaultFormat) = FileFormat;
        bst_set('DefaultFormats',  DefaultFormats);
    end

    % Get file descriptions (one/many)
    if ischar(OutputFiles)
        strFiles = OutputFiles;
    else
        if (length(OutputFiles) == 1)
            strFiles = OutputFiles{1};
        elseif isa(jControl, 'javax.swing.JLabel')
            strFiles = sprintf('%s
', OutputFiles{:});
        else
            strFiles = sprintf('[%d files]', length(OutputFiles));
        end
    end

    % Update the attached control
    if isempty(jControl)
        %disp(strFiles);
    elseif isa(jControl, 'javax.swing.JTextField')
        jControl.setText(strFiles);
    elseif isa(jControl, 'javax.swing.JLabel')
        jControl.setText(['' strFiles]);
    end
end

author: brainstorm-tools


text:
%% Example
% This script is used to initialize the hand guiding functionality on KUKA
% robot remotely

% Mohammad SAFEEA, 9th of June 2017

% TO use this example:
% 1- Start the server on KUKA
% 2- Start the client on Matlab
% 3- Press the white button and hand guide the robot
% 4- To save the coordinates of the robot, long click on the green button
% 5- Repeat three times
% 6- The robot gives you an interval of five seconds, clear the area around
% the robot directly, the robot will reproduce the motion tought.

fprintf('Connecting to server')

ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
t=net_establishConnection( ip );

if ~exist('t','var') || isempty(t)
  warning('Connection could not be establised, script aborted');
  return;
else

startHandGuiding( t )

p1 = getJointsPos( t );

startHandGuiding( t )

p2 = getJointsPos( t );

startHandGuiding( t )

p3 = getJointsPos( t );

startHandGuiding( t )

fprintf('CLear the area round the robot, the robot is going to move');

pause(5)

relVel=0.1;
movePTPJointSpace( t , p1, relVel)
movePTPJointSpace( t , p2, relVel)
movePTPJointSpace( t , p3, relVel)

      %% turn off the server
       net_turnOffServer( t );

       fclose(t);
end

repo: KST-Kuka-Sunrise-Toolbox


text:
 function [xhat, smap, ref] = ir_mri_coil_combine(ykj, varargin)
%function [xhat, smap, ref] = ir_mri_coil_combine(ykj, varargin)
%|
%| Given multiple MRI surface coil images, estimate their sensitivity maps,
%| relative to the sqrt-sum-of-squares (SSoS) map, and then perform the
%| "best" linear combination using the sensitivity maps.
%| todo: currently ignores noise correlations
%|
%| This routine provides an alternative to the usual "sqrt-sum-of-squares"
%| approach to combining multiple coil data.  It produces a complex-valued
%| output image that has approximately a complex gaussian distribution at
%| each voxel.  However, it does not preserve image phase information
%| because the image phase is included in the sensitivity map phase estimates.
%| So it might be logical to just take the real part of the output image.
%| The output is "approximately" gaussian only because of random errors in
%| the sensitivity maps estimated within.
%|
%| The variance is space-dependent because of coil sensitivity variations.
%| todo: return variance map!
%|
%| Signal model: y_kj = s_kj f_j + noise_kj
%|	s_kj	sensitivity map (relative to SSoS reference image)
%|	f_j	unknown underlying object
%|	k: coil index
%|	j: voxel index
%|
%| in
%|	ykj	[(N) ncoil]	noisy complex images (2D or 3D) for each coil
%|
%| options
%|	bodycoil [(N)]		reference image (optional) (default: SSoS)
%|	thresh0			zero reference image values below this fraction
%|				 of peak value (default: 0.05)
%|	thresh1			fraction of reference peak used for median
%|				 initial value in "background" (default: 0.05)
%|
%| out
%|	xhat	[(N)]		coil combination image
%|	smap	[(N) ncoil]	estimated sensitivity maps
%|
%| Uses the 'chol' (Cholesky) option of mri_sensemap_denoise.m which is
%| fast for medium sized images but may use much memory for large images.
%|
%| Copyright 2014-08-18, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(ykj, 'test')
	ir_mri_coil_combine_test
return
end

% defaults
arg.thresh0 = 0.05; % below this set reference image to 0
arg.thresh1 = 0.05; % below this set reference image to median
arg.bodycoil = [];
arg.chol = true;

arg = vararg_pair(arg, varargin);

tmp = size(ykj);
NN = tmp(1:end-1);
ncoil = tmp(end);

% default reference image
if isempty(arg.bodycoil)
	ref = ir_mri_coil_combine_ref(ykj, arg.thresh0);
else
	ref = arg.bodycoil;
end

smap = mri_sensemap_denoise(ykj, 'thresh', arg.thresh1, ...
	'bodycoil', ref, 'chol', 1, 'niter', 1);

coildim = ndims(ykj);
numer = sum(conj(smap) .* ykj, coildim);
denom = sum(abs(smap).^2, coildim);
xhat = div0(numer, denom); % LS solution at each pixel

% ir_mri_coil_combine_ref()
function ref = ir_mri_coil_combine_ref(ykj, thresh0)
coildim = ndims(ykj);
ref = sqrt(sum(abs(ykj).^2, coildim)); % sqrt-sum-of-squares
peak = max(abs(ref(:)));
ref(abs(ref) < peak * thresh0) = 0; % zero low SNR pixels
%ref = ref / peak; % normalize to 1

% ir_mri_coil_combine_test()
% built-in test/example
function ir_mri_coil_combine_test

if 1 % xtrue
	f.dir = [path_find_dir('mri') '/../data/mri/'];
	f.xtrue = [f.dir 'brainweb_t1.jpg'];
	xtrue = single(imread(f.xtrue)');
	xtrue = xtrue(2:end-1,2:end-1); % make it 256^2
	xtrue = downsample2(xtrue, 4); % now 64^2
	[nx, ny] = size(xtrue);
	atrue = 2*pi * (-0.5+([1:nx]'/nx * [1:ny]/ny).^0.5); % smooth phase
	xtrue = xtrue .* exp(1i * atrue); % phase
%	im clf, im('hsv', angle(xtrue), [-pi pi]), cbar, return

	im plc 3 4
	im(1, abs(xtrue), '$|x|$'), cbar
	im(2, 'hsv', angle(xtrue), [-pi pi], '$\angle x$'), cbar
	mask = conv2(abs(xtrue), ones(5), 'same') > 0;
%	im(11, abs(xtrue) + 100*mask), return % check mask
end

if 1 % smap
	ncoil = 4;
	mask4 = repmat(mask, [1 1 ncoil]);
	smap = ir_mri_sensemap_sim('nx', nx, 'ny', ny, 'dx', 192/nx, ...
		'ncoil', ncoil, 'rcoil', 100);
%	smap = smap(:, :, [ncoil 1:(ncoil-1)]); % phase of last coil is cooler
	tmp = sqrt(sum(abs(smap).^2, 3));
	smap = smap / tmp(end/2,end/2); % normalize true smap for simplicity
	im(5, mask4 .* abs(smap), '$|s|$'), cbar
	im(6, 'hsv', mask4 .* angle(smap), [-pi pi], '$\angle s$'), cbar
	im(12, mask .* sqrt(sum(abs(smap).^2, 3)), 's SSoS'), cbar
end

if 1 % y
	ytrue = smap .* repmat(xtrue, [1 1 4]);
	rng(0)
	snr2sigma = @(db, yb) 10^(-db/20) * norm(yb(:)) / sqrt(numel(yb)) / sqrt(2); % for complex noise
	sig = snr2sigma(50, ytrue);
%	sig = 0; % noiseless
	ykj = ytrue + sig * (randn(size(ytrue)) + 1i * randn(size(ytrue)));
	im(7, abs(ykj), '$|y|$'), cbar
	im(8, 'hsv', ir_unwrap(angle(ykj)), [-pi pi], '$\angle y$'), cbar
end

	[xhat shat ref] = ir_mri_coil_combine(ykj);
	ssos = sqrt(sum(abs(ykj).^2, 3));
	xhat = xhat * ir_best_scale(xhat, abs(xtrue));
	ssos = ssos * ir_best_scale(ssos, abs(xtrue));

	pr nrms(xhat(:), col(abs(xtrue)))
	pr nrms(ssos(:), col(abs(xtrue)))
%	pr corrcoef(abs(xtrue(:)), ssos(:))
%	pr corrcoef(abs(xtrue(:)), real(xhat(:)))

%	xhat = reale(xhat);
%	im(3, abs(xhat), '|x hat|'), cbar
%	im(4, angle(xhat), [-pi pi], '< x hat'), cbar
	im(3, real(xhat), 'Re($\hat{x}$)'), cbar
	im(4, imag(xhat), 'Im($\hat{x}$)'), cbar
	im(9, mask4 .* abs(shat), '$|\hat{s}|$'), cbar
	im(10, 'hsv', mask4 .* angle(shat), [-pi pi], '$\angle \hat{s}$'), cbar
	im(11, ssos, 'x SSoS'), cbar

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


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

clf
[X,Y,Z] = meshgrid(gx,gy,gz);
v = [];
i = 0.05
ii = 0.05
for i3 = 1:2:length(gx)
    i3 = gx(i3);
    l = loc(1,:) == i3;
    re3 = reshape(as(l),length(gz),length(gy));
    r = reshape(loc(4,l),length(gz),length(gy));
    l = r > tresh;
    re3(l) = zeros(1,length(find(l)))*nan;
    rect = [i ii 0.15 0.15];
    i = i+0.20;
    if i > 0.9 ; i = 0.05 ; ii = ii + 0.20; end
    axes('position',rect)
    pcolor(gy,-gz,re3)
    caxis([min(as) max(as)])
    colormap(jet)
    title([ 'Long=' num2str(i3) ' km'],'FontSize',8)
    set(gca,'FontSize',8')
    hold on
    l = a.Longitude > i3-dx/2 & a.Longitude < i3+dx/2;
    plot3(a(l,2),-a(l,7),a(l,7)*0,'k.')
    shading interp
    % axis('off')
    hold on
end
%colorbar

author: CelsoReyes


text: % trainCRFScript

author: zouchuhang

repo: LayoutNet


text:
function [xn,yn] = randvertex(x,y,npix)
%RANDVERTEX Adds random noise to the vertices of a polygon.
%	[XN,YN] = RANDVERTEX[X,Y,NPIX] adds uniformly distributed noise to
%	the coordinates of vertices of a polygon. The coordinates of the
%	vertices are input in X and Y, and NPIX is the maximum number of
%	pixel locations by which any pair (X(i),Y(i)) is allowed to deviate.
%	For example, if NPIX = 1, the location of any X(i) will not deviate
%	by more than one pixel location in the x-direction, and similarly for
%	Y(i). Noise is added independently to the two coordinates.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Convert to columns.
x = x(:);
y = y(:);

% Preliminary calculations.
L = length(x);
xnoise = rand(L,1);
ynoise = rand(L,1);
xdev = npix*xnoise.*sign(xnoise - 0.5);
ydev = npix*ynoise.*sign(ynoise - 0.5);

% Add noise and round.
xn = round(x + xdev);
yn = round(y + ydev);

% All pixel locations must be no less than 1.
xn = max(xn,1);
yn = max(yn,1);

author: dipum

repo: dipum-toolbox

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


text:
function [eigvector, eigvalue, elapse] = PCA2(covdata, options)
%PCA	Principal Component Analysis
%
%	Usage:
%       [eigvector, eigvalue] = PCA(data, options)
%       [eigvector, eigvalue] = PCA(data)
%
%             Input:
%               data       - Data matrix. Each row vector of fea is a data point.
%
%     options.ReducedDim   - The dimensionality of the reduced subspace. If 0,
%                         all the dimensions will be kept.
%                         Default is 0.
%
%             Output:
%               eigvector - Each column is an embedding function, for a new
%                           data point (row vector) x,  y = x*eigvector
%                           will be the embedding result of x.
%               eigvalue  - The sorted eigvalue of PCA eigen-problem.
%
%	Examples:
% 			fea = rand(7,10);
% 			[eigvector,eigvalue] = PCA(fea,4);
%           Y = fea*eigvector;
%
%   version 2.2 --Feb/2009
%   version 2.1 --June/2007
%   version 2.0 --May/2007
%   version 1.1 --Feb/2006
%   version 1.0 --April/2004
%
%   Written by Deng Cai (dengcai2 AT cs.uiuc.edu)
%

if (~exist('options','var'))
   options = [];
end

ReducedDim = size(covdata,1);
if isfield(options,'ReducedDim')
    ReducedDim = options.ReducedDim;
end
tmp_T = cputime;

%[nSmp,nFea] = size(data);

ddata = max(covdata, covdata');

dimMatrix = size(ddata,2);
if dimMatrix > 1000 & ReducedDim < dimMatrix/10  % using eigs to speed up!
    option = struct('disp',0);
    [eigvector, eigvalue] = eigs(ddata,ReducedDim,'la',option);
    eigvalue = diag(eigvalue);
else
    [eigvector, eigvalue] = eig(ddata);
    eigvalue = diag(eigvalue);

    [junk, index] = sort(-eigvalue);
    eigvalue = eigvalue(index);
    eigvector = eigvector(:, index);
end

clear ddata;
maxEigValue = max(abs(eigvalue));
eigIdx = find(abs(eigvalue)/maxEigValue < 1e-12);
eigvalue (eigIdx) = [];
eigvector (:,eigIdx) = [];

if ReducedDim < length(eigvalue)
    eigvalue = eigvalue(1:ReducedDim);
    eigvector = eigvector(:, 1:ReducedDim);
end

if isfield(options,'PCARatio')
    sumEig = sum(eigvalue);
    sumEig = sumEig*options.PCARatio;
    sumNow = 0;
    for idx = 1:length(eigvalue)
        sumNow = sumNow + eigvalue(idx);
        if sumNow >= sumEig
            break;
        end
    end
    eigvector = eigvector(:,1:idx);
end

elapse = cputime - tmp_T;

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function fem = line_swap(fem_struct,j,k);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Line_swap is mesh refining tool that swaps an edge that two elements
% share. It reconnects the elements so that the edge they share has been
% flipped. This routine only operates on one line at a time and can only be
% used if the two elements share an edge. This routine can operate on any
% elemnet, even boundry elements, unless the new elements formed by the line
% swap would overlap existing elements. No nodes are moved by this routine.
% The fem.e and fem.ar are the only fields that are updated. If no .nei file
% is present then none will be generated; however, if one is present then
% it will be updated to reflect the new connectivity.
%
% Calls: is_valid_struct.m
%
% Usage: fem = line_swap(fem_struct,j,k);
%
% Variables:
%  fem -- the new split finite element mesh.
%  fem_struct -- the finite element grid structure from the opnml suite.
%          Note: fem_struct.nei will not be generated if not present
%  j -- element number for one element used in the line swap.
%  k -- element number for one element used in the line swap.
%          Note: j and k can be interchanged.
%
% Filename: line_swap.m
% Created by: Ben Holladay
% Date: June 8, 2005
% Last Modified:
%       Oct. 26, 2007 -- Chris Massey, NRL 7322
%       Changed when line swap can not be performed due to mesh tangling,
%       from an error call to a display message and return original fem.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Ensure the apropriate number of input arguements is given, checks that the
%fem_struct is valid.
if nargin == 0
    error('Not enough input arguments; need a fem_struct and two element numbers.');
end
if ~is_valid_struct(fem_struct)
   error('Input argument to line swap must be a valid fem_struct.');
end
if nargin == 1
    error(['Input arguments must include a valid fem_struct and two element ' , ...
    'numbers that share an edge.']);
elseif nargin == 2
    error('Two element numbers must be given.');
elseif nargin == 3
    if j <= 0 | j > size(fem_struct.e,1) | j ~= floor(j)
        error(['Element numbers must positive integers that are less than ',...
            'the maximum number of elements.']);
    end
    if k <= 0 | k > size(fem_struct.e,1) | k ~= floor(k)
        error(['Element numbers must positive integers that are less than ',...
            'the maximum number of elements.']);
    end
else
    error('Too many input arguments.');
end

%Sets the intial fem_struct variables.
x = fem_struct.x;
y = fem_struct.y;
enodes = fem_struct.e;
ar = fem_struct.ar;

%Determine if the elements share a common edge.
elems = enodes([j;k],:);
comp = unique(intersect(elems(1,:),elems(2,:)));
if length(comp) < 2
    error('The elements do not share a common edge.');
end

%Indentify each node for swapping.
[nodes,ti,tj] = unique(elems);
temp1 = setdiff((1:6),ti);
temp1 = elems(temp1);
temp2 = setdiff(nodes,temp1);

%Determine if the angle between the elements is too big for line_swap to
%effectively operate on it.
a2 = (x(enodes([j;k],3))-x(enodes([j;k],2))).^2+(y(enodes([j;k],3))-y(enodes([j;k],2))).^2;
b2 = (x(enodes([j;k],1))-x(enodes([j;k],3))).^2+(y(enodes([j;k],1))-y(enodes([j;k],3))).^2;
c2 = (x(enodes([j;k],2))-x(enodes([j;k],1))).^2+(y(enodes([j;k],2))-y(enodes([j;k],1))).^2;
A = (180/pi)*acos((b2+c2-a2)./(2*sqrt(b2).*sqrt(c2)));
B = (180/pi)*acos((c2+a2-b2)./(2*sqrt(c2).*sqrt(a2)));
C = (180/pi)*acos((a2+b2-c2)./(2*sqrt(a2).*sqrt(b2)));
ang = [A,B,C];
temp3 = find(temp1(1) == elems);
temp4 = find(temp1(2) == elems);
temp5 = ang(temp3);
temp6 = ang(temp4);
temp5 = sum(temp5);
temp6 = sum(temp6);
if temp5 >= 180 | temp6 >= 180
    disp(['The line cannot be swapped because the new elements would '...
        'create overlapping elements.']);
    fem = fem_struct;
    return
elseif temp5 >= 135 | temp6 >= 135
    disp(['Warning: The elements created will contain large angles.']);
end
clear a2 b2 c2 A B C ang

%Swap the lines in the enodes list.
temp3 = setdiff(elems(1,:),temp1);
temp4 = setdiff(elems(2,:),temp1);
temp5 = find(elems(1,:) == temp1(1));
elems(1,temp5) = temp4;
temp6 = find(elems(2,:) == temp1(2));
elems(2,temp6) = temp3;
clear temp5 temp6;

%Recomputes the areas for the new elements.
xnodes = x(elems);
ynodes = y(elems);
temparea = 0.5*(xnodes(:,1).*(ynodes(:,2)-ynodes(:,3))+xnodes(:,2).*...
(ynodes(:,3)-ynodes(:,1))+xnodes(:,3).*(ynodes(:,1)-ynodes(:,2)));
clear xnodes ynodes x y;

%Addes the new elements and areas to the global lists.
enodes([j;k,],:) = elems;
ar([j;k]) = temparea(:);
clear elems temparea;

%Create the output structure.
fem = fem_struct;
fem.e = enodes;
fem.ar = ar;

%Updates the nei if present.
try
    nei = fem_struct.nei;
    nei = [nei,zeros((size(nei,1)),1)];
    tempnei = nei([temp1(:);temp2(:)],:);
    temp3 = find(tempnei(1,:) == temp1(2));
    temp4 = setdiff(1:(size(tempnei,2)),temp3);
    tempnei(1,:) = ([tempnei(1,temp4),0]);
    temp3 = find(tempnei(2,:) == temp1(1));
    temp4 = setdiff(1:(size(tempnei,2)),temp3);
    tempnei(2,:) = ([tempnei(2,temp4),0]);
    temp3 = find(tempnei(3,:) == temp1(1) | tempnei(3,:) == temp1(2));
    if abs(temp3(1)-temp3(2)) == 1
		tempnei(3,:) = [tempnei(3,(1:temp3(1))),temp2(2),tempnei(3,((temp3(2)):(end-1)))];
    else
        tempnei(3,:) = [temp2(2),tempnei(3,(1:end-1))];
    end
    temp3 = find(tempnei(4,:) == temp1(1) | tempnei(4,:) == temp1(2));
    if abs(temp3(1)-temp3(2)) == 1
        tempnei(4,:) = [tempnei(4,(1:temp3(1))),temp2(1),tempnei(4,((temp3(2)):(end-1)))];
    else
        tempnei(4,:) = [temp2(1),tempnei(4,(1:end-1))];
    end
    nei([temp1(:);temp2(:)],:) = tempnei;
    temp3 = find(nei(:,end) ~= 0);
    if isempty(temp3) == 1
        nei = nei(:,(1:end-1));
    end
    temp3 = find(nei(:,end) ~= 0);
    if isempty(temp3) == 1
        nei = nei(:,(1:end-1));
    end
    fem.nei = nei;
catch
end

return

author: CHLNDDEV


text:
function days = month_length_persian ( y, m )

%*****************************************************************************80
%
%% MONTH_LENGTH_PERSIAN returns the number of days in a Persian month.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year in which the month occurred.
%
%    Input, integer M, the number of the month.
%
%    Output, integer DAYS, the number of days
%    in the month.
%
  mdays = [ 31, 31, 31, 31, 31, 31, 30, 30, 30, 30, 30, 29 ];
%
%  Copy the input.
%
  m2 = m;
  y2 = y;
%
%  Get the number of days in the month.
%
  days = mdays(m2);
%
%  If necessary, add 1 day for a leap year.
%
  if ( m2 == 12 && year_is_leap_persian ( y2 ) )
    days = days + 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function jed = ymd_to_jed_common ( y, m, d )

%*****************************************************************************80
%
%% YMD_TO_JED_COMMON converts a Common YMD date to a JED.
%
%  Discussion:
%
%    The "common" calendar is meant to be the calendar which is Julian up to
%    day JED = 2299160, and Gregorian from day JED = 2299161 and after.
%
%    The Julian Ephemeris Date is essentially a count of the number
%    of days that have elapsed since noon, 1 January 4713 BC, at
%    Greenwich, England.  Strictly speaking, the Julian Ephemeris Date
%    is counted from noon, and thus day "0" began at noon on 1 January 4713 BC,
%    and ended at noon on 2 January 4713 BC.
%
%    The Julian Ephemeris Date was devised by Joseph Scaliger in 1583.
%
%    The Julian Ephemeris Date has been adopted by astronomers as
%    a convenient reference for dates.
%
%  Example:
%
%       Y   M     D         JED
%    --------------     -------
%    BC 4713 Jan  1           0
%    AD 1968 May 23     2440000
%    AD 1984 Dec 31     2446065
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, the YMD date.
%
%    Output, real JED, the Julian Ephemeris Date.
%

%
%  Copy the month and year.
%
  y2 = 1582;
  m2 = 10;
  d2 = 4+1;

  cmp = ymd_compare ( y, m, d, y2, m2, d2 );

  if ( cmp == '<' )
    jed = ymd_to_jed_julian ( y, m, d );
    return
  end
%
%  Use the Gregorian calendar for dates strictly after 1752/9/13.
%
  y2 = 1582;
  m2 = 10;
  d2 = 15-1;

  cmp = ymd_compare ( y, m, d, y2, m2, d2 );

  if ( cmp == '>' )
    jed = ymd_to_jed_gregorian ( y, m, d );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'YMD_TO_JED_COMMON - Fatal error!\n' );
  fprintf ( 1, '  Illegal date!\n' );
  error ( 'YMD_TO_JED_COMMON - Fatal error!' );

  return
end


text:
function H = cross( F, G )
%CROSS   Vector cross product.
%   CROSS(F, G) returns the cross product of the CHEBFUN2V objects F and G. If F
%   and G both have two components, then it returns the CHEBFUN2 representing
%       CROSS(F,G) = F(1) * G(2) - F(2) * G(1)
%   where F = (F(1); F(2)) and G = (G(1); G(2)). If F and G have three
%   components then it returns the CHEBFUN2V representing the 3D cross
%   product.

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

% Empty check:
if ( isempty( F ) || isempty( G ) )
    H = chebfun2v;
    return
end

% Get number of components:
Fc = F.components;
Gc = G.components;

% Do curl:
if ( F.nComponents == 2 && G.nComponents == 2 )  % 2D curl
    H = Fc{1} .* Gc{2} - Fc{2} .* Gc{1};
elseif ( F.nComponents == 3 && G.nComponents == 3 ) % 3D curl
    H = [ Fc{2} .* Gc{3} - Fc{3} .* Gc{2} ; ...
          Fc{3} .* Gc{1} - Fc{1} .* Gc{3} ; ...
          Fc{1} .* Gc{2} - Fc{2} .* Gc{1} ];
else
    error('CHEBFUN:CHEBFUN2V:cross:components', ...
        'CHEBFUN2V objects must be both 2- or 3-vectors.');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% [C, E] = curve_ext(Tx, log2(fs), lambda);
%
% Extract a maximum energy, minimum curvature, curve from
% Synchrosqueezed Representation.  Note, energy is given as:
%   abs(Tx).^2
%
% This implements the solution to Eq. (8) of [1].
%
% 1. E. Brevdo, N.S. Fučkar, G. Thakur, and H-T. Wu, "The
% Synchrosqueezing algorithm: a robust analysis tool for signals
% with time-varying spectrum," 2011.
%
% Original Author: Jianfeng Lu (now at NYU CIMS)
% Maintainer: Eugene Brevdo
%
% Inputs:
%   lambda should be >=0.  Default: lambda=0
%
% Outputs:
%   C: the curve locations (indices)
%   E: the (logarithmic) energy of the curve
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/)
%---------------------------------------------------------------------------------

author: HeLiangHIT


text:
%FASTNLMEANSDENOISINGCOLORED  Modification of fastNlMeansDenoising function for colored images
%
%     dst = cv.fastNlMeansDenoisingColored(src)
%     dst = cv.fastNlMeansDenoisingColored(src, 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ Input 8-bit 3-channel or 4-channel image.
%
% ## Output
% * __dst__ Output image with the same size and type as `src`.
%
% ## Options
% * __H__ Parameter regulating filter strength for luminance component. Bigger
%   `H` value perfectly removes noise but also removes image details, smaller
%   `H` value preserves details but also preserves some noise. default 3
% * __HColor__ The same as `H` but for color components. For most images value
%   equals 10 will be enough to remove colored noise and do not distort
%   colors. default 3
% * __TemplateWindowSize__ Size in pixels of the template patch that is used
%   to compute weights. Should be odd. Recommended value 7 pixels. default 7
% * __SearchWindowSize__ Size in pixels of the window that is used to compute
%   weighted average for given pixel. Should be odd. Affect performance
%   linearly: greater `SearchWindowsSize` - greater denoising time.
%   Recommended value 21 pixels. default 21
% * __FlipChannels__ whether to flip the order of color channels in input
%   `src` and output `dst`, between MATLAB's RGB/RGBA order and OpenCV's
%   BGR/BGRA (input: RGB/RGBA->BGR/BGRA, output: BGR/BGRA->RGB/RGBA).
%   default false
%
% The function converts image to CIELAB colorspace and then separately denoise
% L and AB components with given `H` parameters using cv.fastNlMeansDenoising
% function.
%
% See also: cv.fastNlMeansDenoising, cv.fastNlMeansDenoisingMulti,
%  cv.fastNlMeansDenoisingColoredMulti
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
baseDir = pwd;

[f,p] = uigetfile('*.dcm','Select the first image...',fullfile(baseDir,'raw','t2','I0001.dcm'));
if(isnumeric(f)) error('user cancelled.'); end
baseDir = fileparts(p);
if(isempty(baseDir)) baseDir = pwd; end
[tmpImg,im.mmPerVox] = makeCubeIfiles(fullfile(p,f));
imName = {'pd','t2'};
im.img2std = computeCannonicalXformFromIfile(fullfile(p,f));
[f,p] = uigetfile('*.gz','Select t1 file...',baseDir);
if(isnumeric(f)) error('user cancelled.'); end
refNi = niftiRead(fullfile(p,f));
VG.uint8 = uint8(round(mrAnatHistogramClip(double(refNi.data),0.4,0.98).*255));
VG.mat = refNi.qto_xyz;
for(ii=1:2)
    im.img = tmpImg(:,:,[ii:2:end]);
    % We use header info to get us into axial orientation and ensure the
    % left-right direction is correct.
    [im.img, im.mmPerVox, dimOrder, dimFlip] = applyCannonicalXform(im.img, im.img2std, im.mmPerVox);
    [im.img, im.clipVals] = mrAnatHistogramClip(im.img,0.4,0.98);
    VF.uint8 = uint8(round(im.img.*255));
    im.origin = (size(im.img)+1)./2;
    VF.mat = [[diag(im.mmPerVox), [im.origin.*-im.mmPerVox]']; [0 0 0 1]];

    % Now align it to the T1 to get ac-pc space
    rotTrans = spm_coreg(VF,VG);
    % This composite xform will convert inplane voxel space to inplane
    % physical (VF.mat) and then inplane physical to anat physical (the
    % rigid-body rotTrans that we just computed). Since the anat physical
    % space is ac-pc space, that where we want to be.
    im.xform = spm_matrix(rotTrans(:)')*VF.mat;
    fname = fullfile(baseDir, [imName{ii} '.nii.gz']);
    dtiWriteNiftiWrapper(int16(round(im.img.*diff(im.clipVals))), im.xform, fname);
end
clear dt ip VF VG f;

t2 = niftiRead(fullfile(baseDir, 't2.nii.gz'));
handles = guidata(gcf);
newMm = [0.5 0.5 0.5];
[img,xform] = mrAnatResliceSpm(double(t2.data),t2.qto_ijk,dtiGet(handles,'t1bb'),newMm);
img(isnan(img)) = 0;
img = img./max(img(:));
handles = dtiAddBackgroundImage(handles, img, 't2', newMm, diag([1 1 1 1]));
guidata(gcf, handles);

pd = niftiRead(fullfile(baseDir, 'pd.nii.gz'));
handles = guidata(gcf);
% The nifti q-form goes converts voxel space to ac-pc space, but we want
% voxel to t1-anat, so we undo the ac-pc part.
xform = inv(handles.acpcXform)*pd.qto_xyz;
img = double(pd.data);
img = img./max(img(:));
handles = dtiAddBackgroundImage(handles, img, 'pd', t2.pixdim, xform);
guidata(gcf, handles);

%
% Starting from NIFTIS:
%
baseDir = '/biac3/wandell4/data/Achiasma/DL070825_anatomy/';
ni1 = niftiRead(fullfile(baseDir,'raw','t2pd_1.nii.gz');
ni2 = niftiRead(fullfile(baseDir,'raw','t2pd_2.nii.gz');
ref = niftiRead(fullfile(baseDir,'raw','t1.nii.gz');
im1 = ni1.data(:,:,[1:2:40]);
im2 = ni2.data(:,:,[1:2:40]);
xform =  mrAnatRegister(im2,im1);
%bb = [1 1 1; size(im1)];
%newIm2 =  mrAnatResliceSpm(double(im2), xform, bb, ni1.pixdim);
avgT2 = (double(ni1.data(:,:,[1:2:40])) + double(ni2.data(:,:,[1:2:40])))./2;
avgPd = (double(ni1.data(:,:,[2:2:40])) + double(ni2.data(:,:,[2:2:40])))./2;
avgT2 = int16(round(avgT2));
avgPd = int16(round(avgPd));

VG.uint8 = uint8(round(mrAnatHistogramClip(double(ref.data),0.4,0.98).*255));
VG.mat = ref.qto_xyz;
VF.uint8 = uint8(round(mrAnatHistogramClip(double(avgT2),0.4,0.98).*255));
origin = (size(avgT2)+1)./2;
VF.mat = [[diag(ni1.pixdim), [origin.*-ni1.pixdim]']; [0 0 0 1]];

% Now align it to the T1 to get ac-pc space
spm_defaults; global defaults;
estParams = defaults.coreg.estimate;
estParams.sep = [8 4];
rotTrans = spm_coreg(VF,VG,estParams);
acpcXform = spm_matrix(rotTrans(end,:))*VF.mat;

dtiWriteNiftiWrapper(avgT2, acpcXform, fullfile(baseDir, 't2.nii.gz'));
dtiWriteNiftiWrapper(avgPd, acpcXform, fullfile(baseDir, 'pd.nii.gz'));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function H = PlotSwathe(bar,swathe,SwatheOpt)
% =======================================================================
% Plot a line with a shaded swathe (e.g for impulse responses, etc)
% =======================================================================
% H = PlotSwathe(bar,swathe,SwatheOpt)
% -----------------------------------------------------------------------
% INPUT
%   - bar: line to plot
%   - swathe: if a vector, draws symmetric swathe around bar. Otherwise
%       draws asymmetric swathe, with first columnn being the upper limit
%       and second column being the lower limit
% -----------------------------------------------------------------------
% OPTIONAL INPUT
%   - SwatheOpt.linecol: color of line, can be rgb number or string, e.g.
%       'blue'
%   - SwatheOpt.swathecol: color of swathe (can be rgb number or string,
%       e.g 'blue')
%   - SwatheOpt.transp: 1 for transparency (dflt=0)
%   - nticks: number ticks on the horizontal axis
%   - fo: first observation (convention: 1987.00 = 1987Q1)
%   - SwatheOpt.frequency: quarterly ('q') or annual ('y') [dflt='q']
%   - SwatheOpt.swatheonly: set to 1 to plot only the swathe (dflt=0)
% -----------------------------------------------------------------------
% OUPUT
%   - H.bar: handle for the median
%   - H.patch: handle for the patch
%   - H.swathe: H.swathe(1) handle for upp , H.swathe(2) handle for low
%   - H.ax: hande for the axis
%   - H.nobs: number of observations
% -----------------------------------------------------------------------
% EXAMPLE
% x = 3*(1:50);
% swathe = [1*(1:50);4*(1:50)];
% PlotSwathe(x,swathe)
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2015. Updated November 2020
% -----------------------------------------------------------------------

if ~exist('SwatheOpt','var')
    SwatheOpt = PlotSwatheOption;
end

% Get the color of line
if isnumeric(SwatheOpt.linecol) % in case user provides already the rgb code
    linecol = SwatheOpt.linecol;
else % otherwise use rgb to get the code
    linecol = rgb(SwatheOpt.linecol);
end

% Get the color of swathe
if isnumeric(SwatheOpt.swathecol) % in case user provides already the rgb code
    swathecol = SwatheOpt.swathecol;
else % otherwise use rgb to get the code
    swathecol = rgb(SwatheOpt.swathecol);
end

% Number of observations
nobs = length(bar); H.nobs = nobs;

% Make sure data is row-vector
if size(bar,1)>1
    bar =bar';
end

% Check if swathe is vector and set error bands

if isvector(swathe)
    % Check if swathe has the right dimension [1 x T]
    if size(swathe,1)>1
        swathe = swathe';
        low = bar - swathe;
        upp = bar + swathe;
    else
        low = bar - swathe;
        upp = bar + swathe;
    end
else
    % Check if i has the right dimension [2 x T]
    if size(swathe,1)>2
        swathe = swathe';
        upp = swathe(1,:);
        low = swathe(2,:);
    else
        upp = swathe(1,:);
        low = swathe(2,:);
    end
end

% Create the x axis
xaxis = 1:nobs;

%% PLOT
% Initialize patch
edgeColor = swathecol+(1-swathecol)*0.55;
% Set saturation
swathealpha = SwatheOpt.alpha;
% SwatheOpt.transp to make the patch
if SwatheOpt.transp==1
    faceAlpha = swathealpha;
    patchColor = swathecol;
    set(gcf,'renderer','openGL')
else
    faceAlpha = 1;
    patchColor = swathecol+(1-swathecol)*(1-swathealpha);
    set(gcf,'renderer','painters')
end

% Plot the 'bar' line
if SwatheOpt.swatheonly==0
    H.bar = plot(xaxis,bar,'LineWidth',2,'Color',linecol,'Marker',SwatheOpt.marker);
end

% Add the error-bar plot elements
holdStatus = ishold;
if ~holdStatus, hold on,  end

% Make the coordinates for the patch
yP = [low,fliplr(upp)];
xP = [xaxis,fliplr(xaxis)];

% Remove any nans otherwise patch won't work
xP(isnan(yP))=[];
yP(isnan(yP))=[];
H.patch=patch(xP,yP,1,'facecolor',patchColor,'edgecolor','none','facealpha',faceAlpha);

% Make nice edges around the patch.
H.swathe(1)=plot(xaxis,low,'-','Color',edgeColor);
H.swathe(2)=plot(xaxis,upp,'-','Color',edgeColor);
if SwatheOpt.swatheonly==0
    delete(H.bar)
end
if SwatheOpt.swatheonly==0
    H.bar=plot(xaxis,bar,'LineWidth',2,'Color',linecol,'Marker',SwatheOpt.marker);
end
if ~holdStatus, hold off, end

% Make smart x labels
if SwatheOpt.do_dates==1
    DatesPlot(fo,nobs,nticks,SwatheOpt.frequency)
end
H.ax = gca;

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function Test_fixedrank_2factors()
% function Test_LR()
% Test for fixedrankLRquotientfactory geometry (low rank matrix completion)
%
% Paper link: http://www.icml-2011.org/papers/350_icmlpaper.pdf
%

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, July 11, 2013.
% Contributors:
% Change log:

    clear all; clc; close all;

    % Problem
    m = 20;
    n = 20;
    r = 8;
    A = randn(m, r);
    B = randn(n, r);
    C = A*B';

    % Create the problem structure

%     problem.M = fixedrankfactory_2factors(m, n, r);
    %     problem.M = fixedrankfactory_2factors_preconditioned(m, n, r);
        problem.M = fixedrankfactory_2factors_subspace_projection(m, n, r);

    df = problem.M.dim();
    p = 8*df/(m*n);
    %     mask = rand(m, n) <= p;
    mask = ones(m,n );

    problem.cost = @cost;
    function f = cost(X)
        f = .5*norm(mask.*(X.L*X.R' - C), 'fro')^2;
    end

    problem.grad = @(X) problem.M.egrad2rgrad(X, egrad(X));
    function g = egrad(X)
        P = mask.^2 .* (X.L*X.R' - C);
        g.L = P*X.R;
        g.R = P'*X.L;
    end

    problem.hess = @(X, U) problem.M.ehess2rhess(X, egrad(X), ehess(X, U), U);
    function Ress = ehess(X, eta)
        P = (mask.^2).*( X.L*X.R' - C);
        Pdot = (mask.^2).*(eta.L*X.R' + X.L*eta.R');
        Ress.L = Pdot*X.R + P*eta.R;
        Ress.R = Pdot'*X.L + P'*eta.L;
    end

    % % Check numerically whether gradient and Ressian are correct
    checkgradient(problem);
    drawnow;
    pause;
    checkhessian(problem);
    drawnow;
    pause;

    % Initialization
    X0 = [];

    % Options (not mandatory)
    options.maxiter = inf;
    options.maxinner = 30;
    options.maxtime = 120;
    options.tolgradnorm = 1e-5;
    options.Delta_bar = min(m ,n )*r;
    options.Delta0 = min(m ,n )*r/64;

    % Pick an algorithm to solve the problem
    [Xopt costopt info] = trustregions(problem, X0, options);
    % [Xopt costopt info] = steepestdescent(problem, X0, options);
    % [Xopt costopt info] = conjugategradient(problem, X0, options);

    evs = real(hessianspectrum(problem, Xopt));
    evs = real(evs);
    max(evs)/min(evs)
    stairs(sort(evs));
    title(['Eigenvalues of the Hessian of the cost function ' ...
        'at the solution']);
end

author: NicolasBoumal

repo: manopt


text:
function [X_norm, mu, sigma] = featureNormalize(X)
%FEATURENORMALIZE Normalizes the features in X
%   FEATURENORMALIZE(X) returns a normalized version of X where
%   the mean value of each feature is 0 and the standard deviation
%   is 1. This is often a good preprocessing step to do when
%   working with learning algorithms.

mu = mean(X);
X_norm = bsxfun(@minus, X, mu);

sigma = std(X_norm);
X_norm = bsxfun(@rdivide, X_norm, sigma);

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

end

author: anirudhjayaraman

repo: Machine-Learning


text:
%TEST 07: Simple v.s. Area (two methods of building weight matrix)
clear
clc
fprintf('TEST 07: Simple v.s. Area (two methods of building weight matrix)\r');
im = imtest('s4',156);
figure('name','Original Image: Phantom')
imshow(im)
siz = size(im);
angles = 0:3:179;
[W1,p1,~,~] = buildWeightMatrixSimple(im,angles);
[W2,p2,~,~] = buildWeightMatrixArea(im,angles);
%Simple
fprintf('Reconstruct from simple-W:\r')
im_rec1 = tomo_recon_lsqr(W1, p1, siz, 1e-6, 1000);
im_rec1 = uint8(imscale(im_rec1));
figure('name','Simple')
imshow(im_rec1);
%Area
fprintf('Reconstruct from area-W:\r')
im_rec2 = tomo_recon_lsqr(W2, p2, siz, 1e-6, 1000);
im_rec2 = uint8(imscale(im_rec2));
figure('name','Area')
imshow(im_rec2);

author: Sable

repo: mcbench-benchmarks


text:
function e = year_to_epact_julian ( y )

%*****************************************************************************80
%
%% YEAR_TO_EPACT_JULIAN returns the epact of a Julian year.
%
%  Discussion:
%
%    The epact of a year is the age in days of the notional moon on
%    the first day of the year.  If the year begins with a new moon,
%    the epact is zero.  If the new moon occurred the day before,
%    the epact is 1.  There is a unique epact for every golden number.
%
%    Bear in mind that the notional moon is not the one in the sky,
%    but a theoretical one that satisfactorily approximates the behavior
%    of the real one, but which is tame enough to be described by a formula.
%
%  Example:
%
%    Year  Golden Number  Epact
%
%      1 BC     1           8
%      1 AD     2          19
%      2 AD     3           0
%      3 AD     4          11
%      4 AD     5          22
%      5 AD     6           3
%      6 AD     7          14
%      7 AD     8          25
%      8 AD     9           6
%      9 AD    10          17
%     10 AD    11          28
%     11 AD    12           9
%     12 AD    13          20
%     13 AD    14           1
%     14 AD    15          12
%     15 AD    16          23
%     16 AD    17           4
%     17 AD    18          15
%     18 AD    19          26
%     19 AD     1           8
%     20 AD     2          19
%   1066 AD     3           0
%   1900 AD     1           8
%   1919 AD     1           8
%   1938 AD     1           8
%   1957 AD     1           8
%   1976 AD     1           8
%   1995 AD     1           8
%   2014 AD     1           8
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999.
%
%  Parameters:
%
%    Input, integer Y, the year.  The year 0 is illegal input.
%
%    Output, integer E, the epact, between 0 and 28.
%
  if ( y == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'YEAR_TO_EPACT_JULIAN - Fatal error!\n' );
    fprintf ( 1, '  Illegal input Y = 0.\n' );
    error ( 'YEAR_TO_EPACT_JULIAN - Fatal error!' );
  end

  g = year_to_golden_number ( y );

  e = i4_wrap ( 11 * g - 3, 0, 29 );

  return
end


text:
function T = auxstructure(elem)
%% AUXSTRUCTURE auxiliary structure for a 2-D triangulation.
%
%  T = AUXSTRUCTURE(elem) constucts the indices map between elements, edges
%  and nodes, and the boundary information. T is a structure.
%
%  T.neighbor(1:NT,1:3): the indices map of neighbor information of elements,
%  where neighbor(t,i) is the global index of the element oppoiste to the
%  i-th vertex of the t-th element.
%
%  T.elem2edge(1:NT,1:3): the indices map from elements to edges, elem2edge(t,i)
%  is the edge opposite to the i-th vertex of the t-th element.
%
%  T.edge(1:NE,1:2): all edges, where edge(e,i) is the global index of the
%  i-th vertex of the e-th edge, and edge(e,1) < edge(e,2)
%
%  T.bdEdge(1:Nbd,1:2): boundary edges with positive oritentation, where
%  bdEdge(e,i) is the global index of the i-th vertex of the e-th edge for
%  i=1,2. The positive oritentation means that the interior of the domain
%  is on the left moving from bdEdge(e,1) to bdEdge(e,2). Note that this
%  requires elem is positive ordered, i.e., the signed area of each
%  triangle is positive. If not, use elem = fixorder(node,elem) to fix the
%  order.
%
%  T.edge2elem(1:NE,1:4): the indices map from edge to element, where
%  edge2elem(e,1:2) are the global indexes of two elements sharing the e-th
%  edge, and edge2elem(e,3:4) are the local indices of e to edge2elem(e,1:2).
%
%  To save space all the data type in T is uint32. When use them as a input
%  of sparse(i,j,s,m,n), please change them into double type.
%
%  See also auxstructure3.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

totalEdge = uint32(sort([elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])],2));
matlabversion = version;
if str2double(matlabversion(end-5:end-2)) > 2012
    [edge, i2, j] = unique(totalEdge,'rows','legacy');
else
    [edge, i2, j] = unique(totalEdge,'rows');
end
NT = size(elem,1);
elem2edge = uint32(reshape(j,NT,3));
i1(j(3*NT:-1:1)) = 3*NT:-1:1;
i1 = i1';
k1 = ceil(i1/NT);
k2 = ceil(i2/NT);
t1 = i1 - NT*(k1-1);
t2 = i2 - NT*(k2-1);
ix = (i1 ~= i2);
neighbor = uint32(accumarray([[t1(ix),k1(ix)];[t2,k2]],[t2(ix);t1],[NT 3]));
edge2elem = uint32([t1,t2,k1,k2]);
bdElem = t1(t1 == t2);
bdk1 = k1(t1 == t2);
bdEdge = [elem(bdElem(bdk1==1),[2 3]); elem(bdElem(bdk1==2),[3 1]);...
          elem(bdElem(bdk1==3),[1 2])];
T = struct('neighbor',neighbor,'elem2edge',elem2edge,'edge',edge,...
           'edge2elem',edge2elem,'bdElem',bdElem,'bdEdge',bdEdge);

repo: Aero-matlab


text:
function [subj_descriptions, event_descriptions] = list_variables(D, varargin)
% Prints a list of variable names and descriptions for every variable, or specified variables
%
% :Usage:
% ::
%
%    print_summary(D, [optional inputs])
%
% ..
%     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
%
% :Optional Inputs:
%
%   **subj:**
%        followed by a cell array of subject level var names, to only see those vars
%
%   **event:**
%        followed by a cell array of event level var names, to only see those vars
%
%   if either varargin is unspecified, all variables will be printed
%

subj_descriptions = table;
event_descriptions = table;

subj_varnames = D.Subj_Level.names;
if iscolumn(subj_varnames), subj_varnames = subj_varnames'; end % need row

svars = find(strcmp('subj', varargin));
if ~isempty(svars), subj_varnames = varargin{svars+1}; end

% Names and descriptions - Subject level
% --------------------------------------------------------------------
Var_name = subj_varnames';

clear Description

for i=1:length(subj_varnames)

    vname = subj_varnames{i};
    [var,~,~,Description{i, 1}] = get_var(D, vname);

end

if length(subj_varnames) > 0
    subj_descriptions = table(Var_name, Description);
end

fprintf('\n\n --------- SUBJECT LEVEL VARS: data_obj.Subj_Level -------- \n\n');
disp(subj_descriptions);

% Names and descriptions - Event level
% --------------------------------------------------------------------
event_varnames = D.Event_Level.names;
if iscolumn(event_varnames), event_varnames = event_varnames'; end % need row

evars = find(strcmp('event', varargin));
if ~isempty(evars), event_varnames = varargin{evars+1}; end

Var_name = event_varnames';

clear Description

for i=1:length(event_varnames)

    vname = event_varnames{i};
    [var,~,~,Description{i, 1}] = get_var(D, vname);

end

if length(event_varnames) > 0
    event_descriptions = table(Var_name, Description);
end

fprintf('\n\n --------- EVENT LEVEL VARS: data_obj.Event_Level -------- \n\n');
disp(event_descriptions);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%%  Generate a set of data points from a plane
%Copyright (c) 2011, The MathWorks, Inc.

clear all
clc

[CleanX, CleanY] = meshgrid(1:10);
CleanX = reshape(CleanX,100,1);
CleanY = reshape(CleanY,100,1);
CleanZ = 3*CleanX + 4*CleanY;

scatter3(CleanX, CleanY, CleanZ, 'filled')
hold on
foo = fit([CleanX, CleanY], CleanZ, 'poly11')
plot(foo)

%% Add noise vectors to all three dimensions

% Note:  Normally, if we were going to show a regression example, we'd
% create Noisy_Z by adding a noise vector to CleanZ but leave CleanX and
% Clean Y as is.  Here, we're using a fitting technique that is designed to
% create a model where there is noise associated with both the dependent
% and the independent variables.

NoisyX = CleanX + randn(100,1);
NoisyY = CleanY + randn(100,1);
NoisyZ = CleanZ + randn(100,1);

% create a scatter plot

figure1 = figure;
scatter3(NoisyX,NoisyY,NoisyZ, '.k')

%%  Fit a plane to the data using OLS
[foo2, GoF] = fit([NoisyX NoisyY],NoisyZ, 'poly11');

% superimpose the plane on the scatter plot
hold on
h1 = plot(foo2)
set( h1, 'FaceColor', 'r' )

%% Calculate the sum of the Squared Errors

% Calculate the difference between the observed (NoisyX and NoisyY) and
% the "known" value for CleanZ

resid1 = CleanZ - foo(NoisyX, NoisyY);

% Square residuals
resid1_sqrd = resid1.^2;

% Take the sum of the squared residuals
SSE_OLS = sum(resid1_sqrd)

% Create textbox
annotation(figure1,'textbox',...
    [0.147443514644351 0.802739726027397 0.305903765690377 0.0931506849315069],...
    'String',{['SSE OLS = ' num2str(SSE_OLS)]},...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 1]);

%% Use Principal Component Analysis to perform an Orthogonal Regression

% PCA is based on centering and rotation.
% PCA rotates the data such that dimension with the greatest amount of
% variance is parallel to the X axis.  The dimension with the second
% largest amount of variance will be parallel to the Y axis.  This operation
% defines a plane.  The direction with the third largest variance will be
% parallel to the Z axis.  This dimension defines a set of residuals which
% are at right angles to the XY plane.

[coeff,score,roots] = princomp([NoisyX NoisyY NoisyZ]);

basis = coeff(:,1:2);
normal = coeff(:,3);
pctExplained = roots' ./ sum(roots)

% Translate the output from PCA back to the original coordinate system

[n,p] = size([NoisyX NoisyY NoisyZ]);
meanNoisy = mean([NoisyX NoisyY NoisyZ],1);
Predicted = repmat(meanNoisy,n,1) + score(:,1:2)*coeff(:,1:2)';

% Generate a fit object that represents the output from the Orthogonal Regression

[foo3, Gof2] = fit([Predicted(:,1) Predicted(:,2)], Predicted(:,3), 'poly11')

h2  = plot(foo3)
set( h2, 'FaceColor', 'b' )

%% Calculate the Sum of the Squared Errors for the Orthogonal Regression
% Calculate residuals
resid2 = CleanZ - Predicted(:,3);

% Square residuals
resid2_sqrd = resid2.^2;

% Take the sum of the squared residuals
SSE_TLS = sum(resid2_sqrd)

annotation(figure1,'textbox',...
    [0.147443514644351 0.802739726027397 0.305903765690377 0.0931506849315069],...
    'String',{['SSE OLS = ' num2str(SSE_OLS)], ['SSE TLS = ' num2str(SSE_TLS)]},...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 1]);


text:
function mpc = t_case3p_c
%T_CASE3P_C   Six bus hybrid test case, 2 single-phase buses, 4 3-phase buses.
%
%   One bus is a hybrid PV bus (PV on single-phase side). Three phase bus
%   solution should match T_CASE3P_A

%% MATPOWER Case Format : Version 2
mpc.version = '2';

%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 100;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [
	1	3	0	0	0	0	1	1	0.20115203	12.47	1	1.1	0.9;
% 	2	2	0	0	0	0	1	1	0	12.47	1	1.1	0.9;
% 	3	2	0	0	0	0	1	1	0	12.47	1	1.1	0.9;
	4	2	0	0	0	0	1	1	0	12.47	1	1.1	0.9;
% 	5	1	90	30	0	0	1	1	0	12.47	1	1.1	0.9;
% 	6	1	0	0	0	0	1	1	0	12.47	1	1.1	0.9;
% 	7	1	100	35	0	0	1	1	0	12.47	1	1.1	0.9;
% 	8	1	0	0	0	0	1	1	0	12.47	1	1.1	0.9;
% 	9	1	125	50	0	0	1	1	0	12.47	1	1.1	0.9;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
mpc.gen = [
	1	72.3	27.03	300	-300	1.0024291384	100	1	250	10	0	0	0	0	0	0	0	0	0	0	0;
	4	0	0	300	-300	1	100	1	250	10	0	0	0	0	0	0	0	0	0	0	0;
% 	2	163	6.54	300	-300	1.025	100	1	300	10	0	0	0	0	0	0	0	0	0	0	0;
% 	3	85	-10.95	300	-300	1.025	100	1	270	10	0	0	0	0	0	0	0	0	0	0	0;
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [
	1	4	0	0.0576	0	250	250	250	0	0	1	-360	360;
% 	4	5	0.017	0.092	0.158	250	250	250	0	0	1	-360	360;
% 	5	6	0.039	0.17	0.358	150	150	150	0	0	1	-360	360;
% 	3	6	0	0.0586	0	300	300	300	0	0	1	-360	360;
% 	6	7	0.0119	0.1008	0.209	150	150	150	0	0	1	-360	360;
% 	7	8	0.0085	0.072	0.149	250	250	250	0	0	1	-360	360;
% 	8	2	0	0.0625	0	250	250	250	0	0	1	-360	360;
% 	8	9	0.032	0.161	0.306	250	250	250	0	0	1	-360	360;
% 	9	4	0.01	0.085	0.176	250	250	250	0	0	1	-360	360;
];

%%-----  OPF Data  -----%%
%% generator cost data
%	1	startup	shutdown	n	x1	y1	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
	2	1500	0	3	0.11	5	150;
	2	1500	0	3	0.11	5	150;
% 	2	2000	0	3	0.085	1.2	600;
% 	2	3000	0	3	0.1225	1	335;
];

%%-----  3 Phase Model Data  -----%%
%% system data
mpc.freq = 60;      %% frequency, Hz
mpc.basekVA = 1000; %% system kVA base

%% bus data
%	busid	type	basekV	Vm1	Vm2	Vm3	Va1	Va2	Va3
mpc.bus3p = [
	1	1	12.47	1	1	1	0	-120	120;
	2	1	12.47	1	1	1	0	-120	120;
	3	1	4.16	1	1	1	0	-120	120;
	4	1	4.16	1	1	1	0	-120	120;
];

%% buslink data
%	linkid	busid	bus3pid	status
mpc.buslink = [
	1	4	1	1;
];

%% branch data
%	brid	fbus	tbus	status	lcid	len
mpc.line3p = [
	1	1	2	1	1	2000/5280;
	2	3	4	1	1	2500/5280;
];

%% transformer
%	xfid	fbus	tbus	status	R	X	basekVA	basekV
mpc.xfmr3p = [
	1	2	3	1	0.01	0.06	6000	12.47;
];

%% load
%	ldid	ldbus	status	Pd1	Pd2	Pd3	ldpf1	ldpf2	ldpf3
mpc.load3p = [
	1	4	1	1275	1800	2375	0.85	0.9	0.95;
];

%% gen
%	genid	gbus	status	Vg1	Vg2	Vg3	Pg1	Pg2	Pg3	Qg1	Qg2	Qg3
mpc.gen3p = [
% 	1	1	1	1	1	1	2000	2000	2000	0	0	0;
];

%% line construction
%	lcid	R11	R21	R31	R22	R32	R33	X11	X21	X31	X22	X32	X33	C11	C21	C31	C22	C32	C33
mpc.lc = [
	1	0.457541	0.15594 	0.153474	0.466617	0.157996	0.461462	1.078	0.501648	0.384909	1.04813	0.423624	1.06502	15.0671	-4.86241	-1.85323	15.875	-3.09098	14.3254
];

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
filename='Cantilever_hexahedra_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'MMA';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.05;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
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 plotPredictionResults_VASARI(pathResults,nameOutcome,fSetNames,metric,maxOrder,pathFig)
% -------------------------------------------------------------------------
% function plotPredictionResults_VASARI(pathResults,nameOutcome,fSetNames,metrics,maxOrder)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function plots prediction performance estimation results for all the
% different feature set types entered as inputs in the LGG study.
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathResults: Full path to the 'RESULTS' folder where prediction results
%                 are saved.
%                 --> Ex: '/myProject/WORKSPACE/VASARI/RESULTS'
% 2. nameOutcome: String specifying the name of the outcome being displayed
%                 --> Ex: 'progression'
% 3. fSetNames: Cell of strings specifying the name of the type of feature
%               set analyzed.
%               --> Ex: {'VASARI'}
% 4. metric: String specifying the metric to display.
%            --> 'AUC632'
% 5. maxOrder: Integer specifying the maximal multivariable model order.
%              --> Ex: 10
% 6. pathFig: (optional).  Full path to where figure is saved without
%             displaying it. Put '' for displaying the figure and not
%             saving it to 'pathFig' (default).
%             --> Ex: ''
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2017
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2017  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

if nargin < 6
    pathFig = '';
end

startpath = pwd;
cd(pathResults)

signs = {'-r',':b','--g','-.y'};
nFSET = numel(fSetNames);

maxOrderChosen = maxOrder;

if isempty(pathFig)
    figure
else
    h = figure('visible','off');
end
for i = 1:nFSET
    fSET = fSetNames{i};
    results = load(['RESULTS_',fSET,'_',nameOutcome]); results = struct2cell(results); results = results{1};
    nOrders = numel(fieldnames(results));
    if nOrders < maxOrderChosen
        maxOrder = nOrders;
    else
        maxOrder = maxOrderChosen;
    end
    val = zeros(maxOrder,1);
    val_SE = zeros(maxOrder,1);
    for j = 1:maxOrder
        orderName = ['Order',num2str(j)];
        val(j,1) = results.(orderName).(metric);
        val_SE(j,1) = results.(orderName).(['SE_',metric]);
    end
    errorbar(1:maxOrder,val(:,1),val_SE(:,1),signs{i},'LineWidth',3,'MarkerFaceColor',signs{i}(end),'MarkerSize',6)
    hold on
end
set(gca,'FontSize',20)
xlabel('Model Order','FontSize',24)
ylabel('Prediction performance','FontSize',24)
ind = strfind(metric,'632');
if ~isempty(ind)
    metric = [metric,'+'];
end
metric = [metric(1:ind-1),'_{',metric(ind:end),'}'];
ind = strfind(nameOutcome,'Death');
if ~isempty(ind)
    nameOutcome(ind:ind+4) = [];
    nameOutcome = ['Survival',nameOutcome];
end
titleName = [nameOutcome,'(VASARI) -- ',metric];
title(titleName,'FontSize',30,'FontWeight','bold')
legend(fSetNames,'Location','SouthEast')
axis([0 maxOrderChosen+1 0.5 1])
set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10])

if ~isempty(pathFig)
    cd(pathFig)
    saveas(h,titleName,'fig')
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [objfea, roomfea, unaryanglesid] = compObjHypsFeature( rotImg, omap, gc, cn, CAN_POOL, typenum, bufName)
%COMPOBJHYPSFEATURE Summary of this function goes here
%   Detailed explanation goes here
% if ~exist('unaryname','var')
%     unaryname = './object_classifier/objectunary.mat';
% end
global config;

all_hyps = CAN_POOL.sel_hyps;
room_obj = CAN_POOL.room;

num_object = length(all_hyps.objects);
objfea = repmat(struct('rectscore',[],'segmscore',[], ...
    'randomforest',[],'unaryscore',[], ...
    'indicator',[],'om_sum',[],'om_avg', [], 'om_std', [], ...
    'gc_sum',[],'gc_avg', [], 'gc_std', [], ...
    'cn_hist', [], 'cn_entropy', [], 'sz',[]), ...
    num_object, 1);
room_xyz = [min(room_obj.out_points_w,[],1) max(room_obj.out_points_w,[],1)];

%% get image score
checkrect = true(1,1);
checkseg = true(1,1);
objects = all_hyps.objects;
[ cuboid_rect_score, segment_rect_score] = ...
    computeCuboidImgFea( rotImg, {[objects(1); objects]}, checkrect, checkseg, bufName );

for oid = 1:num_object
    objfea(oid).rectscore = max(cuboid_rect_score{1}(oid+1).surfscore);
    if isempty(objfea(oid).rectscore)
        objfea(oid).rectscore = -10;
    end
    objfea(oid).segmscore = segment_rect_score{1}(oid+1);
end

%% get semantic score
% load('./object_classifier/objectunary.mat');
load(config.unaryName);
[ unaryscore, unaryanglesid ] = unarySizeScore( room_xyz, all_hyps.obj_xyz, typenum, unary_size);
forcecuboid = [2 7 9 10 12];
forcerect = [1 3 4 8];
isrectangle = any((all_hyps.obj_xyz(:,4:6)-all_hyps.obj_xyz(:,1:3))<0.01,2);
iscuboid = ~isrectangle;
unaryscore(isrectangle, forcecuboid) = 1000000;
unaryscore(iscuboid, forcerect) = 1000000;

for oid = 1:num_object
    objfea(oid).unaryscore = unaryscore(oid,:);
    objfea(oid).randomforest = all_hyps.scores(oid,:);
end

roomfea.unaryscore = unarySizeScore( room_xyz, [], 0, unary_size);
roomfea.omapScr = room_obj.omapScr;
roomfea.gcScr = room_obj.gcScr;
roomfea.mgScr = room_obj.mgScr;
roomfea.opScr = room_obj.opScr;

%% get vectorized marker
coor = [];
% load('./rectangleDetector/segmentation/uniformvector_lvl6.mat');
[ coor, tri ] = getUniformVector( 6 );
[ indicator ] = insideIndicator( all_hyps.objects, coor );
for oid = 1:num_object
    objfea(oid).indicator = indicator(:,oid);
end

%% get omap and gc feature
[omH, omW, C] = size(omap);
if C~=3
    fprintf('Warning: OMAP should have 3 channels.\n');
end
coor2D = uv2coords(xyz2uvN(coor), omW, omH);
coor2D_indi = sub2ind([omH omW], coor2D(:,2), coor2D(:,1));
[ om_sum, om_avg, om_std ] = omapFeature( omap, indicator, coor2D_indi );
[ gc_sum, gc_avg, gc_std  ] = gcFeature( gc, indicator, coor2D_indi );
[ cn_hist, cn_entropy ] = colorNameFeature( cn, indicator, coor2D_indi );
sz_vec = sum(indicator,1)';
for oid = 1:num_object
    objfea(oid).om_sum = om_sum(oid,:);
    objfea(oid).om_avg = om_avg(oid,:);
    objfea(oid).om_std = om_std(oid,:);

    objfea(oid).gc_sum = gc_sum(oid,:);
    objfea(oid).gc_avg = gc_avg(oid,:);
    objfea(oid).gc_std = gc_std(oid,:);

    objfea(oid).cn_hist = cn_hist(oid,:);
    objfea(oid).cn_entropy = cn_entropy(oid);
    objfea(oid).sz = sz_vec(oid);
end

%%

end

author: zouchuhang


text:
function test_example_bids_emg

% MEM 4gb
% WALLTIME 00:20:00

%% Converting an example EMG dataset for sharing in BIDS
%
% The [BIDS standard](https://bids.neuroimaging.io) does currently not specify how to represent EMG recordings. This example - and the support that is implemented in the **[data2bids](https://github.com/fieldtrip/fieldtrip/blob/release/data2bids.m)** function - should be considered as a preliminary proposal to help researchers with their existing data. This example may also serve to start a discussion on whether and how this data type should be added to the [BIDS specification](http://bids-specification.readthedocs.io/).
%
% EMG is often recorded along with other experiments, such as MEG or functional MRI recordings, but can also be recorded by itself. EMG can for example be used to check compliance with the task, to detect movements and compute reaction times, or to check the activation pattern of specific muscles in a task.
%
% EMG is not part as an individual data type in the BIDS specification. In the **[data2bids](https://github.com/fieldtrip/fieldtrip/blob/release/data2bids.m)** function we are mostly treating it similar to EEG data, with only some differences in the metadata. E.g., the |EEGPlacementScheme|, |CapManufacturer|, and |CapManufacturersModelName| do not apply.  Other than that, we have implemented it the same, with the same data formats and the same way to code events.
%
%% # Example
%
% The example data is available from our [FTP server](ftp://ftp.fieldtriptoolbox.org/pub/fieldtrip/example/bids_emg/)
%
% The example includes the original log files under the “original” directory. It also has a copy of the script to do the conversion under “code”. The reorganized data is under the “bids” directory. According to [the documentation](https://bids-specification.readthedocs.io/en/stable/02-common-principles.html#source-vs-raw-vs-derived-data) the original data can be added to the BIDS dataset under the “sourcedata” directory and code can be added to the “code” directory. This way no information is lost and the conversion/reorganization is fully reproducible.
%
% In general: if your original data is in a different format than the BIDS representation (e.g., DICOMs instead of NIfTI, or data from an EMG system in an unsupported format), you might want to keep the sourcedata, e.g., on a data acquisition collection on the Donders Repository. If it is in the same format (as here), you could simply delete the copy in sourcedata. Your own analyses and the sharing of the data would be based on the BIDS representation.
%
cfg = [];

cfg.InstitutionName             = 'Radboud University';
cfg.InstitutionalDepartmentName = 'Donders Institute for Brain, Cognition and Behaviour';
cfg.InstitutionAddress          = 'Kapittelweg 29, 6525 EN, Nijmegen, The Netherlands';

% required for dataset_description.json
cfg.dataset_description.Name                = 'EMG example';
cfg.dataset_description.BIDSVersion         = 'unofficial extension';

% optional for dataset_description.json
cfg.dataset_description.License             = 'n/a';
cfg.dataset_description.Authors             = 'n/a';
cfg.dataset_description.Acknowledgements    = 'n/a';
cfg.dataset_description.Funding             = 'n/a';
cfg.dataset_description.ReferencesAndLinks  = 'n/a';
cfg.dataset_description.DatasetDOI          = 'n/a';

cfg.method = 'copy'; % the original data is in a BIDS -compliant format and can simply be copied

cfg.bidsroot = fullfile(tempdir, 'bids');  % write to the some scratch directory
cfg.datatype = 'emg';
cfg.sub = 'NT05';

% these are general fields
cfg.Manufacturer           = 'BrainProducts';
cfg.ManufacturersModelName = 'BrainAmp MR plus';

% these are EMG specific fields
cfg.emg.EMGPlacementScheme = 'electrode pairs were placed over various muscles';
cfg.emg.EMGReference       = 'bipolar';
cfg.emg.EMGGround          = 'unclear where the ground electrode was placed on the body';
cfg.emg.EMGChannelCount    = 13; % it might actually be that some channels are recording something else

% there are two recordings of the EMG, during the RS scan and during the GABA scan
% in both EMG recordings the task for the subject was the same

cfg.TaskDescription = 'The subject was lying in the MRI scanner during a resting-state BOLD scan, while EMG was being recorded';
cfg.task = 'resting';
cfg.dataset = dccnpath('/home/common/matlab/fieldtrip/data/ftp/example/bids_emg/original/nt05_EMG_SESS1_RS_RAW.vhdr');
cfg.run = 1;

% write the first EMG recording as run-1
data2bids(cfg);

cfg.TaskDescription = 'The subject was lying in the MRI scanner during a GABA MRS scan, while EMG was being recorded';
cfg.task = 'resting';
cfg.dataset = dccnpath('/home/common/matlab/fieldtrip/data/ftp/example/bids_emg/original/nt05_EMG_SESS1_GABA_RAW.vhdr');
cfg.run = 2;

% write the first EMG recording as run-2
data2bids(cfg);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function value = hartley_determinant ( n )

%*****************************************************************************80
%
%% HARTLEY_DETERMINANT returns the determinant of the HARTLEY matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real VALUE, the determinant.
%
  if ( mod ( n, 4 ) == 1 )
    value =   sqrt ( n ^ n );
  else
    value = - sqrt ( n ^ n );
  end

  return
end


text:
function y = lt( x, y )

if isa(x,'cvxin')||~isa(y,'cvxin')||y.active,
    error( 'CVX error: improper use of the <in> pseudo-operator.' );
end
y.active = true;
y.value  = x;

% 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 cx = cheby_u_poly ( m, n, x )

%*****************************************************************************80
%
%% CHEBY_U_POLY evaluates Chebyshev polynomials U(n,x).
%
%  Differential equation:
%
%    (1-X*X) Y'' - 3 X Y' + N (N+2) Y = 0
%
%  First terms:
%
%    U(0,X) =   1
%    U(1,X) =   2 X
%    U(2,X) =   4 X^2 -   1
%    U(3,X) =   8 X^3 -   4 X
%    U(4,X) =  16 X^4 -  12 X^2 +  1
%    U(5,X) =  32 X^5 -  32 X^3 +  6 X
%    U(6,X) =  64 X^6 -  80 X^4 + 24 X^2 - 1
%    U(7,X) = 128 X^7 - 192 X^5 + 80 X^3 - 8X
%
%  Recursion:
%
%    U(0,X) = 1,
%    U(1,X) = 2 * X,
%    U(N,X) = 2 * X * U(N-1,X) - U(N-2,X)
%
%  Norm:
%
%    Integral ( -1 <= X <= 1 ) ( 1 - X^2 ) * U(N,X)^2 dX = PI/2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of evaluation points.
%
%    Input, integer N, the highest polynomial to compute.
%
%    Input, real X(M,1), the evaluation points.
%
%    Output, real V(M,N+1), the values of the Chebyshev polynomials
%    0 through N at X(1:M).
%
  if ( n < 0 )
    v = [];
    return
  end

  v = zeros ( m, n + 1 );

  v(1:m,1) = 1.0;

  if ( n < 1 )
    return
  end

  x = x(:);

  v(1:m,2) = 2.0 * x(1:m,1);

  for j = 2 : n
    v(1:m,j+1) = 2.0 * x(1:m,1) .* v(1:m,j) - v(1:m,j-1);
  end

  return
end


text:
function test_bug1984

% WALLTIME 00:10:00
% MEM 3gb
% DEPENDENCY

dataset = dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.ds');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% raw data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cfg = [];
cfg.dataset                 = dataset;
cfg.trialdef.eventtype      = 'backpanel trigger';
cfg.trialdef.prestim        = 1;
cfg.trialdef.poststim       = 2;
cfg.trialdef.eventvalue     = 3;                    % trigger value for fully incongruent (FIC)
cfg = ft_definetrial(cfg);
dataFIC = ft_preprocessing(cfg);

cfg = [];
cfg.dataset                 = dataset;
cfg.trialdef.eventtype      = 'backpanel trigger';
cfg.trialdef.prestim        = 1;
cfg.trialdef.poststim       = 1;                    % make these trials shorter
cfg.trialdef.eventvalue     = 9;                    % trigger value for fully congruent (FC)
cfg = ft_definetrial(cfg);
dataFC = ft_preprocessing(cfg);

cfg = [];
cfg.dataset = dataset;
cfg.trl = [
  dataFIC.cfg.trl
  dataFC.cfg.trl
  ];
[dum, order] = sort(cfg.trl(:,1));
cfg.trl = cfg.trl(order,:); % keep them in the original order of the experiment
dataAll = ft_preprocessing(cfg);

cfg = [];
dataAppend = ft_appenddata(cfg, dataFIC, dataFC);

%% compare the two versions

% the order should be different
assert(~isequal(dataAll.trialinfo,  dataAppend.trialinfo));
assert(~isequal(dataAll.sampleinfo, dataAppend.sampleinfo));

% sort them according to the trial order in the original experiment
[dum, order1] = sort(dataAll.sampleinfo(:,1));
[dum, order2] = sort(dataAppend.sampleinfo(:,1));

% after sorting it should be the same
assert(isequal(dataAll.trialinfo(order1,:),  dataAppend.trialinfo(order2,:)));
assert(isequal(dataAll.sampleinfo(order1,:), dataAppend.sampleinfo(order2,:)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% freq data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cfg = [];
cfg.keeptrials = 'yes';
cfg.method = 'mtmfft';
cfg.foi = 1:10;
cfg.tapsmofrq = 2*(cfg.foi);

freqFIC = ft_freqanalysis(cfg, dataFIC);
freqFC  = ft_freqanalysis(cfg, dataFC);
freqAll = ft_freqanalysis(cfg, dataAll);

cfg = [];
freqAppend = ft_appendfreq(cfg, freqFIC, freqFC);

%% compare the two versions

% these should not have sampleinfo as per specification in FT_DATATYPE_FREQ
assert(~isfield(freqAppend, 'sampleinfo'));
assert(~isfield(freqAll,    'sampleinfo'));

% the order should be different
assert(~isequal(freqAll.trialinfo,  freqAppend.trialinfo));
assert(~isequal(freqAll.cumsumcnt,  freqAppend.cumsumcnt));
assert(~isequal(freqAll.cumtapcnt,  freqAppend.cumtapcnt));

% after sorting it should be the same
assert(isequal(freqAll.trialinfo(order1,:),  freqAppend.trialinfo(order2,:)));
assert(isequal(freqAll.cumsumcnt(order1,:),  freqAppend.cumsumcnt(order2,:)));
assert(isequal(freqAll.cumtapcnt(order1,:),  freqAppend.cumtapcnt(order2,:)));


text:
function [kg] = slug2kg(slug)
% Convert units of mass from slugs to kilograms.
% Chad A. Greene 2012
kg = slug*14.5939;


text:
clear;
clc;
close all;

genJSON('LEEDS');

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
function M = spm_eeg_morlet(Rtf, ST, f, ff)
% Generate Morlet wavelets
% FORMAT M = spm_eeg_morlet(Rtf, ST, f, ff)
%
% Rtf - 'wavelet factor', see [1]
% ST  - sample time [ms]
% f   - vector of frequencies [Hz]
% ff  - frequency to fix Gaussian envelope (sigma = Rtf/(2*pi*ff))
%       Default is ff = f, ie.e, a Morlet transform
%       NB: FWHM = sqrt(8*log(2))*sigma_t;
%
% M   - cell vector, where each element contains the filter for each
%       frequency in f
%__________________________________________________________________________
%
% spm_eeg_morlet generates morlet wavelets for specified frequencies f with
% a specified ratio Rtf, see [1], for sample time ST (ms). One obtains the
% wavelet coefficients by convolution of a data vector with the kernels in
% M. See spm_eeg_tf how one obtains instantaneous power and phase estimates
% from the wavelet coefficients.
%
% [1] C. Tallon-Baudry, O. Bertrand, F. Peronnet and J. Pernier, 1998.
% Induced gamma-Band Activity during the Delay of a Visual Short-term
% memory Task in Humans. The Journal of Neuroscience (18): 4244-4254.
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: spm_eeg_morlet.m 7122 2017-06-22 14:54:01Z guillaume $

if nargin < 4
    ff = f;
else
    ff = repmat(ff,1,numel(f));
end

M = cell(1,numel(f));

for i = 1:numel(f)

    % fixed or scale-dependent window
    %----------------------------------------------------------------------
    sigma_t = Rtf/(2*pi*ff(i));

    % this scaling factor is proportional to (Tallon-Baudry, 1998):
    % (sigma_t*sqrt(pi))^(-1/2);
    %----------------------------------------------------------------------
    t = 0:ST*0.001:5*sigma_t;
    t = [-t(end:-1:2) t];
    M{i} = exp(-t.^2/(2*sigma_t^2)) .* exp(2 * 1i * pi * f(i) *t);
    M{i} = M{i} ./ (sqrt(0.5*sum(real(M{i}).^2 + imag(M{i}).^2)));
    M{i} = M{i} - mean(M{i});

end

author: spm


text:
function out = numColumns(f)
%NUMCOLUMNS   Number of columns (or rows) of a CHEBFUN quasimatrix.
%   NUMCOLUMNS(F) returns the number of columns of a column CHEBFUN, or the
%   number of rows of a row CHEBFUN. It is equivalent to min(size(f)).
%
% See also SIZE, NUMEL.

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

if ( isempty(f) )
    % An empty CHENFUN has no columns:
    out = 0;
elseif ( numel(f) == 1 )
    % Possible array-valued:
    out = size(f.funs{1}, 2);
else
    % A quasimatrix:
    out = numel(f);
end

end


text:
function y = perform_huffcoding(x,T,dir)

% perform_huffcoding - perform huffman coding
%
%   y = perform_huffcoding(x,T,dir);
%
%   dir=+1 for coding
%   dir=-1 for decoding
%
%   T is a Huffman tree, computed with compute_hufftree
%
%   Copyright (c) 2008 Gabriel Peyre

if dir==1
    %%% CODING %%%
    C = huffman_gencode(T);
    m = length(C);
    x = round(x);
    if min(x)<1 || max(x)>m
        error('Too small or too large token');
    end
    y = [];
    for i=1:length(x)
        y = [y C{x(i)}];
    end
else
    %%% DE-CODING %%%
    t = T{1};
    y = [];
    while not(isempty(x))
        if x(1)==0
            t = t{1};
        else
            t = t{2};
        end
        x(1) = [];
        if not(iscell(t))
            y = [y t];
            t = T{1};
        end
    end
    y = y(:);
end

%%

function C = huffman_gencode(T)

if not(iscell(T)) % || not(length(T)==2)
    C = {};
    C{T} = -1;
elseif length(T)==1
    C = huffman_gencode(T{1});
    % remove traling -1
    for i=1:length(C)
        C{i} = C{i}(1:end-1);
    end
elseif length(T)==2
    C1 = huffman_gencode(T{1});
    C2 = huffman_gencode(T{2});
    C = {};
    for i=1:length(C1)
        if not(isempty(C1{i}))
            C{i} = [0 C1{i}];
        end
    end
    for i=1:length(C2)
        if not(isempty(C2{i}))
            C{i} = [1 C2{i}];
        end
    end
else
    error('Problem');
end

author: gpeyre

repo: matlab-toolboxes


text:
function pass = test_iszero()
% Test spherefun iszero() command.

f = spherefun( @(x,y,z) 0 + 0*x );
pass(1) = iszero( f );

f = spherefun( @(x,y,z) cos(x) );
pass(2) = ~iszero( f );

end

author: chebfun


text:
function [uu, erriter, num, tt] = test_CMF3D
%
%   Function test_CMF3D
%
%   The matlab function to show how to use the functions CMF3D_Mex and CMF3D_GPU
%
%   Before using the functions CMF3D_mex, you should compile it as follows:
%       >> mex CMF3D_mex.c
%
%   Before using the functions CMF3D_GPU, you should compile the GPU program:
%       >> nvmex -f nvmexopts.bat CMF3D_GPU.cu -IC:\cuda\v4.0\include -LC:\cuda\v4.0\lib\x64 -lcufft -lcudart
%
%   After compilation, you can define all the parameters (penalty, C_s, C_t, para) as follows:
%
%        - penalty: point to the edge-weight penalty parameters to
%                   total-variation function.
%
%          For the case without incorporating image-edge weights,
%          penalty is given by the constant everywhere. For the case
%          with image-edge weights, penalty is given by the pixelwise
%          weight function:
%
%          for example, penalty(x) = b/(1 + a*| grad f(x)|) where b,a > 0 .
%
%        - C_s: point to the capacities of source flows ps
%
%        - C_t: point to the capacities of sink flows pt
%
%        - para: a sequence of parameters for the algorithm
%             para[0,1,2]: rows, cols, heights of the given image
%             para[3]: the maximum iteration number
%             para[4]: the error bound for convergence
%             para[5]: cc for the step-size of augmented Lagrangian method
%             para[6]: the step-size for the graident-projection step to the
%                    total-variation function. Its optimal range is [0.1, 0.17].
%
%
%       Example:
%
%             >> [u, erriter, i, timet] = CMF3D_GPU(single(penalty), single(Cs), single(Ct), single(para));
%
%             >> us = max(u, beta);  % where beta in (0,1)
%
%             >> figure, loglog(erriter,'DisplayName','erriterN');figure(gcf)
%
%             >> isosurface(u,0.5), axis([1 rows 1 cols 1 heights]), daspect([1 1 1]);
%
%
%   Please email Jing Yuan (cn.yuanjing@gmail.com) for any questions,
%   suggestions and bug reports
%
%   The Software is provided "as is", without warranty of any kind.
%
%               Version 1.0
%   https://sites.google.com/site/wwwjingyuan/
%
%   Copyright 2011 Jing Yuan (cn.yuanjing@gmail.com)
%

nfi_data = load_nii('IM_0190_frame_01.nii');
ur = double(nfi_data.img)/255;
%ur = ur(1:100,1:100,1:100);

[rows,cols,heights] = size(ur);

varParas = [rows; cols; heights; 200; 5e-4; 0.35; 0.11];
%                para 0,1,2 - rows, cols, heights of the given image
%                para 3 - the maximum number of iterations
%                para 4 - the error bound for convergence
%                para 5 - cc for the step-size of augmented Lagrangian method
%                para 6 - the step-size for the graident-projection of p

penalty = 0.2*ones(rows, cols, heights);

ulab(1) = 0.2;
ulab(2) = 0.7;

% build up the priori L_2 data terms
fCs = abs(ur - ulab(1));
fCt = abs(ur - ulab(2));

% ----------------------------------------------------------------------
%  Use the function CMF3D_mex to run the algorithm on CPU
% ----------------------------------------------------------------------

%[uu, erriter,num,tt] = CMF3D_mex(single(penalty), single(fCs), single(fCt), single(varParas));

% ----------------------------------------------------------------------
%  Use the function CMF3D_GPU to run the algorithm on GPU
% ----------------------------------------------------------------------

[uu, erriter,num,tt] = CMF3D_GPU(single(penalty), single(fCs), single(fCt), single(varParas));


text:
function setImageRotate(view, imageRotate)
%
% setImageRotate(view, imageRotate)
%
% Sets imageRotate slider values.
if isequal(view.name, 'hidden')
	% instead of a handle to the slider, the ImageRotate field will
	% directly store the rotations for hidden views:
	view.ui.ImageRotate = imageRotate;
else
	% non-hidden view: view.ui.ImageRotate should be a slider struct
	setSlider(view,view.ui.ImageRotate, imageRotate);
end

return


text:
%% Copyright (C) 2016 Lagu
%% Copyright (C) 2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym divisors (@var{x})
%% Get divisors of integer.
%%
%% Example:
%% @example
%% @group
%% x = sym(150);
%% y = divisors(x)
%%   @result{} y = (sym) [1  2  3  5  6  10  15  25  30  50  75  150]  (1×12 matrix)
%% @end group
%% @end example
%% @end defmethod

%% Reference: http://docs.sympy.org/dev/modules/ntheory.html

function y = divisors(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = pycall_sympy__ ('return S(divisors(_ins[0])),', x);
  y = cell2sym(y);
end

%!test
%! assert( isequal( divisors(sym(150)), divisors(sym(-150)) ))

repo: octsympy


text:
function guv = evaluate_log_posterior_grad(this, uv)
%EVALUATE_LOG_POSTERIOR computes the gradient of the log-posterior
%   (negative energy) wrt the flow fields UV
%   Actually only proportional to the log posterior since the variance of neither the
%   spatial nor the data terms is considered
%
%   This is a member function of the class 'hs_optical_flow'.
%
%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-11-30 $
%   $Revision: $
%
% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

% Spatial term
S       = this.spatial_filters;
gu1     = zeros(size(uv,1), size(uv,2));
gv1     = gu1;

for i = 1:length(S)

    u_ = conv2(uv(:,:,1), S{i}, 'valid');
    v_ = conv2(uv(:,:,2), S{i}, 'valid');

    Si = reshape(S{i}(end:-1:1), size(S{i}));

    if isa(this.rho_spatial_u{i}, 'robust_function')
        u_ = -reshape(deriv(this.rho_spatial_u{i}, u_(:)), size(u_));
        v_ = -reshape(deriv(this.rho_spatial_v{i}, v_(:)), size(v_));
    elseif isa(this.rho_spatial_u{i}, 'gsm_density')
        u_ = reshape(evaluate_log_grad(this.rho_spatial_u{i}, u_(:)'), size(u_));
        v_ = reshape(evaluate_log_grad(this.rho_spatial_v{i}, v_(:)'), size(v_));
    else
        error('evaluate_log_posterior: unknown rho function!');
    end;

    gu1 = gu1+conv2(u_, Si, 'full');
    gv1 = gv1+conv2(v_, Si, 'full');
end;

gu2     = zeros(size(uv,1), size(uv,2));
gv2     = gu2;

% Data term
[It Ix Iy] = partial_deriv(this.images, uv, this.interpolation_method);

if isa(this.rho_data, 'robust_function')
    temp   = -reshape(deriv(this.rho_data, It(:)), size(It));

elseif isa(this.rho_data, 'gsm_density')
    temp   = reshape(evaluate_log_grad(this.rho_data, It(:)'), size(It));
else
    error('evaluate_log_posterior: unknown rho function!');
end;

gu2     = temp.*Ix;
gv2     = temp.*Iy;

guv = cat(3, gu2+this.lambda*gu1, gv2+this.lambda*gv1);

author: kristinbranson

repo: JAABA


text:
function [eegPath] = eeg_toolbox_path

% eeg_toolbox_path - locate the eeg_toolbox installation path
%
% [eegPath] = eeg_toolbox_path
%
% This function looks for the installation path of
% eeg_toolbox.m in the matlab path.  If not
% found, update the matlab path (try addpath).
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:53 $

% Licence:  GNU GPL, no implied or express warranties
% History:  12/2003, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

eegPath = fileparts(which('eeg_toolbox'));
if isempty(eegPath),
    msg = sprintf('Cannot find eeg_toolbox on the matlab path.\nPlease use the addpath command.\n\n');
    error(msg);
else
    eegPath = strcat(eegPath,filesep);
end

return

author: PatternRecognition

repo: OpenBMI


text:
function p = mstud_lpdf(x, mu, Cov, nu)
% p = mstud_lpdf(x, mu, Cov, nu)

D = length(x);

x = x(:);
mu = mu(:);

invCov = inv(Cov);

p = gammaln((D+nu)/2);
p = p - gammaln(nu/2);
p = p - D/2 * log(nu * pi);
p = p - 0.5 * log(det(Cov));
p = p + -(D+nu)/2 * log(1 + 1/nu*(x-mu)'*invCov*(x-mu));

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function pdf = semicircular_pdf ( x, a, b )

%*****************************************************************************80
%
%% SEMICIRCULAR_PDF evaluates the Semicircular PDF.
%
%  Formula:
%
%    PDF(X)(A,B) = ( 2 / ( B * PI ) ) * SQRT ( 1 - ( ( X - A ) / B )**2 )
%    for A - B <= X <= A + B
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x < a - b )

    pdf = 0.0;

  elseif ( x <= a + b )

    y = ( x - a ) / b;

    pdf = 2.0 / ( b * pi ) * sqrt ( 1.0 - y * y );

  elseif ( a + b < x )

    pdf = 0.0;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Z,W,P,T] = dosc(X,Y,nocomp,tol)
%  DOSC Direct Orthogonal Signal Correction
%
%  [Z,W,P,T] = dosc(X,Y,nocomp,tol);
%
%  Input
%     X       Matrix of predictor variables (usually spectra)  (I x J)
%     Y       Predicted variables (e.g. concentration)         (I x K)
%     nocomp  number of DOSC components to calculate
%     tol     tolerance used to calculate pseudinverse of X
%             a tolerance of 1E-3 worked well in two cases presented in paper [1]
%
%
%  Output
%     Z       DOSC corrected matrix X                          (I x J)
%     W       Weights used to determine DOSC components        (J x nocomp)
%     P       Loadings used to remove DOSC component from X    (J x nocomp)
%     T       DOSC components                                  (I x nocomp)
%
%
%  Once the calibration is done, new (scaled) x-data can be corrected by
%  newx = x - x*W*P'; Or use mfile dosc_pred
%
% See Reference
% Westerhuis JA, de Jong S and Smilde AK, Direct orthogonal signal correction,
% Chemometrics and Intelligent Laboratory Systems, 56, (2001), 13-25.

% Sijmen de Jong, Oct 99
% Adjusted by Johan Westerhuis
%   ==========================================================================
%   Copyright 2005 Biosystems Data Analysis Group ; Universiteit van Amsterdam
%   ==========================================================================

% project Y onto X      (step 1)
Yhat = X*(pinv(X')'*Y);

% deflate X wrt Yhat    (step 2)
AyX = X-Yhat*(pinv(Yhat)*X);

% find major PCs of AyX
[Ta,D] = eigs(AyX*AyX',nocomp);

% Calculate pseudoinverse of X using tolerance (step 5a)
pinvX = pinv(X',tol)';

% The tolerance can be changed to the number of PCR components (nc) used
% to estimate the pseudo inverse.
% [U,S,V] = svd(X,0);
% Xcorr = U(:,1:nc)*S(1:nc,1:nc)*V(:,1:nc)';
% pinvX = pinv(Xcorr);

W = pinvX*Ta;
T = X*W;

% Calculate loadings to remove DOSC component (step 7a)
P = X'*T*inv(T'*T);

% deflate X wrt to DOSC components (step 6a)
Z = X - T*P';

author: FuSiry

repo: Spectral-preprocessing-algorithm