text:
% Move the roots of a polynom inside the unit circle
%
% Description
%  Stabilize a polynom.
%  Compared to polystab, polystab2 adds 2 elements
%   * Use the LEJA's order to recover the stable polynom from its roots in order
%     to minimize the numerical errors
%   * Even though the result should be stable in a single iteration, numerical
%     errors can still prevent the stability. Therefore, the stabilization step
%     is runned until the final result is finally stable.
%
% Inputs
%  a        : The AR coefficients to stabilize
%  [modmax] : Maximum modulus of the roots.
%             Any root with a bigger modulus than this value will be moved
%             towards the origin.
%             (def. 1)
%
% Outputs
%
% Copyright (c) 2011 University of Crete - Computer Science Department
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function b = polystab2(a, modmax)

    if nargin<2; modmax=1; end

    b = a;

    v = roots(a);

    idx = find(abs(v)>modmax);
    while ~isempty(idx)

        v(idx) = modmax*exp(1j*angle(v(idx)));

        ind = find(a~=0);
        b = a(ind(1))*poly(leja(v)); % Use LEJA's order to minimize num errors

        % Return only real coefficients if input was real:
        if ~any(imag(a))
            b = real(b);
        end

        a = b;
        v = roots(a);

        idx = find(abs(v)>modmax);
    end

    b = b(:);
return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function ndx = read_txt(infilename)
% Creates an Ndx object from information stored in a text file.
% Inputs:
%   infilename: The name of the text file to read.
% Outputs:
%   ndx: An object of type Ndx constructed from the information in
%     the file.

assert(nargin==1)
assert(isa(infilename,'char'))

fid = fopen(infilename);
lines = textscan(fid,'%s%s');
fclose(fid);

models = lines{1};
testsegs = lines{2};
modelset = {};
segset = {};
modelset = union(modelset,models);
segset = union(segset,testsegs);

trialmask = false(length(modelset),length(segset));
[dummy,ii] = ismember(models,modelset);
[dummy,jj] = ismember(testsegs,segset);
trialmask(sub2ind(size(trialmask),ii,jj)) = true;

ndx = Ndx();
ndx.modelset = modelset;
ndx.segset = segset;
ndx.trialmask = trialmask;
assert(ndx.validate())

end

author: nesl


text:
function [F] = spm_vb_regionF (Y,xY,SPM)
% Get log model evidence over a region of data for a GLM
% FORMAT [F] = spm_vb_regionF (Y,xY,SPM)
%
% Y     Matrix of fMRI data (eg. from spm_summarise.m)
% xY    Coordinates etc from region (eg. from spm_voi.m)
% SPM   SPM data structure (this must be loaded in from an
%       SPM.mat file). If this field is not specified this function
%       will prompt you for the name of an SPM.mat file
%
% F     Log model evidence (single number for whole region)
%
% Importantly, the design matrix is normalised so that when you compare
% models their regressors will be identically scaled.
%
% Valid model comparisons also require that the DCT basis set used in high
% pass filtering, as specified in SPM.xX.K.X0, is the same for all models
% that are to be compared.
%
% W. Penny, G. Flandin, and N. Trujillo-Barreto. (2007). Bayesian Model
% Comparison of Spatially Regularised General Linear Models. Human
% Brain Mapping, 28(4):275-293.
%__________________________________________________________________________
% Copyright (C) 2005-2014 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_vb_regionF.m 6079 2014-06-30 18:25:37Z spm $

try
    SPM;
catch
    [Pf, sts] = spm_select(1,'^SPM\.mat$','Select SPM.mat');
    if ~sts, return; end
    swd = spm_file(Pf,'fpath');
    load(fullfile(swd,'SPM.mat'))
end

s=length(SPM.Sess);
if s>1
    disp('spm_vb_regionF: only works for single session data');
    return
end

% Normalise design matrices to have columns of unit
% variance - except for last
X=SPM.xX.X;
X=norm_design(X);

% Make residual forming matrix (to later remove drifts)
X0=SPM.xX.K.X0;
N=size(X0,1);
iX0=pinv(X0);
R=eye(N)-X0*iX0;

% Remove drifts
Y=R*Y;

% Normalise Y to unit dev
sigma=std(Y);
N=size(Y,1);
Y=Y./(ones(N,1)*sigma);

xyz  = xY.XYZmm;
N    = size(xyz,2);
xyz  = SPM.xVol.M\[xyz;ones(1,size(xyz,2))];
xyz  = xyz(1:3,:);
vxyz = spm_vb_neighbors(xyz',1);

priors.W='Voxel - Shrinkage';
priors.A='Voxel - Shrinkage';

block=spm_vb_init_volume(X,1);
block=spm_vb_set_priors(block,priors,vxyz);
block.maxits = 64;
block.update_alpha  = 1;
block.update_w      = 1;
block.update_lambda = 1;
block.update_F      = 1;
block=spm_vb_glmar(Y,block);

F=block.F;

function [Xout] = norm_design (Xin)
% Normalise design matrix to have unit variance columns
% FORMAT [Xout] = norm_design (Xin)

X=Xin(:,1:end-1);
s=std(X,[],1);
iS=diag(1./s);
Xout=[X*iS,Xin(:,end)];

author: spm


text:
function [Q] = CouetteIC2D(x, y, time)

% function [Q] = CouetteIC2D(x, y, time)
% Purpose: evaluat solution for Couette flow

gamma = 1.4;

rad = sqrt(x.^2 + y.^2); theta = atan2(y, x);

Q = zeros(size(x,1), size(x,2), 4);

utheta = (-rad + 16./rad)/75;
p = 1 + (1./(75^2))*( (rad.^2)/2 - 32*log(rad) - 128./rad.^2 );

Q(:,:,1) = 1;
Q(:,:,2) = -sin(theta).*utheta;
Q(:,:,3) =  cos(theta).*utheta;
Q(:,:,4) = p/(gamma-1) + 0.5*(Q(:,:,2).^2 + Q(:,:,3).^2)./Q(:,:,1);
return

repo: Aero-matlab


text:
initialize;
configjson.imageLocation = fullfile(pwd, 'demo_img');
if(~exist(fullfile(pwd, 'demo_result')))
    mkdir(fullfile(pwd, 'demo_result'));
end
configjson.outputLocation = fullfile(pwd, 'demo_result');

imageFile='demo_img/2007_009084.jpg';
objectProposal='randomPrim';

%%% Proposals for all images in imageLocation

runObjectProposals(objectProposal);

% Proposals for one image
% 1)

proposals=runObjectProposals(objectProposal,imageFile);

% 2)

im=imread(imageFile);
proposals=runObjectProposals(objectProposal,im);

%%% Certain number of proposals

proposals=runObjectProposals(objectProposal,im,100);

author: Cloud-CV

repo: object-proposals


text:
function [Ei,Angle] = Associate(PopObj,W)
% Associate each solution with a reference direction

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    Angle      = acos(1-pdist2(PopObj,W,'cosine'));
    [Angle,Ei] = min(Angle',[],1);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function varargout = out_figure_image( hFig, imgFile, imgLegend)
% OUT_FIGURE_IMAGE: Save window contents as a bitmap image.
%
% USAGE: img = out_figure_image(hFig)           : Extract figure image and return it
%              out_figure_image(hFig, imgFile)  : Extract figure image and save it to imgFile
%              out_figure_image(hFig, 'Viewer') : Extract figure image and open it with the image viewer
%              out_figure_image(hFig, figFile)  : Save figure as a Matlab .fig file after removing some callbacks
%              out_figure_image(hFig, 'Figure') : Copy figure and removes callbacks
%              out_figure_image(hFig)           : Extract figure image and save it to a user selected file
%              out_figure_image(hFig, ..., imgLegend)
%              out_figure_image(hFig, ..., 'time')

% @=============================================================================
% 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, 2008-2019

global GlobalData;
drawnow;

% No legend: plot the time
if (nargin < 3)
    imgLegend = 'time';
elseif isempty(imgLegend)
    imgLegend = '';
end

%% ===== GET FILENAME =====
% If image filename is not specified
if (nargin <= 1) && (nargout == 0)
    % === Build a default filename ===
    % Get default directories and format
    LastUsedDirs = bst_get('LastUsedDirs');
    DefaultFormats = bst_get('DefaultFormats');
    if isempty(DefaultFormats.ImageOut)
        DefaultFormats.ImageOut = 'tif';
    end
    % Get the default filename (from the window title)
    wndTitle = get(hFig, 'Name');
    if isempty(wndTitle)
        imgDefautFile = 'img_default';
    else
        imgDefautFile = file_standardize(wndTitle, 0, '_', 1);
        imgDefautFile = strrep(imgDefautFile, '__', '_');
    end
    % Add extension
    imgDefautFile = [imgDefautFile, '.', lower(DefaultFormats.ImageOut)];
    imgDefautFile = strrep(imgDefautFile, '_.', '.');

    % === Ask user filename ===
    % Ask confirmation for the figure filename
    imgDefaultFile = bst_fullfile(LastUsedDirs.ExportImage, imgDefautFile);
    [imgFile, FileFormat] = java_getfile('save', 'Save figure as...', imgDefaultFile, 'single', 'files', ...
        {{'.tif'}, 'TIFF image, compressed (*.tif)',      'TIF'; ...
         {'.jpg'}, 'JPEG image (*.jpg)',                  'JPG'; ...
         {'.bmp'}, 'Bitmap file (*.bmp)',                 'BMP'; ...
         {'.fig'}, 'Matlab figure (*.fig)',               'FIG'; ...
         {'.png'}, 'Portable Network Graphics (*.png)',   'PNG'; ...
         {'.hdf'}, 'Hierarchical Data Format (*.hdf)',    'HDF'; ...
         {'.pbm'}, 'Portable bitmap (*.pbm)',             'PBM'; ...
         {'.pgm'}, 'Portable Graymap (*.pgm)',            'PGM'; ...
         {'.ppm'}, 'Portable Pixmap (*.ppm)',             'PPM';}, DefaultFormats.ImageOut);
    if isempty(imgFile)
        return
    end
    % Save new default export path
    LastUsedDirs.ExportImage = bst_fileparts(imgFile);
    bst_set('LastUsedDirs', LastUsedDirs);
    % Save default export format
    DefaultFormats.ImageOut = FileFormat;
    bst_set('DefaultFormats',  DefaultFormats);
elseif (nargout == 0)
    [fPath,fBase,fExt] = bst_fileparts(imgFile);
    if ~isempty(fExt)
        FileFormat = upper(fExt(2:end));
    else
        FileFormat = '';
    end
else
    FileFormat = 'ARRAY';
    imgFile = '';
end

%% ===== PREPARE FIGURE =====
% Get figure Type
[tmp,iFig,iDS] = bst_figures('GetFigure', hFig);
if ~isempty(iDS)
    FigureId = GlobalData.DataSet(iDS).Figure(iFig).Id;
else
    FigureId = [];
end
% If figure is a registered data figure
if ~isempty(iDS) && ~isempty(imgLegend)
    % If 3DAxes => Add a time legend
    if strcmpi(imgLegend, 'time')
        isAddTime = ~isempty(FigureId) && (strcmpi(FigureId.Type, '3DViz') || strcmpi(FigureId.Type, 'Topography'));
        if isAddTime && ~isempty(GlobalData.UserTimeWindow.CurrentTime)
            if (GlobalData.UserTimeWindow.CurrentTime > 2)
                imgLegend = sprintf('%4.3fs ', GlobalData.UserTimeWindow.CurrentTime);
            else
                imgLegend = sprintf('%dms ', round(GlobalData.UserTimeWindow.CurrentTime * 1000));
            end
        else
            imgLegend = '';
        end
    end
    % Create legend
    if ~isempty(imgLegend)
        % Get figure color
        figColor = get(hFig, 'Color');
        figPos   = get(hFig, 'Position');
        % Find opposite colors
        if (sum(figColor .^ 2) > 0.8)
            textColor = [0 0 0];
        else
            textColor = [.3 1 .3];
        end
        % Create axes
        hAxesLabel = axes(...
            'Parent',        hFig, ...
            'Tag',           'AxesTimestamp', ...
            'Units',         'Pixels', ...
            'Position',      [0, 0, figPos(3), 30], ...
            'Color',         'none', ...
            'XColor',        figColor, ...
            'YColor',        figColor, ...
            'XLim',          [0,1], ...
            'YLim',          [0,1], ...
            'YGrid',         'off', ...
            'XGrid',         'off', ...
            'XMinorGrid',    'off', ...
            'XTick',         [], ...
            'YTick',         [], ...
            'TickLength',    [0,0], ...
            'XTickLabel',    [], ...
            'Box',           'off', ...
            'Interruptible', 'off', ...
            'BusyAction',    'queue');
        % Create label
        hLabel = text(...
            10 / figPos(3), .5, imgLegend, ...
            'FontUnits',   'points', ...
            'FontWeight',  'bold', ...
            'FontSize',    15, ...
            'Color',       textColor, ...
            'Interpreter', 'none', ...
            'Tag',         'LabelTimestamp', ...
            'Parent',      hAxesLabel);
    end
else
    imgLegend = '';
end
% For time series figures: hide buttons
if ~isempty(FigureId) && ismember(FigureId.Type, {'DataTimeSeries', 'ResultsTimeSeries', 'Spectrum'})
    % Find existing buttons
    if bst_get('isJavacomponent')
        hButtons = findobj(hFig, 'Type', 'hgjavacomponent');
    else
        hButtons = findobj(hFig, 'Type', 'uicontrol');
        if (length(hButtons) > 1)
            hButtons(cellfun(@(c)isempty(strfind(c, 'Button')), get(hButtons, 'Tag'))) = [];
        elseif (length(hButtons) == 1) && isempty(strfind(get(hButtons, 'Tag'), 'Button'))
            hButtons = [];
        end
    end
    isVisible = get(hButtons, 'Visible');
    % Hide them
    set(hButtons, 'Visible', 'off');
else
    hButtons = [];
end
% Focus on figure (captures the contents the topmost figure)
pause(.01);
drawnow;
figure(hFig);
drawnow;

%% ===== SAVE FIGURE =====
if (~isempty(FileFormat) && strcmpi(FileFormat, 'FIG')) || strcmpi(imgFile, 'Figure')
    % === PREPARE FIGURE ===
    % Check figure type
    if ~isempty(FigureId) && ((strcmpi(FigureId.Type, 'Topography') && strcmpi(FigureId.SubType, '2DLayout')) || strcmpi(FigureId.Type, 'Connect') || strcmpi(FigureId.Type, 'MriViewer'))
        bst_error('This figure cannot be exported as a Matlab .fig file.', 'Save figure', 0);
        return;
    end
    % Ask use to remove callbacks or not
    if strcmpi(imgFile, 'Figure')
        isRemoveCallbacks = 0;
    else
        isRemoveCallbacks = java_dialog('confirm', ...
            ['If you save the figure with all the callback functions, you will' 10 ...
             'keep the interactivity of the figure but you will need to have ' 10 ...
             'the same version of Brainstorm running to open it again.' 10 10 ...
             'Remove callbacks from the figure?' 10], 'Save figure');
    end

    % Save existing callbacks
    CloseRequestFcn_bak       = get(hFig, 'CloseRequestFcn');
    KeyPressFcn_bak           = get(hFig, 'KeyPressFcn');
    WindowButtonDownFcn_bak   = get(hFig, 'WindowButtonDownFcn');
    WindowButtonMotionFcn_bak = get(hFig, 'WindowButtonMotionFcn');
    WindowButtonUpFcn_bak     = get(hFig, 'WindowButtonUpFcn');
    ResizeFcn_bak             = get(hFig, 'ResizeFcn');
    if isprop(hFig, 'KeyReleaseFcn')
        KeyReleaseFcn_bak = get(hFig, 'KeyReleaseFcn');
    end
    if isprop(hFig, 'WindowScrollWheelFcn')
        WindowScrollWheelFcn_bak = get(hFig, 'WindowScrollWheelFcn');
    end

    % Remove callbacks that crash without Brainstorm
    set(hFig, 'CloseRequestFcn', 'closereq');
    set(hFig, 'KeyPressFcn', []);
    if isprop(hFig, 'KeyReleaseFcn')
        set(hFig, 'KeyReleaseFcn', []);
    end
    % Remove other callbacks
    if isRemoveCallbacks
        set(hFig, 'WindowButtonDownFcn', []);
        set(hFig, 'WindowButtonMotionFcn', []);
        set(hFig, 'WindowButtonUpFcn', []);
        set(hFig, 'ResizeFcn', []);
        if isprop(hFig, 'WindowScrollWheelFcn')
            set(hFig, 'WindowScrollWheelFcn', []);
        end
    end
    % Find colorbar
    hColorbar = findobj(hFig, '-depth', 1, 'Tag', 'Colorbar');
    if ~isempty(hColorbar)
        ColorbarButton_bak = get(hColorbar, 'ButtonDownFcn');
        set(hColorbar, 'ButtonDownFcn', []);
    end
    % Display figure toolbar
    set(hFig, 'MenuBar', 'figure', 'ToolBar', 'auto');

    % Copy figure
    if strcmpi(imgFile, 'Figure')
        % Copy graphic objects
        hNewFig = copyobj(hFig, 0);
        % Copy callbacks
        set(hNewFig, 'CloseRequestFcn',       'closereq');
        set(hNewFig, 'KeyPressFcn',           []);
        set(hNewFig, 'WindowButtonDownFcn',   WindowButtonDownFcn_bak);
        set(hNewFig, 'WindowButtonMotionFcn', WindowButtonMotionFcn_bak);
        set(hNewFig, 'WindowButtonUpFcn',     WindowButtonUpFcn_bak);
        set(hNewFig, 'ResizeFcn',             ResizeFcn_bak);
        if isprop(hNewFig, 'KeyReleaseFcn')
            set(hNewFig, 'KeyReleaseFcn', []);
        end
        if isprop(hNewFig, 'WindowScrollWheelFcn')
            set(hNewFig, 'WindowScrollWheelFcn', WindowScrollWheelFcn_bak);
        end
        % Copy appdata
        AppData = getappdata(hFig);
        fields = fieldnames(AppData);
        for i = 1:length(fields)
            setappdata(hNewFig, fields{i}, AppData.(fields{i}));
        end
        % Reposition figure
        set(hNewFig, 'Position', get(hNewFig, 'Position') - [0, 70, 0, 0]);
    % Save figure
    else
        saveas(hFig, imgFile, 'fig');
    end

    % === RESTORE FIGURE ===
    % Restore common callbacks
    set(hFig, 'CloseRequestFcn', CloseRequestFcn_bak);
    set(hFig, 'KeyPressFcn',     KeyPressFcn_bak);
    if isprop(hFig, 'KeyReleaseFcn')
        set(hFig, 'KeyReleaseFcn', KeyReleaseFcn_bak);
    end
    % Restore colobar
    if ~isempty(hColorbar)
        set(hColorbar, 'ButtonDownFcn', ColorbarButton_bak);
    end
    % Restore optional callbaks
    if isRemoveCallbacks
        set(hFig, 'WindowButtonDownFcn',   WindowButtonDownFcn_bak);
        set(hFig, 'WindowButtonMotionFcn', WindowButtonMotionFcn_bak);
        set(hFig, 'WindowButtonUpFcn',     WindowButtonUpFcn_bak);
        set(hFig, 'ResizeFcn',             ResizeFcn_bak);
        if isprop(hFig, 'WindowScrollWheelFcn')
            set(hFig, 'WindowScrollWheelFcn',  WindowScrollWheelFcn_bak);
        end
    end
    % Remove figure toolbar
    set(hFig, 'MenuBar', 'none', 'ToolBar', 'none');

% ===== SAVE IMAGE =====
else
    % If figure contains a video in an ActiveX control: must be extracted with screencapture()
    isForceScreencapture = ~isempty(FigureId) && strcmpi(FigureId.Type, 'Video') && ~isempty(iDS) && ismember(GlobalData.DataSet(iDS).Figure(iFig).Handles.PlayerType, {'VLC', 'WMPlayer'});
    % Headless display: we must print the figure
    if (GlobalData.Program.GuiLevel == -1)
        frameGfx.cdata = print(hFig, '-noui', '-r0', '-RGBImage');
    % Get figure bitmap
    elseif (bst_get('MatlabVersion') >= 804) && ~isForceScreencapture
        % Matlab function getframe() was finally fixed in R2014b
        frameGfx = getframe(hFig);
        % MRI Viewer: We need to update the figure to redraw all the java objects
        if ~isempty(FigureId) && strcmpi(FigureId.Type, 'MriViewer')
            figure_mri('ResizeCallback', hFig);
        end
    else
        % frameGfx = getscreen(hFig);
        figPos = get(hFig, 'Position');
        decoSize = bst_get('DecorationSize');
        % Get the screen definition
        ScreenDef = bst_get('ScreenDef');
        % Single screen (TODO: Handle the cases where the two screens are organized in different ways)
        if (length(ScreenDef) == 1)
            ZoomFactor = ScreenDef(1).zoomFactor;
        elseif (figPos(1) < ScreenDef(1).matlabPos(1) + ScreenDef(1).matlabPos(3))
            ZoomFactor = ScreenDef(1).zoomFactor;
        else
            ZoomFactor = ScreenDef(2).zoomFactor;
        end
        % Altman's screencapture doesn't work with ZoomFactor higher than 1
        if (ZoomFactor > 1.01)
            figFix   = round(figPos .* ZoomFactor);
            decoSize = round(decoSize .* ZoomFactor);
            capturePos = [figFix(1) - figPos(1) + decoSize(1), ...
                          -figFix(4) + figPos(4) + decoSize(4), ...
                          figFix(3) - decoSize(3) + 1, ...
                          figFix(4) - decoSize(4)];
        else
            capturePos = [decoSize(1), ...
                          decoSize(4), ...
                          figPos(3), ...
                          figPos(4)];
        end
        % Capture figure
        frameGfx.cdata = screencapture(hFig, capturePos);
    end
    % Save image file or return it in argument
    if (nargout == 0) && ~isempty(imgFile)
        if strcmpi(imgFile, 'Viewer')
            view_image(frameGfx.cdata);
        else
            out_image(imgFile, frameGfx.cdata);
        end
    else
        varargout{1} = frameGfx.cdata;
    end
end

%% ===== RESTORE IMAGE ======
% Delete created label
if ~isempty(imgLegend)
    delete(hAxesLabel);
end
% Show the buttons again
if ~isempty(hButtons)
    for i = 1:length(hButtons)
        set(hButtons(i), 'Visible', isVisible{i});
    end
end

author: brainstorm-tools


text:
function [data, label] = normalize_data(data, label, method)

% Data normalization
%
% Inputs:
%       data                data of size dxn, where d is dimension and n is number of sets
%       label               label data of size 1xn, n is number of labels.
% Output:
%       data                normalized data of size dxn, where d is dimension and n is number of sets
%       label               label data of size 1xn, n is number of labels. (no change!)
%
%
% Created by H.Kasai on July 04, 2017

    switch method
        case 'mean_std'
            data = data - repmat(mean(data),[size(data,1),1]);
            data = data./repmat(sqrt(sum(data.^2)),[size(data,1),1]);
        case 'mean'
            data = data - repmat(mean(data),[size(data,1),1]);
        case 'dimension_mean'
            data = data - repmat(mean(data,2),[1, size(data,2)]);
        case 'none'
            % Do nothing.
        case 'std' % normalize each column data to l2-norm
            %data = normc(data);
            data = data./repmat(sqrt(sum(data.^2)),[size(data,1),1]);
        case 'sample_mean_per_class'
            classes = unique(label);
            class_num = length(classes);
            for j = 1 : class_num
                train_idx = find(label == classes(j));
                class_mean = mean(data(:,train_idx), 2);
                data(:,train_idx) = data(:,train_idx) - repmat(class_mean, 1, length(train_idx));
            end
        otherwise

    end

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function test_ft_badsegment

% WALLTIME 00:10:00
% MEM 2gb

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

ft_debug off
interactive = false;

%%

cfg = [];
cfg.dataset = dataset;
cfg.channel = 'MEG';
cfg.trl = [1 7*3*300 0]; % first 7 trials, each of 3 seconds
cfg.continuous = 'yes';
cfg.demean = 'yes';
datacnt = ft_preprocessing(cfg);

cfg = [];
cfg.length = 1;
datatrl = ft_redefinetrial(cfg, datacnt);

%%

cfg = [];
cfg.method = 'template';
neighbours = ft_prepare_neighbours(cfg, datatrl);

%%

if interactive
  % this can be used to determine the thresholds
  cfg = [];
  cfg.method = 'summary';
  cfg.neighbours = neighbours;
  cfg.metric = 'var';
  data_rejected = ft_rejectvisual(cfg, datatrl);
end

%%
% each method has its own optimal threshold
% FT_REJECTVISUAL can be used to determine the appropriate value

cfg = [];
cfg.feedback = 'yes';
cfg.neighbours = neighbours;

cfg.metric = 'var';
cfg.threshold = 1.5e-25;

% cfg.metric = 'std';
% cfg.threshold = 4.2e-13;
%
% cfg.metric = 'max';
% cfg.threshold = 1.5e-12;
%
% cfg.metric = 'min';
% cfg.threshold = -12e-13;
%
% cfg.metric = 'range';
% cfg.threshold = 3000e-15;
%
% cfg.metric = 'zvalue';
% cfg.threshold = 3;
%
% cfg.metric = 'neighbexpvar';
% cfg.threshold = 0.5;
%
% cfg.metric = 'neighbstdratio';
% cfg.threshold = 0.5;
%
% cfg.metric = 'neighbcorr';
% cfg.threshold = 0.6;
%
% these options can be varied for neighbstdratio and neighbcorr
% cfg.nbdetect = 'median';
% cfg.nbdetect = 'any';
% cfg.nbdetect = 'most';
% cfg.nbdetect = 'all';

outcfg = ft_badsegment(cfg, datatrl);

%%

cfg = [];
cfg.artfctdef.badsegment.artifact = outcfg.artfctdef.badsegment.artifact;
cfg.continuous = 'yes';
cfg.blocksize = 21;
cfg.colorgroups = ones(size(datatrl.label));
cfg.linecolor = 'br';
cfg.layout = 'CTF151.lay';
ft_databrowser(cfg, datatrl);

%%

cfg = [];
cfg.artfctdef.badsegment.artifact = outcfg.artfctdef.badsegment.artifact;
data_rejected = ft_rejectartifact(cfg, datatrl);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [normed_data] = normalization_all(data)
% Calculating mean and standard deviation on whole image, then
% normalization, to keep consistent

mean_val=mean(data,1);
sigma_val=std(data,0,1);

normed_data=bsxfun(@minus,data,mean_val);
normed_data=bsxfun(@rdivide,normed_data,sigma_val);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function exact = hb_exact_read ( input_unit, nrow, nrhs, rhscrd, rhsfmt, rhstyp )

%*****************************************************************************80
%
%% HB_EXACT_READ reads the exact solution vectors in an HB file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992.
%
%  Parameters:
%
%    Input, integer INPUT_UNIT, the unit from which data is read.
%
%    Input, integer NROW, the number of rows or variables.
%
%    Input, integer NRHS, the number of right hand sides.
%
%    Input, integer RHSCRD, the number of lines in the file for
%    right hand sides.
%
%    Input, character ( len = 20 ) RHSFMT, the format for reading values
%    of the right hand side.
%
%    Input, character ( len = 3 ) RHSTYP, the right hand side type.
%    First character is F for full storage or M for same as matrix.
%    Second character is G if starting "guess" vectors are supplied.
%    Third character is X if exact solution vectors are supplied.
%    Ignored if NRHS = 0.
%
%    Output, real EXACT(NROW,NRHS), the exact solution vectors.
%
  exact = [];

  if ( 0 < rhscrd )

    if ( rhstyp(3) == 'X' )

      [ p, code, w, m ] = s_to_format ( rhsfmt );

      line_num = 1 + floor ( ( nrow - 1 ) / p );

      for irhs = 1 : nrhs

        jhi = 0;

        for i = 1 : line_num

          line = fgetl ( input_unit );
          jlo = jhi + 1;
          jhi = min ( jlo + p - 1, nrow );

          khi = 0;

          for j = jlo : jhi
            klo = khi + 1;
            khi = min ( klo + w - 1, length ( line ) );
            s = line(klo:khi);
            exact(j,irhs) = s_to_r8 ( s );
          end
        end
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function channel = serverBind(hostName, portNum)
import java.io.*;
import java.net.*;
import java.nio.*;
import java.nio.channels.*;

soRcvBuffSize	= 1024*1024*2;		%2 MByte
soSndBuffSize	= 1024*1024*2;
	try
		%==================================================
		% Configure the server and bind it to a port.
		%==================================================
		isa		= InetSocketAddress(hostName,portNum);
		channel = ServerSocketChannel.open();
		sock	= channel.socket();
		sock.setReceiveBufferSize(soRcvBuffSize);
		sock.setReuseAddress(true);
		channel.configureBlocking(false);
		sock.bind(isa);

	catch
		s	= lasterror;

		try
			channel.close();
		catch
			%do nothing
		end

		rethrow(s);
	end
end

author: Sable

repo: mcbench-benchmarks


text:
function [mesh_merge] = meshSuperquadrics(x, varargin)

n = size(x, 1);

% with tapering or not
taper = false;
% color = 'r';
arclength = 0.02;

for k = 1 : size(varargin, 2)
    if strcmp(varargin{k}, 'Taper')
        taper = varargin{k + 1};
    end
%     if strcmp(varargin{k}, 'Color')
%         color = varargin{k + 1};
%     end
    if strcmp(varargin{k}, 'Arclength')
        arclength= varargin{k + 1};
    end
end

mesh = cell(n, 1);
v_num = zeros(n, 1);
f_num = zeros(n, 1);

for i = 1 : n
    [x_mesh, y_mesh, z_mesh] = singleMeshSuperquadrics(x(i, :), ...
        'Arclength', arclength, 'Taper', taper);
    [mesh{i}.f, mesh{i}.v] = mesh2tri(x_mesh, y_mesh, z_mesh,'b'); % f b x
    v_num(i) = size(mesh{i}.v, 1);
    f_num(i) = size(mesh{i}.f, 1);
    if i > 1
        mesh{i}.f = mesh{i}.f + sum(v_num(1 : i - 1));
    end
end
mesh_merge.v = zeros(sum(v_num), 3);
mesh_merge.f = zeros(sum(f_num), 3);
mesh_merge.v(1 : v_num(1), :) = mesh{1}.v;
mesh_merge.f(1 : f_num(1), :) = mesh{1}.f;
for i = 2 : n
    mesh_merge.v(...
        sum(v_num(1 : i - 1)) + 1 : sum(v_num(1 : i - 1)) + v_num(i), :) = ...
        mesh{i}.v;
    mesh_merge.f(...
        sum(f_num(1 : i - 1)) + 1 : sum(f_num(1 : i - 1)) + f_num(i), :) = ...
        mesh{i}.f;
end
end

author: ChirikjianLab

repo: Marching-Primitives


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

% Script for generating a 2d SABR density using the Wu expansion

% base scenario
clear; clc;
t=0; T=1;
f = 0.03;
alpha = 0.075;
beta = 0.5;
nu = 0.05;
rho = 0;

% All expansions 0th to 2nd order of Wu
%sabrdensity = @(x,y) psabr30(t,T,f,x,alpha,y,beta,nu,rho);
% sabrdensity = @(x,y) psabr31(t,T,f,x,alpha,y,beta,nu,rho);
sabrdensity = @(x,y) psabr32(t,T,f,x,alpha,y,beta,nu,rho);

lcoeff = 1; ucoeff = 1;
lowerbound = f-lcoeff*f;
upperbound = f + ucoeff*f;

xvals = lowerbound:.001:upperbound;
yvals = 0.05:.001:.1;

zvals2d = density2d(xvals,yvals,sabrdensity);

% Density 2d
figure; surf(yvals,xvals,zvals2d);


text:
function ind = subsind(gB,subs)

% select by isIndexed
if check_option(subs,'indexed')
  ind = gB.isIndexed;
  subs = delete_option(subs,'indexed');
else
  ind = true(length(gB),1);
end

% select for minerals
isMineralName = cellfun(@ischar,subs);
if any(isMineralName)
  ind = gB.hasPhase(subs{isMineralName});
end
subs = subs(~isMineralName);

% select for grains
isGrain = cellfun(@(x) isa(x,'grain2d'),subs);
if any(isGrain)
  ind = ind & gB.hasGrain(subs{isGrain});
end
subs = subs(~isGrain);

% other indexing
for i = 1:length(subs)

  if isa(subs{i},'logical')

    sub = any(subs{i}, find(size(subs{i}')==max(size(ind)),1));

    ind = ind & reshape(sub,size(ind));

  elseif isnumeric(subs{i})

    if any(subs{i} <= 0 | subs{i} > length(gB))
      error('Out of range; index must be a positive integer or logical.')
    end

    ind = subs{i};
    return

  elseif isa(subs{i},'polygon')

    ind = ind & inpolygon(gB,subs{i})';

  end
end
end

author: mtex-toolbox

repo: mtex

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


text:
function Tree = forwardPropRAE(allKids, W1,W2,W3,W4,b1,b2,b3, Wcat, bcat, alpha_cat, updateWcat, beta, words_embedded, label, hiddenSize, sl, freq, f, f_prime)

collapsed_sentence = 1:sl;
Tree = tree2;
Tree.pp = zeros((2*sl-1),1);
Tree.nodeScores = zeros(2*sl-1,1);
Tree.nodeNames = 1:(2*sl-1);
Tree.kids = zeros(2*sl-1,2);
Tree.numkids = ones(2*sl-1,1);
Tree.node_y1c1 = zeros(hiddenSize,2*sl-1);
Tree.node_y2c2 = zeros(hiddenSize,2*sl-1);
Tree.freq = zeros(2*sl-1,1);
Tree.nodeFeatures = [words_embedded zeros(hiddenSize, sl-1)];
Tree.nodeFeatures_unnormalized = [words_embedded zeros(hiddenSize, sl-1)];
Tree.nodeDelta_out1 = zeros(hiddenSize,2*sl-1);
Tree.nodeDelta_out2 = zeros(hiddenSize,2*sl-1);
Tree.parentDelta = zeros(hiddenSize,2*sl-1);

if updateWcat
    % Classification Error
    nodeUnder = ones(2*sl-1,1);

    for i = sl+1:2*sl-1
        kids = allKids(i,:);
        n1 = nodeUnder(kids(1));
        n2 = nodeUnder(kids(2));
        nodeUnder(i) = n1+n2;
    end

    [cat_size, ~] = size(Wcat);
    Tree.catDelta = zeros(cat_size, 2*sl-1);
    Tree.catDelta_out = zeros(hiddenSize,2*sl-1);

    % classifier on single words
    sm = sigmoid(Wcat*words_embedded + bcat);
    lbl_sm = (1-alpha_cat).*(label(:,ones(sl,1)) - sm);
    Tree.nodeScores(1:sl) = 1/2*(lbl_sm.*(label(:,ones(sl,1)) - sm));
    Tree.catDelta(:, 1:sl) = -(lbl_sm).*sigmoid_prime(sm);

    for i = sl+1:2*sl-1

        kids = allKids(i,:);

        c1 = Tree.nodeFeatures(:,kids(1));
        c2 = Tree.nodeFeatures(:,kids(2));

        % Eq. (2) in the paper: p = f(W(1)[c1; c2] + b(1))
        p = f(W1*c1 + W2*c2 + b1);

        % See last paragraph in Section 2.3
        p_norm1 = p./norm(p);

        % Eq. (7) in the paper (for special case of 1d label)
        sm = sigmoid(Wcat*p_norm1 + bcat);

        lbl_sm = beta * (1-alpha_cat).*(label - sm);
        Tree.catDelta(:, i) = -(lbl_sm).*sigmoid_prime(sm);
        J = 1/2*(lbl_sm'*(label - sm));

        Tree.nodeFeatures(:,i) = p_norm1;
        Tree.nodeFeatures_unnormalized(:,i) = p;
        Tree.nodeScores(i) = J;
        Tree.numkids = nodeUnder;

    end

    Tree.kids = allKids;

else

    % Reconstruction Error
    for j=1:sl-1

        [~, size2] = size(words_embedded);
        c1 = words_embedded(:,1:end-1);
        c2 = words_embedded(:,2:end);

        freq1 = freq(1:end-1);
        freq2 = freq(2:end);

        p = f(W1*c1 + W2*c2 + b1(:,ones(size2-1,1)));
        p_norm1 = bsxfun(@rdivide,p,sqrt(sum(p.^2,1)));

        y1_unnormalized = f(W3*p_norm1 + b2(:,ones(size2-1,1)));
        y2_unnormalized = f(W4*p_norm1 + b3(:,ones(size2-1,1)));

        y1 = bsxfun(@rdivide,y1_unnormalized,sqrt(sum(y1_unnormalized.^2,1)));
        y2 = bsxfun(@rdivide,y2_unnormalized,sqrt(sum(y2_unnormalized.^2,1)));

        y1c1 = alpha_cat*(y1-c1);
        y2c2 = alpha_cat*(y2-c2);

        % Eq. (4) in the paper: reconstruction error
        J = 1/2*sum((y1c1).*(y1-c1) + (y2c2).*(y2-c2));

        % finding the pair with smallest reconstruction error for constructing tree
        [J_min J_minpos] = min(J);

        Tree.node_y1c1(:,sl+j) = y1c1(:,J_minpos);
        Tree.node_y2c2(:,sl+j) = y2c2(:,J_minpos);
        Tree.nodeDelta_out1(:,sl+j) = f_prime(y1_unnormalized(:,J_minpos)) * y1c1(:,J_minpos);
        Tree.nodeDelta_out2(:,sl+j) = f_prime(y2_unnormalized(:,J_minpos)) * y2c2(:,J_minpos);

        words_embedded(:,J_minpos+1)=[];
        words_embedded(:,J_minpos)=p_norm1(:,J_minpos);
        Tree.nodeFeatures(:, sl+j) = p_norm1(:,J_minpos);
        Tree.nodeFeatures_unnormalized(:, sl+j) = p(:,J_minpos);
        Tree.nodeScores(sl+j) = J_min;
        Tree.pp(collapsed_sentence(J_minpos)) = sl+j;
        Tree.pp(collapsed_sentence(J_minpos+1)) = sl+j;
        Tree.kids(sl+j,:) = [collapsed_sentence(J_minpos) collapsed_sentence(J_minpos+1)];
        Tree.numkids(sl+j) = Tree.numkids(Tree.kids(sl+j,1)) + Tree.numkids(Tree.kids(sl+j,2));

        freq(J_minpos+1) = [];
        freq(J_minpos) = (Tree.numkids(Tree.kids(sl+j,1))*freq1(J_minpos) + Tree.numkids(Tree.kids(sl+j,2))*freq2(J_minpos))/(Tree.numkids(Tree.kids(sl+j,1))+Tree.numkids(Tree.kids(sl+j,2)));

        collapsed_sentence(J_minpos+1)=[];
        collapsed_sentence(J_minpos)=sl+j;

    end
end
end

author: jacoxu


text:
function [ImTrans,tau] = preAlign(ImData)
% for pre-alignment
numFrame = size(ImData,3);
ImTrans = ImData;
IDcenter = round(numFrame/2);
tau = zeros(6,numFrame);
numLevel = ceil(log2(max([size(ImData,1),size(ImData,2)])/50))+1;
for i = IDcenter-1:-1:1
    disp(['frame: ',num2str(i)]);
    [ImTrans(:,:,i),tau(:,i)] = regMGNC(ImTrans(:,:,i+1),ImData(:,:,i),tau(:,i+1),numLevel);
end
for i = IDcenter+1:+1:numFrame
    disp(['frame: ',num2str(i)]);
    [ImTrans(:,:,i),tau(:,i)] = regMGNC(ImTrans(:,:,i-1),ImData(:,:,i),tau(:,i-1),numLevel);
end
end

author: andrewssobral

repo: lrslibrary


text:
function [lbm] = amu2lbm(amu)
% Convert mass from atomic mass units to pounds-mass.
% Chad Greene 2012
lbm = amu*3.660864489409e-27;


text:
function eventLog = ma_executeDVManeuver_finite_steered(dVVectNTW, thruster, initialState, eventNum, celBodyData)
%ma_executeDVManeuver_finite_steered Summary of this function goes here
%   Detailed explanation goes here
    ut = initialState(1);
    rVectUT = initialState(2:4);
    vVectUT = initialState(5:7);
    bodyID = initialState(8);
    mass0 = sum(initialState(9:12));

    bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);
    gmu = bodyInfo.gm;

    isp = thruster.isp;
    thrust = thruster.thrust;
    dvMag = norm(dVVectNTW);
    massFlowRate = getMdotFromThrustIsp(thrust, isp);
    burnDur = getBurnDuration(mass0, massFlowRate, isp, dvMag);

    if(burnDur == 0)
        eventLog = initialState;
        eventLog(13) = eventNum;
        return;
    end

    dVVectNTW = reshape(dVVectNTW,3,1);
    rhs = @(t,x) propBurnRHS(t, x, gmu, massFlowRate, dVVectNTW, thrust);
    tspan = [ut, ut + burnDur];
    x0 = [rVectUT, vVectUT, mass0];
    options = odeset('RelTol',1E-6,'AbsTol',1E-6);
    [T,Y] = ode45(rhs,tspan,x0,options);

    dM = mass0 - Y(:,7);

    eventLog = zeros(length(T), length(initialState));
    eventLog(:,1) = T;
    eventLog(:,2:7) = Y(:,1:6);
    eventLog(:,8) = bodyID;
    eventLog(:,9) = initialState(9);
    eventLog(:,13) = eventNum;
    switch(thruster.propType)
        case 'fuelOxMass'
            eventLog(:,10) = initialState(10) - dM;
            eventLog(:,11) = initialState(11);
            eventLog(:,12) = initialState(12);
        case 'monoMass'
            eventLog(:,10) = initialState(10);
            eventLog(:,11) = initialState(11) - dM;
            eventLog(:,12) = initialState(12);
        case 'xenonMass'
            eventLog(:,10) = initialState(10);
            eventLog(:,11) = initialState(11);
            eventLog(:,12) = initialState(12) - dM;
    end
end

function xdot = propBurnRHS(~, x, gmu, massFlowRate, dVVectNTW, thrust)
    rVect = x(1:3);
    r = norm(rVect);
    vVect = x(4:6);
    m = x(7);

    gAccel = -(gmu/r^3)*rVect; %km^3/s^2 / km^2 = km/s^2

    dVVectECI = getNTW2ECIdvVect(dVVectNTW, rVect', vVect');
    dvNorm = normVector(dVVectECI);
    tForce = (thrust * dvNorm)/1000; %kN = tons * m/s^2/1000 = km/s^2

    xdot(1:3) = vVect;
    xdot(4:6) = gAccel + tForce/m;
    xdot(7) = -massFlowRate;

    xdot = xdot';
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [cIX_,gIX_,num] = ArtifactAnalysis(cIX,gIX,CellXYZ,absIX)
U = unique(gIX);
cIX_abs = absIX(cIX);

Std = zeros(length(U),3);
Range = zeros(length(U),3);
TF = zeros(length(gIX),1);
for i = 1:length(U),
    IX = gIX==U(i);
    YXZ = CellXYZ(cIX_abs(IX),:);
    for j = 1:3,
        % remove outliers
        data = YXZ(:,j);
        percntiles = prctile(data,[5 95]); %5th and 95th percentile
        outlierIX = data < percntiles(1) | data > percntiles(2);
        data(outlierIX) = [];
        if ~isempty(data),
            % calculate standard deviation
            Std(i,j) = std(data);
            Range(i,j) = max(data)-min(data);
        end
    end
    thres_std = 0.5;
    thres_range = 100;
    tf1 = Std(i,1)<thres_std;
    tf2 = Std(i,2)<thres_std;
    tf3 = Std(i,3)<thres_std;
    tf = tf1 || tf2 || tf3 && max(Range(i,:))>thres_range;
    TF(IX) = tf;
end

cIX_ = cIX(find(TF));
gIX_ = gIX(find(TF));
num = length(unique(gIX_));
end

author: xiuyechen

repo: FishExplorer


text:
function PlotContour3D(Nout, u, levels)

% function PlotContour3D(Nout, u, levels)
% Purpose: generic routine to plot contours for triangulated data

Globals3D;

% create fine output grid
[ro,so,to] = EquiNodes3D(Nout);
Vo = Vandermonde3D(N, ro, so, to);
interp = Vo*invV;

xo = interp*x; yo = interp*y; zo = interp*z; uo = interp*u;
xo = xo(:); yo = yo(:); zo = zo(:); uo = uo(:);

% form symmetric tetrahedralization of local nodes
startrow = zeros(Nout+1, Nout+1);

sk = 1;
for i=0:Nout
  for j=0:Nout-i
    startrow(j+1, i+1) = sk;
    sk = sk + Nout+1-i-j;
  end
end

% contruct tetrahedralization
ltet = zeros(N*N*N,4);

sk = 1;
for i=0:Nout
  for j=0:Nout-i
    for k=0:Nout-i-j-1
      % Add Tet 1
      ltet(sk,1) = startrow(j+1, i+1)+k;
      ltet(sk,2) = startrow(j+1, i+1)+k+1;
      ltet(sk,3) = startrow(j+2, i+1)+k;
      ltet(sk,4) = startrow(j+1, i+2)+k;
      sk = sk+1;

      if(k < Nout-i-j-1)
	% Add Tet 2
        ltet(sk,1) = startrow(j+1, i+1)+k+1;
        ltet(sk,2) = startrow(j+2, i+1)+k;
        ltet(sk,3) = startrow(j+1, i+2)+k;
        ltet(sk,4) = startrow(j+1, i+2)+k+1;
        sk = sk+1;

        % Add Tet 3
        ltet(sk,1) = startrow(j+1, i+1)+k+1;
        ltet(sk,2) = startrow(j+2, i+1)+k+1;
        ltet(sk,3) = startrow(j+2, i+1)+k;
        ltet(sk,4) = startrow(j+1, i+2)+k+1;
        sk = sk+1;

        % Add Tet 4
        ltet(sk,1) = startrow(j+1, i+2)+k;
        ltet(sk,2) = startrow(j+2, i+2)+k;
        ltet(sk,3) = startrow(j+1, i+2)+k+1;
        ltet(sk,4) = startrow(j+2, i+1)+k;
	    sk = sk+1;

        % Add Tet 5
        ltet(sk,1) = startrow(j+2, i+1)+k;
        ltet(sk,2) = startrow(j+2, i+1)+k+1;
        ltet(sk,3) = startrow(j+2, i+2)+k;
        ltet(sk,4) = startrow(j+1, i+2)+k+1;
	    sk = sk+1;
      end

      if(k < Nout-i-j-2)
	% Add Tet 6
        ltet(sk,1) = startrow(j+1, i+2)+k+1;
        ltet(sk,2) = startrow(j+2, i+1)+k+1;
        ltet(sk,3) = startrow(j+2, i+2)+k;
        ltet(sk,4) = startrow(j+2, i+2)+k+1;
        sk = sk+1;
      end
    end
  end
end

Nltet = sk-1;

% create global tet mesh
Npo = size(ro,1);
tet = zeros(K*Nltet, 4);
ids = 1:Nltet;
for k=1:K
  tet(ids, :) = ltet + (k-1)*Npo;
  ids = ids+Nltet;
end

Nlevels = length(levels);

minu = min(uo(tet), [], 2);
maxu = max(uo(tet), [], 2);

v1 = tet(:,1); v2 = tet(:,2); v3 = tet(:,3); v4 = tet(:,4);
u1 = uo(v1);    u2 = uo(v2);    u3 = uo(v3);    u4 = uo(v4);
x1 = xo(v1);    x2 = xo(v2);    x3 = xo(v3);    x4 = xo(v4);
y1 = yo(v1);    y2 = yo(v2);    y3 = yo(v3);    y4 = yo(v4);
z1 = zo(v1);    z2 = zo(v2);    z3 = zo(v3);    z4 = zo(v4);

trix = []; triy = []; triz = []; triu = [];;

for n=1:Nlevels
  lev = levels(n);

  % find candidate elements
  ks = find(minu<= lev & maxu >= lev);

  Nks = length(ks);

  if(Nks>0)

    fc = zeros(Nks,6);

    u1ks = u1(ks); u2ks = u2(ks); u3ks = u3(ks); u4ks = u4(ks);

    % find local coordinate at each of 6 edges
    c1 = (lev-u1ks)./(u2ks-u1ks);  fc(:,1) = (min(u1ks,u2ks)<lev & max(u1ks,u2ks)>lev);
    c2 = (lev-u2ks)./(u3ks-u2ks);  fc(:,2) = (min(u2ks,u3ks)<lev & max(u2ks,u3ks)>lev);
    c3 = (lev-u3ks)./(u1ks-u3ks);  fc(:,3) = (min(u3ks,u1ks)<lev & max(u3ks,u1ks)>lev);
    c4 = (lev-u1ks)./(u4ks-u1ks);  fc(:,4) = (min(u1ks,u4ks)<lev & max(u1ks,u4ks)>lev);
    c5 = (lev-u2ks)./(u4ks-u2ks);  fc(:,5) = (min(u2ks,u4ks)<lev & max(u2ks,u4ks)>lev);
    c6 = (lev-u3ks)./(u4ks-u3ks);  fc(:,6) = (min(u3ks,u4ks)<lev & max(u3ks,u4ks)>lev);

    % find triangle intersection
    tris = find(sum(fc, 2)==3); Ntris = length(tris);  tfc = fc(tris,:); ids = find(tfc');
    ktris = ks(tris);

    % trim list
    tc1 = c1(tris); tc2 = c2(tris); tc3 = c3(tris); tc4 = c4(tris); tc5 = c5(tris); tc6 = c6(tris);

    tx1 = x1(ktris);   ty1 = y1(ktris);   tz1 = z1(ktris);
    tx2 = x2(ktris);   ty2 = y2(ktris);   tz2 = z2(ktris);
    tx3 = x3(ktris);   ty3 = y3(ktris);   tz3 = z3(ktris);
    tx4 = x4(ktris);   ty4 = y4(ktris);   tz4 = z4(ktris);

    xc = zeros(Ntris, 6); yc = zeros(Ntris, 6); zc = zeros(Ntris, 6);
    xc(:,1) = (1-tc1).*tx1 + tc1.*tx2; yc(:,1) = (1-tc1).*ty1 + tc1.*ty2; zc(:,1) = (1-tc1).*tz1 + tc1.*tz2;
    xc(:,2) = (1-tc2).*tx2 + tc2.*tx3; yc(:,2) = (1-tc2).*ty2 + tc2.*ty3; zc(:,2) = (1-tc2).*tz2 + tc2.*tz3;
    xc(:,3) = (1-tc3).*tx3 + tc3.*tx1; yc(:,3) = (1-tc3).*ty3 + tc3.*ty1; zc(:,3) = (1-tc3).*tz3 + tc3.*tz1;
    xc(:,4) = (1-tc4).*tx1 + tc4.*tx4; yc(:,4) = (1-tc4).*ty1 + tc4.*ty4; zc(:,4) = (1-tc4).*tz1 + tc4.*tz4;
    xc(:,5) = (1-tc5).*tx2 + tc5.*tx4; yc(:,5) = (1-tc5).*ty2 + tc5.*ty4; zc(:,5) = (1-tc5).*tz2 + tc5.*tz4;
    xc(:,6) = (1-tc6).*tx3 + tc6.*tx4; yc(:,6) = (1-tc6).*ty3 + tc6.*ty4; zc(:,6) = (1-tc6).*tz3 + tc6.*tz4;
    xc = xc'; yc = yc'; zc = zc';

    ids = reshape(ids, 3, Ntris);

    trix = [trix, xc(ids)];
    triy = [triy, yc(ids)];
    triz = [triz, zc(ids)];
    triu = [triu, lev*ones(size(ids))];

    % find quadrilateral intersection
    quads = find(sum(fc, 2)==4); Nquads = length(quads);  qfc = fc(quads,:); ids = find(qfc');
    kquads = ks(quads);

    % quad list
    qc1 = c1(quads); qc2 = c2(quads); qc3 = c3(quads); qc4 = c4(quads); qc5 = c5(quads); qc6 = c6(quads);

    qx1 = x1(kquads);   qy1 = y1(kquads);   qz1 = z1(kquads);
    qx2 = x2(kquads);   qy2 = y2(kquads);   qz2 = z2(kquads);
    qx3 = x3(kquads);   qy3 = y3(kquads);   qz3 = z3(kquads);
    qx4 = x4(kquads);   qy4 = y4(kquads);   qz4 = z4(kquads);

    xc = zeros(Nquads, 6); yc = zeros(Nquads, 6); zc = zeros(Nquads, 6);
    xc(:,1) = (1-qc1).*qx1 + qc1.*qx2; yc(:,1) = (1-qc1).*qy1 + qc1.*qy2; zc(:,1) = (1-qc1).*qz1 + qc1.*qz2;
    xc(:,2) = (1-qc2).*qx2 + qc2.*qx3; yc(:,2) = (1-qc2).*qy2 + qc2.*qy3; zc(:,2) = (1-qc2).*qz2 + qc2.*qz3;
    xc(:,3) = (1-qc3).*qx3 + qc3.*qx1; yc(:,3) = (1-qc3).*qy3 + qc3.*qy1; zc(:,3) = (1-qc3).*qz3 + qc3.*qz1;
    xc(:,4) = (1-qc4).*qx1 + qc4.*qx4; yc(:,4) = (1-qc4).*qy1 + qc4.*qy4; zc(:,4) = (1-qc4).*qz1 + qc4.*qz4;
    xc(:,5) = (1-qc5).*qx2 + qc5.*qx4; yc(:,5) = (1-qc5).*qy2 + qc5.*qy4; zc(:,5) = (1-qc5).*qz2 + qc5.*qz4;
    xc(:,6) = (1-qc6).*qx3 + qc6.*qx4; yc(:,6) = (1-qc6).*qy3 + qc6.*qy4; zc(:,6) = (1-qc6).*qz3 + qc6.*qz4;
    xc = xc'; yc = yc'; zc = zc';

    ids1 = reshape(ids, 4, Nquads);
    ids2 = ids1([1;4;2;3],:);
    ids3 = ids1([1;3;4;2],:);

    % find length perimeters bounded by 4 nodes
    perims1 = QuadPerimeter3D(xc(ids1),yc(ids1),zc(ids1));
    perims2 = QuadPerimeter3D(xc(ids2),yc(ids2),zc(ids2));
    perims3 = QuadPerimeter3D(xc(ids3),yc(ids3),zc(ids3));

    % locate minimum perimeter quad
    idsA = find(perims1<min(perims2, perims3));
    idsB = find(perims2<min(perims1, perims3));
    idsC = find(perims3<min(perims1, perims2));

    % use only minimum perimeter quads
    ids  = [ids1(:,idsA),ids2(:,idsB),ids3(:,idsC)];

    % divide into two triangles (ok since quad must be convex)
    ids1 = [ids([1;2;3],:),ids([1;3;4],:)];

    % add triangles to list
    trix = [trix, xc(ids1)];
    triy = [triy, yc(ids1)];
    triz = [triz, zc(ids1)];
    triu = [triu, lev*ones(size(ids1))];

  end
end

ha = patch(trix,triy,triz,triu); shading interp
set(ha, 'EdgeColor', 'none')

camlight; camlight; material shiny; lighting gouraud; view(3); axis equal;
axis off;

return;


text:
function [V,mrGrayFileName]=analyze2mrGrayFunctional(analHeaderFileName,endianType);
% function [V,mrGrayFileName]=analyze2mrGray(analHeaderFileName,endianType);
% Converts Analyze7.5 images to Stanford VISTA lab / mrGray vAnatomy.dat format.
% Analyze headers hold far more information that mrGray volumes. We just toss this information out when we write the mrGray file
% but this routine returns the Analyze header info in 'V'.
% Uses the SPM99 spm_vol command to parse the header. Uses misc. VISTA lab routines to write the volume.
% See also mrGray2Analyze
% If you pass in the (optional) endianType it will use this when reading 16bit data files. Otherwise it assumes 'ieee-be' (non-intel).

% ARW 081701
% ARW Added >16 bit data handling 082101

% Check whether we're reading big or little endian data.
% Default is 'native'; DKP/ Arizona data are 'ieee-be'
if (~exist('endianType','var'))
    endianType='ieee-be';
end

% Read the header
V=spm_vol(analHeaderFileName);

% Now read in the image data
fid=fopen(V.fname,'rb');
if (~fid)
    error (['Could not read header file: ',analHeaderFileName]);
end

im=fread(fid,inf,'uint8');
fclose(fid);
imLen=length(im);

if(imLen~=prod(V.dim(1:3)))
    disp('Wrong data size for 8-bit. Trying 16bit instead');
    fid=fopen(V.fname,'rb',endianType);

    im=fread(fid,inf,'uint16');
    fclose(fid);
    imLen=length(im);
    if(imLen~=prod(V.dim(1:3)))
        error ('Not 16 bit either. Quitting');
    end

    % Now we have to convert this to 8-bits
    % We're just going to scale for now..
    % Could rescale to (a) The max value, (b) V.dim(4) (the bit depth) , (c) some value we cleverly compute to throw away,say, 5% of the outliers in the histogram
    % Choose (a) for now.
    im=im-min(im(:));
    im=fix(im./max(im(:))*255);

end

img=reshape(im,V.dim(1:3));
mmPerVox=diag(V.mat(1:3,1:3));
mmPerVox=mmPerVox;

% Need to rearrange the image. Analyze data run with sag and axial directions flipped. Also upside down in mrGray...
newimg=zeros(V.dim(3),V.dim(2),V.dim(1));
for thisIm=1:(V.dim(2))
   newimg(:,thisIm,:)=flipud(squeeze(img(:,thisIm,:))');
end

% Now write out a mrGray file
mrGrayFileName = writeVolAnat(newimg, mmPerVox);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [varargout]=patch_dual(varargin)

% function [Vd,Fd,Fds]=patch_dual(V,F,fixBoundaryOption)
% ------------------------------------------------------------------------
% Computes the dual of the input tesselation defined by the vertices V and
% faces F.
%
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2013/04/28
% 2015/07/07 Updated for GIBBON
% 2018/04/06 Added color output
% 2018/04/06 Improved handling of edge cells
%------------------------------------------------------------------------
%%

switch nargin
    case 2
        V=varargin{1};
        F=varargin{2};
        fixBoundaryOption=1;
    case 3
        V=varargin{1};
        F=varargin{2};
        fixBoundaryOption=varargin{3};
end

%Cope with 2D input
if size(V,2)==2
    V(:,3)=0;
end

%%

if iscell(F)
    %     numFacesperSet=cellfun(@(x) size(x,1),F);

    N=[];
    vertexFaceConnectivity=[];
    numFacesPrevous=0;
    for q=1:1:numel(F)
        %The current face set
        f=F{q};

        %Get indices of face neighbourhood
        [vertexFaceConnectivity_f]=tesIND(f,V,0);

        %Offset face indices by previous
        vertexFaceConnectivity_f(vertexFaceConnectivity_f>0)=vertexFaceConnectivity_f(vertexFaceConnectivity_f>0)+numFacesPrevous;

        %Get face-vertex connectivity
        vertexFaceConnectivity=[vertexFaceConnectivity vertexFaceConnectivity_f];

        %Derive vertex normals
        N=[N; patchNormal(f,V)];

        numFacesPrevous=numFacesPrevous+size(f,1);
    end
    vertexFaceConnectivity=sort(vertexFaceConnectivity,2);
    [~,J,~] = find(vertexFaceConnectivity);
    vertexFaceConnectivity=vertexFaceConnectivity(:,min(J):end);
    vertexFaceConnectivity=full(vertexFaceConnectivity);

    L=vertexFaceConnectivity>0;

    Nv=ones(size(V,1),size(N,2));
    for q=1:1:size(N,2)
        nf=N(:,q);
        nv=nan(size(vertexFaceConnectivity));
        nv(L)=nf(vertexFaceConnectivity(L));
        nv=gnanmean(nv,2);
        Nv(:,q)=nv;
    end
else
    %Get indices of face neighbourhood
    [vertexFaceConnectivity]=tesIND(F,V,0);

    %Get input tesselation vertex-normals
    [~,~,Nv]=patchNormal(F,V);
end

%Create central face coordinates
if iscell(F)
    VF=[];
    for q=1:1:numel(F)
        VF=[VF; patchCentre(F{q},V)];
    end
else
    VF=patchCentre(F,V);
end

if fixBoundaryOption>0

    [Eb]=patchBoundary(F);
    if ~isempty(Eb)
        switch fixBoundaryOption
            case 1 %Include original boundary vertices
                ind_Eb=unique(Eb(:));

                %Create mid-edge vertices
                VE=zeros(size(Eb,1),size(V,2));
                for q=1:1:size(V,2)
                    X=V(:,q);
                    VE(:,q)=mean(X(Eb),2);
                end

                Vd=[VF; VE; ];

                %Add new points in IND_F
                vertexFaceConnectivity(Eb(:,1),end+1)=(1:size(VE,1))'+size(VF,1);
                ind1=vertexFaceConnectivity(ind_Eb,end);
                vertexFaceConnectivity(Eb(:,2),end+1)=(1:size(VE,1))'+size(VF,1);
                ind2=vertexFaceConnectivity(ind_Eb,end);
                vertexFaceConnectivity(ind_Eb,end+1)=(1:1:numel(ind_Eb))'+size(VF,1)+size(VE,1);

                VB=(Vd(ind1,:)+Vd(ind2,:))/2; %Use this average for now to avoid inverted elements VB=V(ind_Eb,:);

                %Collect point sets
                Vd=[VF; VE; VB];

                vertexFaceConnectivity=sort(vertexFaceConnectivity,2,'descend');
                indMax=find(sum(vertexFaceConnectivity,1)==0,1);
                vertexFaceConnectivity=vertexFaceConnectivity(:,1:indMax);

                replaceVb=1; %Replace on so average is replaced later by real boundary vertices
            case 2 % Do not include original boundary vertices
                %Create mid-edge vertices
                VE=zeros(size(Eb,1),size(V,2));
                for q=1:1:size(V,2)
                    X=V(:,q);
                    VE(:,q)=mean(X(Eb),2);
                end

                %Collect point sets
                Vd=[VF; VE; ];

                %Add new points in IND_F
                vertexFaceConnectivity(Eb(:,1),end+1)=(1:size(VE,1))'+size(VF,1);
                vertexFaceConnectivity(Eb(:,2),end+1)=(1:size(VE,1))'+size(VF,1);

                vertexFaceConnectivity=sort(vertexFaceConnectivity,2,'descend');
                indMax=find(sum(vertexFaceConnectivity,1)==0,1);
                vertexFaceConnectivity=vertexFaceConnectivity(:,1:indMax);

                replaceVb=0;
        end
    else
        Vd=VF;
        replaceVb=0;
    end
else
    Vd=VF;
    replaceVb=0;
end

%% Creating face array
[I,J,v] = find(vertexFaceConnectivity);

Xfd=accumarray({I,J},Vd(v,1),size(vertexFaceConnectivity),[],NaN);
Yfd=accumarray({I,J},Vd(v,2),size(vertexFaceConnectivity),[],NaN);
Zfd=accumarray({I,J},Vd(v,3),size(vertexFaceConnectivity),[],NaN);

Xfd_mean=gnanmean(Xfd,2);
Yfd_mean=gnanmean(Yfd,2);
Zfd_mean=gnanmean(Zfd,2);

Xfd=Xfd-Xfd_mean(:,ones(1,(size(Xfd,2))));
Yfd=Yfd-Yfd_mean(:,ones(1,(size(Yfd,2))));
Zfd=Zfd-Zfd_mean(:,ones(1,(size(Zfd,2))));
m=sqrt(Xfd.^2+Yfd.^2+Zfd.^2);

dx=Xfd.*Xfd(:,ones(size(Xfd,2),1));
dy=Yfd.*Yfd(:,ones(size(Yfd,2),1));
dz=Zfd.*Zfd(:,ones(size(Zfd,2),1));
d=(1./m).*(1./m(:,ones(size(dx,2),1))).*(dx+dy+dz);
A=real(acos(d));

%Make sure angle conforms to mesh side
%Cross product
Xx=Yfd.*Zfd(:,ones(size(Xfd,2),1))-Zfd.*Yfd(:,ones(size(Xfd,2),1));
Yx=Zfd.*Xfd(:,ones(size(Xfd,2),1))-Xfd.*Zfd(:,ones(size(Xfd,2),1));
Zx=Xfd.*Yfd(:,ones(size(Xfd,2),1))-Yfd.*Xfd(:,ones(size(Xfd,2),1));
mx=sqrt(Xx.^2+Yx.^2+Zx.^2);

%Correct wrap
dx=Xx.*Nv(:,ones(size(Xfd,2),1));
dy=Yx.*Nv(:,2*ones(size(Yfd,2),1));
dz=Zx.*Nv(:,3*ones(size(Zfd,2),1));
d=(1./mx).*(dx+dy+dz);
logic1=d>0;
logicAngle= (A<(pi) & logic1);
A(logicAngle)=(2*pi)-A(logicAngle);

%Sort based on angle
[~,J_sort]=sort(A,2);
I_sort=(1:1:size(J_sort,1))';
I_sort=I_sort(:,ones(1,size(J_sort,2)));
IND_sort = sub2ind(size(J_sort),I_sort,J_sort);

% Creating faces matrix
Fds=vertexFaceConnectivity(IND_sort);

% %Determine face order
% Xfdn=Xfd; Yfdn=Yfd; Zfdn=Zfd;
% for q=1:1:size(Xfd,1)
%     Vc=[Xfd(q,:); Yfd(q,:); Zfd(q,:)]';
%     logicValid=~any(isnan(Vc),2);
%     [R]=pointSetPrincipalDir(Vc(logicValid,:)); %Fit local coordinate system with 3rd direction pointing outward of local planar-ish region
%     Vcn=Vc*R; %Rotate coordinate set to prepare for 2D Delaunay based triangulation
%     Xfdn(q,:)=Vcn(:,1);
%     Yfdn(q,:)=Vcn(:,2);
%     Zfdn(q,:)=Vcn(:,3);
% end
%
% %Fix face order
% theta_n=atan2(Yfdn,Xfdn);
% [~,J_sort]=sort(theta_n,2);
% I_sort=(1:1:size(J_sort,1))';
% I_sort=I_sort(:,ones(1,size(J_sort,2)));
% IND_sort = sub2ind(size(J_sort),I_sort,J_sort);
%
% % Creating faces matrix
% Fds=vertexFaceConnectivity(IND_sort);

%% Splitting up into seperate face types and fix face normals
n_sum=sum(Fds>0,2);
faceTypesNum=unique(n_sum);
faceTypesNum=faceTypesNum(faceTypesNum>0);
Fd=cell(1,numel(faceTypesNum));
C=(1:1:size(V,1))';
Cd=cell(1,numel(faceTypesNum));
for q=1:1:numel(faceTypesNum)
    %Get faces
    logicFacesNow=(n_sum==faceTypesNum(q)); %logic for current faces
    F_now=Fds(logicFacesNow,1:faceTypesNum(q)); %The current face set

    %Flip face orientation if required
    Nv_now=Nv(logicFacesNow,:); %Appropriate face normals based on input mesh
    [N_now]=patchNormal(F_now,Vd); %Current face normals
    logicFlip=dot(Nv_now,N_now,2)<0; %dot product is negative if face is inverted
    F_now(logicFlip,:)=fliplr(F_now(logicFlip,:)); %Flip faces

    indFaces=find(logicFacesNow);
    indFlip=indFaces(logicFlip);
    Fds(indFlip,1:faceTypesNum(q))=fliplr(Fds(indFlip,1:faceTypesNum(q)));

    Fd{q}=F_now; %Store faces in cell array
    Cd{q}=C(logicFacesNow);
end

if replaceVb==1 && fixBoundaryOption~=2
    Vd(end-numel(ind_Eb)+1:end,:)=V(ind_Eb,:);
end

varargout{1}=Vd; %Vertices
varargout{2}=Fd; %Cell array of faces
varargout{3}=Fds; %Array of face indices (with zeros)
varargout{4}=Cd; %Colors

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [ y, ierror ] = y_check_julian ( y )

%*****************************************************************************80
%
%% Y_CHECK_JULIAN checks a Julian year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, the year, which must not be 0.
%
%    Output, integer IERROR, is 0 if Y is legal, and 1 otherwise.
%
  if ( y ~= 0 )
    ierror = 0;
  else
    ierror = 1;
  end

  return
end


text:
function uinit = initSobol(u0,LB,UB,PLB,PUB,Ninit,optimState,options)
%INITSOBOL Initialize with space-filling design in the box

nvars = numel(u0);

try
    % Try to generate a quasirandom Sobol space-filling sequence

    MaxSeed = 997;
    if all(isfinite(u0))                        % Seed depends on U0
        strseed = num2str(u0(1:min(10,end)));
        seed = mod(prod(uint64(strseed)),MaxSeed)+1;
    else
        seed = randi(MaxSeed);                  % Random seed
    end
    r = i4_sobol_generate(nvars,Ninit,seed)';
    uinit = bsxfun(@plus, PLB, bsxfun(@times, r, PUB-PLB));
catch
    % Latin hypercube sampling if generation of Sobol sequence failed
    uinit = initLHS(u0,LB,UB,PLB,PUB,Ninit,optimState,options);
end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
%
% FILE NAME:    miscModelCurationScript_20190323.m
%
% PURPOSES: This script conducts a number of tasks:
%           1. Add Metadata to the annotation field of Human1
%           2. Reformat EC-number in eccodes field
%           3. Remove rxnComps field
%           4. Turn `version` into a blank field
%           5. Initialize rxnConfidenceScores field with zero
%

%% Load the model

if ~exist('ihuman','var')
    load('humanGEM.mat');  % version 1.0.0
end

%% Add Metadata to the annotation field of Human1

% Generate the Metadata as a structure
annotation.defaultLB=-1000;
annotation.defaultUB=1000;
annotation.taxonomy='9606';
annotation.note='Human genome-scale metabolic models are important tools for the study of human health and diseases, by providing a scaffold upon which different types of data can be analyzed. This is the latest version of human-GEM, which is a genome-scale model of the generic human cell. The objective of human-GEM is to serve as a community model for enabling integrative and mechanistic studies of human metabolism.';
annotation.sourceUrl='https://github.com/SysBioChalmers/human-GEM';
annotation.authorList='Jonathan Robinson, Hao Wang, Pierre-Etienne Cholley, Pınar Kocabaş';
annotation.email='nielsenj@chalmers.se';
annotation.organization='Chalmers University of Technology';

ihuman.annotation = annotation;
ihuman.description = 'Generic genome-scale metabolic model of Homo sapiens';

%% Reformat EC-number

% Multiple EC numbers should be separated by a semicolon and a blank space
% "; " according to RAVEN issue #184. Some incorrectly formatted eccodes
% elements, which are separated by " or " as reported in #93, are fixed here.
eccodes = regexprep(ihuman.eccodes,'\s*or\s*',';');

% Consistengly add a blank space after each semicolon
eccodes = regexprep(eccodes,';','; ');

% Update to the model
ihuman.eccodes = eccodes;

%% Emtpy version field

% Turn `version` into a blank field to retain a simple and clear work flow
ihuman.version = '';

%% Remove rxnComps field

% Take away the rxnComps field, according to RAVEN #184
ihuman = rmfield(ihuman, 'rxnComps');

%% Initialize rxnConfidenceScores field

% Assign elements in rxnConfidenceScores field with 0, based on #48
ihuman.rxnConfidenceScores(:) = 0;

%% Save the model files

writeHumanYaml(ihuman, 'humanGEM');
movefile('Human-GEM.yml','../../model/');
save('../../model/Human-GEM.mat','ihuman');

author: SysBioChalmers

repo: Human-GEM


text:
function value = r8_beta_sample ( aa, bb )

%*****************************************************************************80
%
%% R8_BETA_SAMPLE generates a beta random deviate.
%
%  Discussion:
%
%    This procedure returns a single random deviate from the beta distribution
%    with parameters A and B.  The density is
%
%      x^(a-1) * (1-x)^(b-1) / Beta(a,b) for 0 < x < 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2013
%
%  Author:
%
%    Original FORTRAN77 version by Barry Brown, James Lovato.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Russell Cheng,
%    Generating Beta Variates with Nonintegral Shape Parameters,
%    Communications of the ACM,
%    Volume 21, Number 4, April 1978, pages 317-322.
%
%  Parameters:
%
%    Input, real AA, the first parameter of the beta distribution.
%    0.0 < AA.
%
%    Input, real BB, the second parameter of the beta distribution.
%    0.0 < BB.
%
%    Output, real VALUE, a beta random variate.
%
  if ( aa <= 0.0 )
    fprintf ( 1, ' \n' );
    fprintf ( 1, 'R8_BETA_SAMPLE - Fatal error!\n' );
    fprintf ( 1, '  AA <= 0.0\n' );
    error ( 'R8_BETA_SAMPLE - Fatal error!\n' );
  end

  if ( bb <= 0.0 )
    fprintf ( 1, ' \n' );
    fprintf ( 1, 'R8_BETA_SAMPLE - Fatal error!\n' );
    fprintf ( 1, '  BB <= 0.0\n' );
    error ( 'R8_BETA_SAMPLE - Fatal error!\n' );
  end
%
%  Algorithm BB
%
  if ( 1.0 < aa && 1.0 < bb )

    a = min ( aa, bb );
    b = max ( aa, bb );
    alpha = a + b;
    beta = sqrt ( ( alpha - 2.0 ) / ( 2.0 * a * b - alpha ) );
    gamma = a + 1.0 / beta;

    while ( 1 )

      u1 = r8_uniform_01_sample ( );
      u2 = r8_uniform_01_sample ( );
      v = beta * log ( u1 / ( 1.0 - u1 ) );
      w = a * exp ( v );

      z = u1 ^ 2 * u2;
      r = gamma * v - log ( 4.0 );
      s = a + r - w;

      if ( 5.0 * z <= s + 1.0 + log ( 5.0 ) )
        break
      end

      t = log ( z );
      if ( t <= s )
        break
      end

      if ( t <= ( r + alpha * log ( alpha / ( b + w ) ) ) )
        break
      end

    end
%
%  Algorithm BC
%
  else

    a = max ( aa, bb );
    b = min ( aa, bb );
    alpha = a + b;
    beta = 1.0 / b;
    delta = 1.0 + a - b;

    k1 = delta * ( 1.0 / 72.0 + b / 24.0 ) / ( a / b - 7.0 / 9.0 );
    k2 = 0.25 + ( 0.5 + 0.25 / delta ) * b;

    while ( 1 )

      u1 = r8_uniform_01_sample ( );
      u2 = r8_uniform_01_sample ( );

      if ( u1 < 0.5 )

        y = u1 * u2;
        z = u1 * y;

        if ( k1 <= 0.25 * u2 + z - y )
          continue
        end

      else

        z = u1 ^ 2 * u2;

        if ( z <= 0.25 )

          v = beta * log ( u1 / ( 1.0 - u1 ) );
          w = a * exp ( v );

          if ( aa == a )
            value = w / ( b + w );
          else
            value = b / ( b + w );
          end

          return

        end

        if ( k2 < z )
          continue
        end

      end

      v = beta * log ( u1 / ( 1.0 - u1 ) );
      w = a * exp ( v );

      if ( log ( z ) <= ...
        alpha * ( log ( alpha / ( b + w ) ) + v ) - log ( 4.0 ) )
        break
      end

    end

  end

  if ( aa == a )
    value = w / ( b + w );
  else
    value = b / ( b + w );
  end

  return
end


text:
function trapezoidalODE(positions, velocities, firstOrder, dt, nSteps, edges)

f = figure;

% Want to solve dx/dt = firstOrder * x
x = [positions ; velocities];
n = size(positions,1);

for i=1:nSteps
    % TODO:  Step x forward in time here!
    x = inv(eye(2 * n) - dt * firstOrder / 2) * (eye(2 * n) + dt * firstOrder / 2) * x;

    positions = x(1:n,:);
    plotGraph(positions,edges,f);
    title(sprintf('Trapezoidal step %d',i));
    drawnow;
end

author: Doraemonzzz


text:
function y = Fault_decision_tau(u)
global st13

% Using model defined by the structure st13 to make decision about fault in the couple tau

u(3)=1e-10*u(3).*(u(end).^6);u=u(1:3);
sker=st13.x2sup+(abs(u))'*abs(u)*ones(st13.Nlsup,1)-2*st13.xsup*abs(u);
y=(st13.w)'*exp(-sker./(2*(st13.sigma).^2))+st13.b;

if y>=0.8
    y=1;
else
    y=0;
end
end


text:
function figbolden(pos,pv,m)
%FIGBOLDEN embolden, resize and recolour the current figure =(POS,PV,M)
%
% Inputs: pos = [xmin ymin width height] gives the lower left corner position and the window size in pixels
%               [width height] leaves the lower left corner alone
%               [width] has a standard aspect ratio of 4:3
%               [-width/height] leaves the area unchanged but fixes the aspect ratio
%         pv is a cell name containing attribute-value pairs.
%            default = {'FontName' 'Ariel'; 'FontSize' 16; 'LineWidth' 2; 'MarkerSize' 8}
%         m is a mode string:
%                'l' lists the changes made
%                'd' uses default pv settings (use [] for second argument)
%                'c' change default colours to improve contrast on a white background
%                    g->[0,0.7,0],c->[0,0.7,0.7],y->[0.83,0.83,0]
%                'x' suppresses all changes
%
% Bug: gives an error message if log axes have been used

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

ps={'Title' 'XLabel' 'YLabel' 'Children'};
if nargin<3
    m='';
end
if nargin<2 || any(m=='d')
    pv={'FontName' 'Arial'; 'FontSize' 16; 'LineWidth' 2; 'MarkerSize' 8};
end
pp={'Symbol';'Wingdings'};      % protected fonts
if nargin<1
    pos=[];
end
mlist=any(m=='l');  % list changes
mnotx=~any(m=='x'); % do changes
scsz=get(0,'screensize');
if length(pos)
    po=get(gcf,'position');
    if length(pos)>2            % position is specified
        po(1:2)=pos(1:2);
        pos(1:2)=[];      % remove xmin,ymin
    end
    if length(pos)>1
        po(3:4)=pos(1:2);
    else
        if pos(1)>0
            po(3:4)=[1 0.75]*pos(1);
        else
            po(3:4)=[-pos(1) 1]*sqrt(-po(3)*po(4)/pos(1)); % preserve area
        end
    end
    set(gcf,'position',po);
end
if any(m=='c')
    cc='gcy';
    cv=[0,0.7,0;0,0.7,0.7;0.83,0.83,0];
    for ic=1:3
        hlist=findobj(gcf,'color',cc(ic));
        if length(hlist)
            if mnotx
                set(hlist,'color',cv(ic,:));
            end
            if mlist
            for ih=1:length(hlist)
                 fprintf(['change %f Color: ''%c'' -> [%g %g %g]\n'],hlist(ih),cc(ic),cv(ic,:));
            end
            end
        end
    end
end
hlist=get(gcf,'children');
while length(hlist)
    pl=get(hlist(1));
    %fprintf('list length = %d, handle = %f\n',length(hlist),hlist(1));
    for i=1:size(pv,1)
        if isfield(pl,pv{i,1})
            if i>1 || all(~strcmpi(get(hlist(1),pv{i,1}),pp))
                pval=get(hlist(1),pv{i,1});
                if ~all(size(pval)==size(pv{i,2})) || ~all(pval(:) == pv{i,2}(:))
                    if mnotx
                        set(hlist(1),pv{i,1},pv{i,2})
                    end
                    if mlist
                        if ischar(pval)
                            fch='s';
                        else
                            fch='g';
                        end
                        fprintf(['change %f %s: %' fch ' -> %' fch '\n'],hlist(1),pv{i,1},pval,pv{i,2});
                    end
                end
            end
        end
    end
    for i=1:length(ps)
        if isfield(pl,ps{i})
            hlist=[hlist; get(hlist(1),ps{i})];
            %fprintf('add %f:%s\n',hlist(1),ps{i});
        end
    end
    hlist(1)=[];
end

repo: ekfmonoslam


text:
%kicomposite 'Composites two images '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros icomposite.pane file
%
% Parameters:
% MultiChoice: toggle 'Flag', default: 1: 'Flag toggle selection'
%    Choices are:
%   1: 'aPLUSb'
%   2: 'aOVERb'
%   3: 'aOVERb'
%   4: 'aINb'
%   5: 'bINa'
%   6: 'aOUTb'
%   7: 'bOUTa'
%   8: 'aATOPb'
%   9: 'bATOPa'
%   10: 'aXORb'
% InputFile: i1 'Input 1', required: 'First Input data object'
% InputFile: i2 'Input 2', required: 'input file'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kicomposite({i1, i2}, {'toggle',1;'i1','';'i2','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% icomposite - no short description specified
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [logp, yhat, res] = tapas_unitsq_sgm(r, infStates, ptrans)
% Calculates the log-probability of response y=1 under the unit-square sigmoid model
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Transform zeta to its native space
ze = exp(ptrans(1));

% Initialize returned log-probabilities as NaNs so that NaN is
% returned for all irregualar trials
n = size(infStates,1);
logp = NaN(n,1);
yhat = NaN(n,1);
res  = NaN(n,1);

% Weed irregular trials out from inferred states and responses
x = infStates(:,1,1);
x(r.irr) = [];
y = r.y(:,1);
y(r.irr) = [];

% Avoid any numerical problems when taking logarithms close to 1
logx = log(x);
log1pxm1 = log1p(x-1);
logx(1-x<1e-4) = log1pxm1(1-x<1e-4);
log1mx = log(1-x);
log1pmx = log1p(-x);
log1mx(x<1e-4) = log1pmx(x<1e-4);

% Calculate log-probabilities for non-irregular trials
reg = ~ismember(1:n,r.irr);
logp(reg) = y.*ze.*(logx -log1mx) +ze.*log1mx -log((1-x).^ze +x.^ze);
yhat(reg) = x;
res(reg) = (y-x)./sqrt(x.*(1-x));

return;

author: translationalneuromodeling

repo: tapas


text:
function P_hat = simpleEVD(X, r)

%%%Function to implement the Simple-EVD algorithm and returns a basis
%%%matrix for the new subspace -- edit these comments before Gitting

%%%                         Inputs                          %%%
%%%         X - Emperical Covariance data matrix (m X m)    %%%
%%%         r - Target rank of output                       %%%

%%%                         Outputs                         %%%
%%%         P_hat - Basis matrix for output (m X r)         %%%

[P_hat, ~] = svds(X, r);
end


text:
function version = gpstuff_version(str)
% Return version of currently used gpstuff
%
% gpstuff_version() returns string including release number and release
%                   date
%
% gpstuff_version('date') returns the release date (string)
%
% gpstuff_version('release') returns the release number (double)
%
% Copyright (c) 2013 Ville Tolvanen

if nargin==0
  version='4.2_2013-06-14';
elseif isequal(str, 'date')
  version='2013-06-14';
elseif isequal(str, 'release')
  version=4.2;
end

end

author: gpstuff-dev

repo: gpstuff

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


text:
classdef SQP < ALGORITHM
% <single> <real> <large/none> <constrained/none>
% Sequential quadratic programming

%------------------------------- Reference --------------------------------
% P. T. Boggs and J. W. Tolle, Sequential quadratic programming, Acta
% Numerica, 1995, 4(1): 1-51.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate random solution
            X     = Problem.Initialization(1);
            ro    = 0.5;
            eta   = 0.1;
            lam   = zeros(1,length(X.con));
            Bk    = eye(Problem.D);
            sigma = 0.8;
            dfk   = Problem.CalObjGrad(X.dec);
            Ai    = Problem.CalConGrad(X.dec);

            %% Optimization
            while Algorithm.NotTerminated(X)
                [~,dk,mu,~] = qpsubp(dfk',Bk,[],[],-Ai,-X.con');
                tau = max(norm(mu,inf),norm(lam,inf));
                if sigma*(tau+0.05) >= 1
                    sigma = 1/(tau+2*0.05);
                end
                for mk = 0 : 20
                    temp = eta*ro^mk*dphi1(X,sigma,dfk,dk);
                    X1   = Problem.Evaluation(X.dec+ro^mk*dk');
                    if phi1(X1,sigma)-phi1(X,sigma) < temp
                        break;
                    end
                end
                [dfk0,Ai0] = deal(dfk,Ai);
                dfk        = Problem.CalObjGrad(X1.dec);
                Ai         = Problem.CalConGrad(X1.dec);
                lam = pinv(-Ai')*dfk';
                sk  = (X1.dec-X.dec)';
                yk  = dlax(dfk,-Ai,lam) - dlax(dfk0,-Ai0,lam);
                if sk'*yk>0.2*sk'*Bk*sk
                    omega = 1;
                else
                    omega = 0.8*sk'*Bk*sk/(sk'*Bk*sk-sk'*yk);
                end
                zk = omega*yk + (1-omega)*Bk*sk;
                Bk = Bk + zk*zk'/(sk'*zk) - (Bk*sk)*(Bk*sk)'/(sk'*Bk*sk);
                X  = X1;
            end
        end
    end
end

function p = phi1(X,sigma)
    p = X.obj + 1/sigma*norm(max(X.con',0),1);
end

function dp = dphi1(X,sigma,df,d)
    dp = df*d - 1/sigma*norm(max(X.con',0),1);
end

function dl = dlax(df,Ai,lam)
    dl = df' - Ai'*lam;
end

author: BIMK

repo: PlatEMO


text:
% test for image transformation

name = 'barb';
n = 128;

M = load_image(name, n);

% find sampling location for an image
[Y,X] = meshgrid(1:n,1:n);
pos = [X(:),Y(:)]';

% perform the transform
alpert_vm = [2 2];
% this is very important for images
options.part_type = '2axis';
w = perform_alpert_transform_2d(M(:),pos,alpert_vm, 1, options);

% perform some thresholding
T = 20; % this threshold is ok for images in range [1...256], otherwise scale it
w1 = w .* (abs(w)>T);

% reconstruction
v1 = perform_alpert_transform_2d(w1,pos,alpert_vm, -1, options);
M1 = reshape(v1,n,n);

clf;
subplot(1,2,1);
imagesc(M);
title('Original'); axis image; axis off;
subplot(1,2,2);
imagesc(M1);
title('Approximated'); axis image; axis off;
colormap gray(256);

author: gpeyre

repo: matlab-toolboxes


text:
function [voice accom midi] = main(x, config, midi)
%MAIN   Main function for voice/accompaniment separation.
%   Format: [voice accom midi] = main(x, config, midi)
%   Inputs:
%       x:      Mixed signal. Must be mono-channel, sampled at 16 kHz.
%       config (optional):
%               Struct of configurations. If omitted or empty, default
%               settings will be used.
%       midi (optional):
%               Referential pitch contour, in midi number scale, one number
%               per frame (40 ms frame length, 20 ms frame shift).
%               Non-positive frames will be totally assigned to the
%               accompaniment.
%   Outputs:
%       voice:  Separated voice signal.
%       accom:  Separated accompaniment signal.
%       midi (optional):
%           Extracted or provided pitch contour. Positive numbers are f0's
%           in midi numbers, 0's stand for accompaniment, and -1's stand
%           for unvoiced frames.

    if (nargin < 2 || isempty(config))
        config.configMfcc = configMfcc;
        config.configEsi = configEsi;
        config.configSep = configSep;
        load HMM
        config.mfccHmm = mfccHmm0;
        config.esiHmm = esiHmm0;
    end

    % Replace zeros with very small noise to avoid nasty NaN problems
    zeroInd = find(x == 0);
    x(zeroInd) = randn(length(zeroInd), 1) * 1e-10;

    if (nargin < 3)
        midi = extractPitch(config.configMfcc, config.configEsi, x, config.mfccHmm, config.esiHmm);
    end
    [voice accom] = separate(config.configSep, x, midi);
end

author: MaigoAkisame


text:
function [out, T] = BF_ParEval(x,x_z,s,beVocal)
% BF_ParEval    Sneaky eval function for parfor loops
%
% parfor loops don't allow the 'eval' function.
% This sneaky sneaky gets around that.
%
%---INPUTS:
% x and y are possible elements of the string s to be evaluated
% Stores any text output as T

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

if beVocal
    % Any text output from operations is printed to screen
    out = eval(s);
else
    % Text output from operations is suppressed and stored in T
    [T, out] = evalc(s);
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [data_mod, index_rem]=feat_sel_sim(data, measure, p)
%
% Feature selection method using similarity measure and fuzzy entroropy
% measures based on the article:
%
% P. Luukka, (2011) Feature Selection Using Fuzzy Entropy Measures with
% Similarity Classifier, Expert Systems with Applications, 38, pp.
% 4600-4607
%
% Function call:
% [data_mod, index_rem]=feat_sel_sim(data, measure, p)
%
% OUTPUTS:
% data_mod      data without removed feature
% index_rem     index of removed feature in original data
%
%INPUTS:
% data          data matrix, contains class values
% measure       fuzzy entropy measure, either 'luca' or 'park'
%               currently coded
% p             parameter of Lukasiewicz similarity measure
%               p in (0, \infty) as default p=1.
if nargin<3
   p=1;
end
if nargin<2
    measure='luca'
end

l=max(data(:,end));     % #-classes
m=size(data,1);         % #-samples
t=size(data,2)-1;       % #-features
dataold=data;
tmp=[];
% forming idealvec using arithmetic mean
idealvec=zeros(l,t);
for k=1:l
    idealvec_s(k,:)=mean(data(find(data(:,end)==k),1:t));
end

%scaling data between [0,1]
data_v=data(:,1:t);
data_c=data(:,t+1); %labels
mins_v = min(data_v);
Ones = ones(size(data_v));
data_v = data_v+Ones*diag(abs(mins_v));
for k=1:l
    tmp=[tmp;abs(mins_v)];
end
tmp;
idealvec_s = idealvec_s+tmp;
maxs_v = max(data_v);
data_v = data_v*diag(maxs_v.^(-1));
idealvec_s=idealvec_s./repmat(maxs_v,l,1);
data = [data_v, data_c];
% sample data
datalearn_s=data(:,1:t);

% similarities
sim=zeros(t,m,l);
for j=1:m
    for i=1:t
        for k=1:l
            sim(i,j,k)=(1-abs(idealvec_s(k,i)^p-datalearn_s(j,i))^p)^(1/p);
        end
    end
end

% reduce number of dimensions in sim
sim=reshape(sim,t,m*l)';

% possibility for two different entropy measures
if measure=='luca'
% moodifying zero and one values of the similarity values to work with
% De Luca's entropy measure
    delta=1E-10;
    sim(find(sim==0))=delta;
    sim(find(sim==1))=1-delta;
    H=sum(-sim.*log(sim)-(1-sim).*log(1-sim));

elseif measure=='park'
    H=sum(sin(pi/2*sim)+sin(pi/2*(1-sim))-1);
end

% find maximum feature
[i, index_rem]=max(H);

% removing feature from the data
data_mod=[dataold(:,1:index_rem-1) dataold(:,index_rem+1:end)];
j=1;

return

author: Sable

repo: mcbench-benchmarks


text:
function [k, sk, n2] = rbfKernCompute(kern, x, x2)

% RBFKERNCOMPUTE Compute the RBF kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the radial basis function
% kernel given inputs associated with rows and columns.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 : the input matrix associated with the columns of the kernel.
% RETURN k : the kernel matrix computed at the given points.
% RETURN sk : unscaled kernel matrix (i.e. only the exponential part).
%
% FORMAT
% DESC computes the kernel matrix for the radial basis function
% kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : the kernel matrix computed at the given points.
% RETURN sk : unscaled kernel matrix (i.e. only the exponential part).
%
% SEEALSO : rbfKernParamInit, kernCompute, kernCreate, rbfKernDiagCompute
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006
%
% MODIFICATIONS : Mauricio Alvarez, 2009, David Luengo, 2009

% KERN

if nargin < 3
  n2 = dist2(x, x);
else
  n2 = dist2(x, x2);
end

wi2 = (.5 .* kern.inverseWidth);
sk = exp(-n2*wi2);
k = kern.variance*sk;
if isfield(kern, 'isNormalised') && (kern.isNormalised == true)
    k = k * sqrt(kern.inverseWidth/(2*pi));
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [xy, idx] = paduapts(n, dom)
%PADUAPTS   Padua points.
%   XY = PADUAPTS(N) returns a two-column vector containing the x and y
%   co-ordinates of the degree N first-kind Padua points on [-1 1] x [-1,1].
%
%   XY = PADUAPTS(N, [a, b, c, d]) returns the degree N Padua points on the
%   domain [a b] x [c d].
%
%   The ordering and defniiton of the points is consistent with Padua2DM [1].
%
%   [XY, IDX] = PADUAPTS(...) returns also a logical matrix IDX which denotes
%   the entries of the Chebyshev (N+1)x(N+2) tensor product grid that form XY.
%
% See also CHEBPTS.

% References:
%   [1]: Marco Caliari, Stefano De Marchi, Alvise Sommariva, Marco Vianello
%        "Padua2DM: fast interpolation and cubature at the Padua points in
%         Matlab/Octave."

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

% Default domain:
if ( nargin < 2 )
    dom = [-1, 1, -1, 1];
end

% Trivial case:
if ( n == 0 )
    xy = dom([1, 3]);
    return
end

% 1D Chebyshev grids:
xn1 = chebpts(n+1, dom(1:2), 2);
xn2 = chebpts(n+2, dom(3:4), 2);

% Order is flipped from Chebfun for consistency with [1]:
xn1 = xn1(end:-1:1);
xn2 = xn2(end:-1:1);

% Full tensor grid.
[x1, x2] = meshgrid(xn1, xn2);

% Extract every other term.
idx = true(1, (n+1)*(n+2));
idx(1:2:end) = false;
if ( mod(n, 2) )
    idx = reshape(idx, n+2, n+1);
else
    idx = reshape(idx, n+1, n+2).';
end

% Padua points.
xy = [x1(idx), x2(idx)];

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ Pt ] = transformPoints( P, bb, bw )
%TRANSFORMPOINTS Transforms coordinates of detected landmarks into the
%normalized image frame.
%   Detailed explanation goes here
%
% INPUT:
%   P       ...     [2 x M (double)] Detected landmarks coordinates in
%                   normalized image frame. M is the number of components
%   bb      ...     [1 x 4 (double)] bbox of normalized image frame in
%                   original image (i.e. the extended bbox of face
%                   detection)
%   bw      ...     [2 x 1 (double)] base window size
%
% OUTPUT:
%   Pt      ...     [2 x M (double)] Tranformed coordinates of detected
%                   landmarks
%
% 05-08-11 Michal Uricar

    Origin = [bb(1); bb(2)];

    % compute the scale factor
    scalefac = [ (bb(3)-bb(1)+1)/bw(1);
                 (bb(4)-bb(2)+1)/bw(2) ];

    M = size(P, 2);

    Pt = P .* repmat(scalefac, 1, M) + repmat(Origin, 1, M);

end

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function x_store = my_enframe(x, frame_size, frame_shift, useLastPartialFrame)
if nargin<4
    useLastPartialFrame = 0;    % whether to use the last few samples that are not as long as a frame
    % if set to >0, it represents the minimum percentage of a frame
    % required to generate the last frame
end

[nSample, nCh] = size(x);

% directly create an array of desired format. This will make the call to
% enframe_core faster
if nSample < 2^16
    frame_idx = (uint16(1) : uint16(nSample))';
elseif nSample<2^32
    frame_idx = (uint32(1) : uint32(nSample))';
else
    frame_idx = (1 : nSample)';
end

idx = enframe_core(frame_idx,frame_size,frame_shift, useLastPartialFrame);  % get the index of samples in each frame

if isempty(idx)
    x_store = [];
else
    if useLastPartialFrame>0 && idx(end,end)==0          % if there is 0 in the index, the last index must be zero.
        idx2 = max(1,idx);
        x_store = x(idx2,:);            % index all channels simultaneously for speed
        x_store(idx==0,:) = 0;          % sometimes, the last frames are padded with 0
    else
        x_store = x(idx,:);
    end
    x_store = reshape(x_store, frame_size, size(idx,2), nCh);
    x_store = permute(x_store, [1 3 2]);
end
end

author: singaxiong

repo: SignalGraph


text:
function H = hankel(c,r)

%Disciplined convex/geometric programming information for HANKEL:
%   HANKEL imposes no convexity restrictions on its arguments. Instead
%   of using the HANKEL function, however, consider creating a matrix
%   variable using the 'hankel' or 'upper_hankel' keyword; e.g.
%       variable X(5,5) hankel;
%       variable Y(4,4) upper_hankel;

%
% Check arguments
%

if nargin < 2,
    r = zeros(size(c));
else
    temp = cvx_subsref( r, 1 ) - cvx_subsref( c, numel(c) );
    if ~cvx_isnonzero( temp ),
        warning('MATLAB:hankel:AntiDiagonalConflict',['Last element of ' ...
               'input column does not match first element of input row. ' ...
               '\n         Column wins anti-diagonal conflict.'])
    end
end

%
% Compute indices and construct data vector
%

r  = vec( r );
c  = vec( c );
nc = length( c );
nr = length( r );
x  = [ c ; cvx_subsref( r, 2 : nr, 1 ) ];

%
% Construct matrix
%

cidx = ( 1 : nc )';
ridx = 0 : nr - 1;
H    = cidx(:,ones(nr,1)) + ridx(ones(nc,1),:);
H    = reshape( cvx_subsref( x, H( : ) ), size( H ) );

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

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook