text:
function varargout = process_segment_fastsurfer( varargin )
% PROCESS_SEGMENT_FASTSURFER: Run the segmentation of a T1 MRI with FastSurfer.
%
% USAGE:     OutputFiles = process_segment_fastsurfer('Run',     sProcess, sInputs)
%         [isOk, errMsg] = process_segment_fastsurfer('Compute', iSubject, iAnatomy=[default], nVertices, isInteractive)
%                          process_segment_fastsurfer('ComputeInteractive', iSubject, iAnatomy)

% @=============================================================================
% 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, 2021-2023

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Segment MRI with FastSurfer';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Import', 'Import anatomy'};
    sProcess.Index       = 33;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/SegFastSurfer';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'import'};
    sProcess.OutputTypes = {'import'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 0;
    % Title
    sProcess.options.install.Comment = [
        'FastSurfer must be installed on the computer, the command
' ...
        'run_fastsurfer.sh must be available in the system path.
' ...
        'Click "Online tutorial" for installation instructions.'];
    sProcess.options.install.Type    = 'label';
    % Doc
    sProcess.options.optdoc.Comment = [
        '
' ...
        'Network specific arguments:
' ...
        '--clean_seg : Flag to clean up FastSurferCNN segmentation
' ...
        '--no_cuda   : Flag to disable CUDA usage in FastSurferCNN
' ...
        '--batch     : Batch size for inference. Default: 16.
' ...
        '              Lower this to reduce memory requirement.
' ...
        '--order     : Interp. for mri_convert T1 before segmentation
' ...
        '              (0=nearest, 1=linear(default), 2=quadratic, 3=cubic)

' ...
        'Surface pipeline arguments:
' ...
        '--fstess    : Use mri_tesselate instead of marching cube (default)
' ...
        '--fsqsphere : Use FreeSurfer default instead of novel
' ...
        '              spectral spherical projection for qsphere
' ...
        '--fsaparc   : Use FS aparc segmentations + DL prediction
' ...
        '--surfreg   : sphere.reg registration with FreeSurfer
' ...
        '--parallel  : Run both hemispheres in parallel
' ...
        '--threads   : Set openMP and ITK threads to ...

' ...
        'Other:
' ...
        '--py        : Python version to use. Default: python3.6
' ...
        '--seg_only  : Only run FastSurferCNN
' ...
        '--surf_only : Only run the surface pipeline recon_surf

'];
    sProcess.options.optdoc.Type    = 'label';
    % Option: Extra command-line parameters
    sProcess.options.param.Comment = 'Command-line options:';
    sProcess.options.param.Type    = 'text';
    sProcess.options.param.Value   = '--batch 4 --surfreg --parallel --threads 4';
    % Option: Subject name
    sProcess.options.subjectname.Comment = 'Subject name:';
    sProcess.options.subjectname.Type    = 'subjectname';
    sProcess.options.subjectname.Value   = '';
    % File selection options
    SelectOptions = {...
        '', ...                            % Filename
        '', ...                            % FileFormat
        'open', ...                        % Dialog type: {open,save}
        'Select FastSurfer output folder...', ...     % Window title
        'ExportAnat', ...                  % LastUsedDir: {ImportData,ImportChannel,ImportAnat,ExportChannel,ExportData,ExportAnat,ExportProtocol,ExportImage,ExportScript}
        'single', ...                      % Selection mode: {single,multiple}
        'dirs', ...                        % Selection mode: {files,dirs,files_and_dirs}
        {{'.folder'}, 'FastSurfer database folder', 'FastSurfer'}, ... % Available file formats
        []};                               % DefaultFormats: {ChannelIn,DataIn,DipolesIn,EventsIn,AnatIn,MriIn,NoiseCovIn,ResultsIn,SspIn,SurfaceIn,TimefreqIn}
    % Option: FastSurfer output folder
    sProcess.options.outdir.Comment = 'Output folder:';
    sProcess.options.outdir.Type    = 'filename';
    sProcess.options.outdir.Value   = SelectOptions;
    % Default outdir
    sProcess.options.defdir.Comment = 'Default output folder: $HOME/.brainstorm/tmp/fastsurfer';
    sProcess.options.defdir.Type    = 'label';
    % Option: Number of vertices
    sProcess.options.nvertices.Comment = 'Number of vertices (cortex): ';
    sProcess.options.nvertices.Type    = 'value';
    sProcess.options.nvertices.Value   = {15000, '', 0};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    OutputFiles = {};
    % Number of vertices
    nVertices = sProcess.options.nvertices.Value{1};
    if isempty(nVertices) || (nVertices < 50)
        bst_report('Error', sProcess, [], 'Invalid number of vertices.');
        return
    end
    % Get subject name
    SubjectName = file_standardize(sProcess.options.subjectname.Value);
    if isempty(SubjectName)
        bst_report('Error', sProcess, [], 'Subject name is empty.');
        return;
    end
    % Get subject
    [sSubject, iSubject] = bst_get('Subject', SubjectName);
    if isempty(iSubject)
        bst_report('Error', sProcess, [], ['Subject "' SubjectName '" does not exist.']);
        return
    end
    % Get extra parameters
    param = sProcess.options.param.Value;
    outdir = sProcess.options.outdir.Value{1};
    % Call processing function
    [isOk, errMsg] = Compute(iSubject, [], nVertices, 0, param, outdir);
    % Handling errors
    if ~isOk
        bst_report('Error', sProcess, [], errMsg);
    elseif ~isempty(errMsg)
        bst_report('Warning', sProcess, [], errMsg);
    end
    % Return an empty structure
    OutputFiles = {'import'};
end

%% ===== COMPUTE SEGMENTATION =====
function [isOk, errMsg] = Compute(iSubject, iAnatomy, nVertices, isInteractive, param, outdir)
    errMsg = '';
    isOk = 0;

    % ===== INSTALL FASTSURFER =====
    bst_progress('text', 'Testing FastSurfer installation...');
    % Get FreeSurfer path from environment variable
    FreeSurferDir = getenv('FREESURFER_HOME');
    SubjectsDir = getenv('SUBJECTS_DIR');
    if isempty(FreeSurferDir) || ~file_exist(FreeSurferDir) || isempty(SubjectsDir) || ~file_exist(SubjectsDir)
        errMsg = ['FreeSurfer is not installed on your computer.' 10 ...
                  'The environment variables FREESURFER_HOME and SUBJECTS_DIR must be set.' 10 10 ...
                  'See the online Brainstorm tutorial about FastSurfer for help:' 10 ...
                  'https://neuroimage.usc.edu/brainstorm/Tutorials/SegFastSurfer'];
        return
    end
    % Get FastSurfer path with system call
    [res,FastSurferDir] = system('which run_fastsurfer.sh');
    if (res ~= 0) || isempty(FastSurferDir)
        errMsg = ['FastSurfer is not installed on your computer.' 10 ...
                  'The executable run_fastsurfer.sh must be available in the system path.' 10 ...
                  'Click on "Help" in the process options for installation instructions.'];
        return
    end
    FastSurferDir = bst_fileparts(strtrim(FastSurferDir));

    % ===== GET SUBJECT =====
    % Get subject
    [sSubject, iSubject] = bst_get('Subject', iSubject);
    if isempty(sSubject)
        errMsg = 'Subject does not exist.';
        return
    end
    % Check if a MRI is available for the subject
    if isempty(sSubject.Anatomy)
        errMsg = ['No MRI available for subject "' SubjectName '".'];
        return
    end
    % Get default MRI if not specified
    if isempty(iAnatomy)
        iAnatomy = sSubject.iAnatomy;
    end

    % ===== CHECK OUTPUT FOLDER =====
    % Folder for FastSurfer output
    if ~isempty(outdir)
        TmpDir = [];
        procDir = outdir;
    else
        TmpDir = bst_get('BrainstormTmpDir', 0, 'fastsurfer');
        procDir = TmpDir;
    end
    if ~file_exist(procDir)
        mkdir(procDir);
    end
    % Check if subject dir already exists
    subjid = strrep(sSubject.Name, '@', '');
    subjDir = bst_fullfile(procDir, subjid);
    if file_exist(subjDir)
        if isInteractive
            isDelete = java_dialog('confirm', ...
                    ['Subject folder already exists:
' subjDir '

' ...
                    'Delete existing folder?

'], 'FastSurfer segmentation');
            if isDelete
                file_delete(subjDir, 1, 3);
            end
        else
            bst_report('Warning', 'process_segment_fastsurfer', [], ['Subject folder already exists: ' subjDir]);
        end
    end

    % ===== DELETE EXISTING SURFACES =====
    % Confirm with user that the existing surfaces will be removed
    if isInteractive && ~isempty(sSubject.Surface)
        isDel = java_dialog('confirm', ['Warning: There are already surfaces in this subject.' 10 ...
            'Running FastSurfer will remove all the existing surfaces.' 10 10 ...
            'Delete the existing files?' 10 10], 'FastSurfer segmentation');
        if ~isDel
            errMsg = 'Process aborted by user.';
            % Delete temporary files
            if ~isempty(TmpDir)
                file_delete(TmpDir, 1, 1);
            end
            return;
        end
    end

    % ===== VERIFY FIDUCIALS IN MRI =====
    bst_progress('text', 'Saving temporary files...');
    % Load MRI file
    T1FileBst = sSubject.Anatomy(iAnatomy).FileName;
    sMri = in_mri_bst(T1FileBst);
    % If the SCS transformation is not defined: compute MNI transformation to get a default one
    if isempty(sMri) || ~isfield(sMri, 'SCS') || ~isfield(sMri.SCS, 'NAS') || ~isfield(sMri.SCS, 'LPA') || ~isfield(sMri.SCS, 'RPA') || (length(sMri.SCS.NAS)~=3) || (length(sMri.SCS.LPA)~=3) || (length(sMri.SCS.RPA)~=3) || ~isfield(sMri.SCS, 'R') || isempty(sMri.SCS.R) || ~isfield(sMri.SCS, 'T') || isempty(sMri.SCS.T)
        % Issue warning
        errMsg = 'Missing NAS/LPA/RPA: Computing the MNI normalization to get default positions.';
        % Compute MNI normalization
        [sMri, errNorm] = bst_normalize_mni(T1FileBst);
        % Handle errors
        if ~isempty(errNorm)
            errMsg = [errMsg 10 'Error trying to compute the MNI normalization: ' 10 errNorm 10 ...
                'Missing fiducials: the surfaces cannot be aligned with the MRI.'];
        end
    end

    % ===== SAVE MRI AS NII =====
    % Save MRI in .nii format
    NiiFile = bst_fullfile(procDir, [subjid, '.nii']);
    out_mri_nii(sMri, NiiFile);
    % If a "world transformation" was not available in the MRI in the database, it was set to a default when saving to .nii
    % Let's reload this file to get the transformation matrix, it will be used when importing the results
    if ~isfield(sMri, 'InitTransf') || isempty(sMri.InitTransf) || isempty(find(strcmpi(sMri.InitTransf(:,1), 'vox2ras')))
        % Load again the file, with the default vox2ras transformation
        [tmp, vox2ras] = in_mri_nii(NiiFile, 0, 0, 0);
        % Prepare the history of transformations
        if ~isfield(sMri, 'InitTransf') || isempty(sMri.InitTransf)
            sMri.InitTransf = cell(0,2);
        end
        % Add this transformation in the MRI
        sMri.InitTransf(end+1,[1 2]) = {'vox2ras', vox2ras};
        % Save modification on hard drive
        bst_save(file_fullpath(T1FileBst), sMri, 'v7');
    end

    % ===== RUN FASTSURFER =====
    bst_progress('text', 'Running FastSurfer...    (see command window)');
    strCall = ['cd ' FastSurferDir '; run_fastsurfer.sh --t1 "' NiiFile '" --sid "' subjid '" --sd "' procDir '" ' param];
    disp(['BST> System call: ' strCall]);
    status = system(strCall)
    % Error handling
    if (status ~= 0)
        errMsg = ['FastSurfer failed.', 10, 'Check the Matlab command window for more information.'];
        return
    end

    % ===== IMPORT OUTPUT FOLDER =====
    % Import FreeSurfer anatomy folder
    isExtraMaps = 0;
    isKeepMri = 1;
    isVolumeAtlas = 1;
    FsDir = bst_fullfile(procDir, subjid);
    errMsg = import_anatomy_fs(iSubject, FsDir, nVertices, 0, [], isExtraMaps, isVolumeAtlas, isKeepMri);
    if ~isempty(errMsg)
        return;
    end
    % Delete temporary files
    if ~isempty(TmpDir)
        file_delete(TmpDir, 1, 1);
    end
    % Return success
    isOk = 1;
end

%% ===== COMPUTE/INTERACTIVE =====
function ComputeInteractive(iSubject, iAnatomy) %#ok<DEFNU>
    % Get inputs
    if (nargin < 2) || isempty(iAnatomy)
        iAnatomy = [];
    end
    % Get default parameters
    sProcess = GetDescription();
    optdoc = sProcess.options.optdoc.Comment;
    defParam = sProcess.options.param.Value;
    defDir = bst_fullfile(bst_get('BrainstormTmpDir'), 'fastsurfer');
    % Ask for extra processing parameters
    [param, isCancel] = java_dialog('input', optdoc, 'FastSurfer command-line options', [], defParam);
    if isCancel
        return;
    end
    % Ask for output folder
    [outdir, isCancel] = java_dialog('input', 'FastSurfer output folder:', 'FastSurfer MRI segmentation', [], defDir);
    if isCancel
        return;
    end
    % Ask for number of vertices
    [nVertices, isCancel] = java_dialog('input', 'Number of vertices on the cortex surface:', 'FastSurfer segmentation', [], '15000');
    if isempty(nVertices) || isCancel
        return
    end
    nVertices = str2double(nVertices);
    % Open progress bar
    bst_progress('start', 'FastSurfer', 'FastSurfer MRI segmentation...');
    % Run FastSurfer
    isInteractive = 1;
    [isOk, errMsg] = Compute(iSubject, iAnatomy, nVertices, isInteractive, param, outdir);
    % Error handling
    if ~isOk
        bst_error(errMsg, 'FastSurfer MRI segmentation', 0);
    elseif ~isempty(errMsg)
        java_dialog('msgbox', ['Warning: ' errMsg]);
    end
    % Close progress bar
    bst_progress('stop');
end


text:
%-------------------------------------------------------------------------%
%   University of Illinois at Urbana Champaign
%   Department of Mechanical Science and Engineering
%   Laboratory of Photonics Research on Bio/Nano Environment
%   Written by : Tung Yuen Lau
%   Advisor: Prof. Kimani Toussaint
%   Start date: FEB 16 2012
%   Descriptiont: This is a simple function to detect edges in a 3D matrix
%   Format: image_pyramid(3D matrix of stacked images)
%-------------------------------------------------------------------------%
function final = canny3D(im, filsize, sigma, th_up, th_low)

im = double(im);

hfil = floor(filsize/2);
[w,h,d] = size(im);
[x y z] = meshgrid(-hfil:hfil,-hfil:hfil,-hfil:hfil);
fil_x = exp(-x.^2/(2*sigma^2))./(sigma*sqrt(2*pi)); clear x;
fil_y = exp(-y.^2/(2*sigma^2))./(sigma*sqrt(2*pi)); clear y;
fil_z = exp(-z.^2/(2*sigma^2))./(sigma*sqrt(2*pi)); clear z;
f = fil_x .* fil_y .* fil_z; clear fil_x; clear fil_y; clear fil_z;
f = f/sum(abs(f(:)));
imfil = imfilter(im,f,'replicate'); clear f;
[imfil_x , imfil_y, imfil_z] = gradient(imfil); clear imfil;

%%
%   Thinning (non-maximum suppression)
im_sub = nonmax_sup(imfil_x,imfil_y,imfil_z,th_up, th_low);
clear im_th; clear imfil_theta_z;

final = im_sub;

author: Sable

repo: mcbench-benchmarks


text:
function a2hbc_examples()
%%%%%%%%%%%%%%%%%%%%%%%%
% Examples using a2hbc %
%%%%%%%%%%%%%%%%%%%%%%%%
%
a2hbc_root_path = [fileparts(mfilename('fullpath')) filesep];

% Here you have several examples using a2hbc from the command line.

%% MIT format

% To warm up you can start using the automatic mode.

a2hbc( ...
    'recording_name', [ a2hbc_root_path '..' filesep 'example recordings' filesep '208.dat'], ...
    'recording_format', 'MIT', ...
    'op_mode', 'auto');

keyboard

% then you can try the semi-assisted mode, the algorithm only ask you for
% help when it doubts about a cluster label.

a2hbc( ...
    'recording_name', [ a2hbc_root_path '..' filesep 'example recordings' filesep '208.dat'], ...
    'recording_format', 'MIT', ...
    'op_mode', 'slightly-assisted');   % <----

keyboard

% or the assisted mode. In this case, the algorithms always ask you to
% label each cluster found.

a2hbc( ...
    'recording_name', [ a2hbc_root_path '..' filesep 'example recordings' filesep '208.dat'], ...
    'recording_format', 'MIT', ...
    'op_mode', 'assisted');           % <----

keyboard

% if you want you can control some parameters of the algorithm, as the
% number of clusters to look for.

a2hbc( ...
    'recording_name', [ a2hbc_root_path '..' filesep 'example recordings' filesep '208.dat'], ...
    'recording_format', 'MIT', ...
    'NumOfClusters', 3, ...          % <----
    'op_mode', 'assisted');

keyboard

% After the processing if you want to continue tunning up some parameters
% you can try this in order to use the GUI.

a2hbc( ...
    'recording_name', [ a2hbc_root_path '..' filesep 'example recordings' filesep '208.dat'], ...
    'recording_format', 'MIT', ...
    'NumOfClusters', 3, ...
    'InteractiveMode', true, ...      % <----
    'op_mode', 'assisted');

keyboard

%% Other formats, ISHNE

%% Other formats, AHA

author: marianux

repo: ecg-kit

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


text:
function [] = emailme(x)
%% Send email message when simulation finish to my Gmail account
% Based on:
% http://obasic.net/how-to-send-e-mail-from-matlab
%
% Coded by: Manuel Diaz

%% Set Parameters of my Gmail Account
setpref('Internet', 'E_mail', 'manuel.ade@gmail.com');
setpref('Internet', 'SMTP_Username', 'manuel.ade@gmail.com');
setpref('Internet', 'SMTP_Password', 'xxxxxx');
setpref('Internet', 'SMTP_Server', 'smtp.gmail.com');
props = java.lang.System.getProperties;
props.setProperty('mail.smtp.auth','true');
props.setProperty('mail.smtp.socketFactory.class', 'javax.net.ssl.SSLSocketFactory');
props.setProperty('mail.smtp.socketFactory.port', '465');

if x == 'yes'; % Send message
sendmail('manuel.ade@gmail.com', 'Mail from XXX', ...
    ['Hello Mr. Diaz! This is a mesage from XXX!' 10 ...
     'Your simulation has completed succesfully']);
end

repo: Aero-matlab


text:
function isOpen = hlp_pll_openPool(numWorkers,profileName)
% open a parallel toolbox pool

if nargin<2
    profileName = 'local';
end

sz = matlabpool('size');
isOpen = (sz>0);
% check if matlabpool is open...
if ~isOpen
    % ... if not, open one
    matlabpool(profileName,numWorkers);
elseif sz~=numWorkers
    % wrong number of workers
    % ...restart pool
    matlabpool('close');
    matlabpool(profileName,numWorkers);
end

repo: REST


text:
function density = calc_packing_density(dim,N,min_euclidean_dist)
%CALC_PACKING_DENSITY Density of packing given by minimum distance
%
%Syntax
% density = calc_packing_density(dim,N,min_euclidean_dist);
%
%Description
% DENSITY = CALC_PACKING_DENSITY(dim,N,MIN_EUCLIDEAN_DIST) sets DENSITY to
% be the density of a packing of S^dim by N equal spherical caps with centers
% having minimum Euclidean distance MIN_EUCLIDEAN_DIST.
%
% The argument dim must be a positive integer.
% The argument N must be a positive integer or an array of positive integers.
% The argument MIN_EUCLIDEAN_DIST must be an array of real nubers of the same array size as N.
% The result DENSITY will be an array of the same size as N.
%
%Notes
% The packing density is defined to be the sum of the areas of the spherical
% caps of the packing, divided by the area of the unit sphere S^dim.
%
% The spherical radius of the caps in the packing is half the minimum spherical
% distance between center points. The spherical radius for N == 1 is a special
% case. It is defined to be pi.
%
%Examples
% > N=2:6
%  N =
%       2     3     4     5     6
%
% > dist=eq_min_dist(2,N)
%  dist =
%      2.0000    1.4142    1.4142    1.4142    1.4142
%
% > calc_packing_density(2,N,dist)
%  ans =
%      1.0000    0.4393    0.5858    0.7322    0.8787
%
%See also
% EQ_MIN_DIST, AREA_OF_CAP, AREA_OF_SPHERE, EQ_PACKING_DENSITY

% Copyright 2004-2005 Paul Leopardi for the University of New South Wales.
% $Revision 1.10 $ $Date 2005-06-01 $
% Documentation files renamed
% $Revision 1.10 $ $Date 2005-05-27 $
% Function e2s changed name to euc2sph_dist
% $Revision 1.00 $ $Date 2005-02-12 $
%
% For licensing, see COPYING.
% For references, see AUTHORS.
% For revision history, see CHANGELOG.

%
% Check number of arguments
%
error(nargchk(3,3,nargin));
%
s_cap = euc2sph_dist(min_euclidean_dist)/2;
%
% The spherical radius for N == 1 is a special case. It is pi.
%
s_cap(N == 1) = pi;
density = N.*area_of_cap(dim, s_cap)./area_of_sphere(dim);
%
% end function

repo: gspbox


text:
%This Matlab script can be used to reproduce Figure 2.6 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Number of BS antennas
M = 100;

%Set the angle of the UE
theta = pi/6;

%Set the ASD
ASD = 10;

%Define the antenna spacing (in number of wavelengths)
antennaSpacing = 1/2; %Half wavelength distance

%Compute spatial correlation matrix with local scattering model and
%different angular distributions
R_Gaussian = functionRlocalscattering(M,theta,ASD,antennaSpacing,'Gaussian');
R_Uniform = functionRlocalscattering(M,theta,ASD,antennaSpacing,'Uniform');
R_Laplace = functionRlocalscattering(M,theta,ASD,antennaSpacing,'Laplace');

%Channel correlation matrix with uncorrelated fading
R_uncorrelated = eye(M);

%Extract the eigenvalues and place them in decreasing order
eigenvalues_Gaussian = flipud(eig(R_Gaussian));
eigenvalues_Uniform = flipud(eig(R_Uniform));
eigenvalues_Laplace = flipud(eig(R_Laplace));
eigenvalues_uncorr = flipud(eig(R_uncorrelated));

%Replace negative eigenvalues with a small positive number
%(since the correlation matrices should be Hermitian)
eigenvalues_Gaussian(eigenvalues_Gaussian<0) = 1e-16;
eigenvalues_Uniform(eigenvalues_Uniform<0) = 1e-16;
eigenvalues_Laplace(eigenvalues_Laplace<0) = 1e-16;

%% Plot the simulation results
figure;
hold on; box on;

plot(1:M,10*log10(eigenvalues_Laplace),'r--','LineWidth',1);
plot(1:M,10*log10(eigenvalues_Uniform),'b-.','LineWidth',1);
plot(1:M,10*log10(eigenvalues_Gaussian),'k','LineWidth',1);
plot(1:M,10*log10(eigenvalues_uncorr),'k:','LineWidth',1);

xlabel('Eigenvalue number in decreasing order');
ylabel('Normalized eigenvalue [dB]');
legend('Laplace','Uniform','Gaussian','Location','SouthEast');
ylim([-50 10]);

author: emilbjornson

repo: massivemimobook


text:
function tapas_hgf_binary_mab_plotTraj(r)
% Plots the estimated or generated trajectories for the binary HGF perceptual model for multi-armed
% bandit situations.
%
% Usage example:  est = tapas_fitModel(responses, inputs); tapas_hgf_binary_plotTraj(est);
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 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/>.

% Optional plotting of standard deviations (true or false)
plotsd = true;

% Optional plotting of responses (true or false)
ploty = true;

% Set up display
scrsz = get(0,'screenSize');
outerpos = [0.2*scrsz(3),0.2*scrsz(4),0.8*scrsz(3),0.8*scrsz(4)];
figure(...
    'OuterPosition', outerpos,...
    'Name', 'HGF trajectories');

% Set up colors
colors = [1 0 0; 0.67 0 1; 0 0.67 1; 0.67 1 0];

% Number of bandits
b = r.c_prc.n_bandits;

% Number of trials
n = size(r.u,1);

% Time axis
if size(r.u,2) > 1
    t = r.u(:,end)';
else
    t = ones(1,n);
end

ts = cumsum(t);
ts = [0, ts];

% Number of levels
try
    l = r.c_prc.n_levels;
catch
    l = (length(r.p_prc.p)+1)/5;
end

% Upper levels
for j = 1:l-2

    % Subplots
    subplot(l,1,j);

    if plotsd == true
        upperprior = r.p_prc.mu_0(l-j+1) +sqrt(r.p_prc.sa_0(l-j+1));
        lowerprior = r.p_prc.mu_0(l-j+1) -sqrt(r.p_prc.sa_0(l-j+1));
        upper = [upperprior; r.traj.mu(:,l-j+1,1)+sqrt(r.traj.sa(:,l-j+1,1))];
        lower = [lowerprior; r.traj.mu(:,l-j+1,1)-sqrt(r.traj.sa(:,l-j+1,1))];

        plot(0, upperprior, 'ob', 'LineWidth', 1);
        hold all;
        plot(0, lowerprior, 'ob', 'LineWidth', 1);
        fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
             'b', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    end
    plot(ts, [r.p_prc.mu_0(l-j+1); r.traj.mu(:,l-j+1,1)], 'b', 'LineWidth', 2);
    hold all;
    plot(0, r.p_prc.mu_0(l-j+1), 'ob', 'LineWidth', 2); % prior
    xlim([0 ts(end)]);
    title(['Posterior expectation of x_' num2str(l-j+1)], 'FontWeight', 'bold');
    ylabel(['\mu_', num2str(l-j+1)]);
end

% Second level
subplot(l,1,l-1)

if plotsd == true
    for j=1:b
        upperprior = r.p_prc.mu_0(2) +sqrt(r.p_prc.sa_0(2));
        lowerprior = r.p_prc.mu_0(2) -sqrt(r.p_prc.sa_0(2));
        upper = [upperprior; r.traj.mu(:,2,j)+sqrt(r.traj.sa(:,2,j))];
        lower = [lowerprior; r.traj.mu(:,2,j)-sqrt(r.traj.sa(:,2,j))];

        plot(0, upperprior, 'o', 'Color', colors(j,:), 'LineWidth', 1);
        hold all;
        plot(0, lowerprior, 'o', 'Color', colors(j,:), 'LineWidth', 1);
        fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
             colors(j,:), 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    end
end
for j=1:b
    plot(ts, [r.p_prc.mu_0(1); r.traj.mu(:,2,j)], 'Color', colors(j,:), 'LineWidth', 2);
    hold all;
    plot(0, r.p_prc.mu_0(1), 'o', 'Color', colors(j,:), 'LineWidth', 2); % prior
end
xlim([0 ts(end)]);
title('Posterior expectations of x_2', 'FontWeight', 'bold');
ylabel('\mu_2');

% Input level
subplot(l,1,l);

for j=1:b
    plot(ts, [tapas_sgm(r.p_prc.mu_0(2), 1); tapas_sgm(r.traj.mu(:,2,j), 1)], 'Color', colors(j,:), 'LineWidth', 2);
    hold all;
    plot(0, tapas_sgm(r.p_prc.mu_0(2), 1), 'o', 'Color', colors(j,:), 'LineWidth', 2); % prior
end
plot(ts(2:end), r.u(:,1), '.', 'Color', [0 0 0]); % inputs
plot(ts(2:end), r.traj.wt(:,1), 'k') % implied learning rate
if (ploty == true) && ~isempty(find(strcmp(fieldnames(r),'y'))) && ~isempty(r.y)
    y = r.y(:,1);
    if ~isempty(find(strcmp(fieldnames(r),'irr')))
        y(r.irr) = NaN; % weed out irregular responses
        plot(ts(r.irr)+1,  1.08.*ones([1 length(r.irr)]), 'x', 'Color', [1 0.7 0], 'Markersize', 11, 'LineWidth', 2); % irregular responses
    end
    for j=1:b
        plot(find(y==j), 1.08*ones([1 length(find(y==j))]), '.', 'Color', colors(j,:)); % responses
    end
    title(['Response y, input u (black dots), learning rate (fine black line), and posterior expectation of reward s(\mu_2) ', ...
           '(colour coded) for \rho=', num2str(r.p_prc.rho(2:end)), ', \kappa=', ...
           num2str(r.p_prc.ka(2:end)), ', \omega=', num2str(r.p_prc.om(2:end))], ...
      'FontWeight', 'bold');
    ylabel('y, u, s(\mu_2)');
    axis([0 ts(end) -0.15 1.15]);
else
    title(['Input u (black dots), learning rate (fine black line), and posterior expectation of input s(\mu_2) ', ...
           '(red) for \rho=', num2str(r.p_prc.rho(2:end)), ', \kappa=', ...
           num2str(r.p_prc.ka(2:end)), ', \omega=', num2str(r.p_prc.om(2:end))], ...
      'FontWeight', 'bold');
    ylabel('u, s(\mu_2)');
    axis([0 ts(end) -0.1 1.1]);
end
plot(ts(2:end), 0.5, 'k');
xlabel('Trial number');
hold off;

author: translationalneuromodeling

repo: tapas


text:
function fakeaxislabels(real,disp,whichax,hax,formatstr)

if ~exist('whichax','var'),
  whichax = 'x';
end
if ~exist('hax','var') || isempty(hax),
  hax = gca;
end
if ~exist('formatstr','var'),
  formatstr = '%f';
end
if strcmpi(whichax,'y')
  tick = 'ytick';
  ticklabel = 'yticklabel';
else
  tick = 'xtick';
  ticklabel = 'xticklabel';
end

labels = cell(1,length(disp));
for i = 1:length(disp),
  labels{i} = sprintf(formatstr,disp(i));
end

set(hax,tick,real,ticklabel,labels);

author: kristinbranson

repo: JAABA


text:
function ball_grid_display ( ng, xy )

%*****************************************************************************80
%
%% BALL_GRID_DISPLAY displays grid points inside a ball.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NG, the number of grid points inside the ball.
%
%    Input, real XY(3,NG), the grid points.
%
  scatter3 ( xy(1,:), xy(2,:), xy(3,:), 'b.' );
  axis equal
  title ( sprintf ( '%d grid points inside a ball', ng ) )
  grid on

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Y_pred, loglikelihoodX] = LDApredict(model,Gamma,X,classification,intercept)
% For an inferred LDA model and state timecourses, computes the
% likelihood manifold in X space and computes the likelihood function for
% each class of labels.
if nargin<4
    classification=true;
end
if nargin<5
    intercept=true;
end
[T2,K] = size(Gamma);
[T,nDimX] = size(X);
nDimY = size(model.state(1).W.Mu_W,2)-nDimX;
betas_mu = cell(K,1);
betas_sigma = cell(K,1);
for k = 1:K
    betas_mu{k} = model.state(k).W.Mu_W((nDimX+1):end,1:nDimX);
    if strcmp(model.train.covtype,'full') || strcmp(model.train.covtype,'uniquefull') ...
            || strcmp(model.train.covtype,'sharedfull')
        S = model.train.S==1;
        betas_sigma{k} = model.state(k).W.S_W(logical(S(:)),logical(S(:)));
    else
        betas_sigma{k} = model.state(k).W.S_W(1:nDimX,(nDimX+1):end,(nDimX+1):end);
    end
end

% Iterate through labels, noting that the distribution in x space can be
% computed as a sum for each condition:
if intercept
    numconds = size(betas_mu{1},1)-1;
else
    numconds = size(betas_mu{1},1);
end
betamu_givenY = zeros(numconds,K,nDimX);
betasig_givenY = zeros(numconds,K,nDimX,nDimX);
for testcond = 1:numconds
    % compute each state's distribution in x space - assuming intercept has
    % been used
    for k=1:K
        if intercept %account for mean term in generative model:
            %betamu_givenY(testcond,k,:) = betas_mu{k}(1,:) + betas_mu{k}(1 + testcond,:);
            betamu_givenY(testcond,k,:) = betas_mu{k}(1 + testcond,:);
        else
            betamu_givenY(testcond,k,:) = betas_mu{k}(testcond,:);
        end
        if strcmp(model.train.covtype,'full') || strcmp(model.train.covtype,'uniquefull') || ...
            strcmp(model.train.covtype,'sharedfull')
            betasig_givenY(testcond,k,:,:) = betas_sigma{k}([1:nDimY:nDimY*nDimX],[1:nDimY:nDimY*nDimX])+...
                betas_sigma{k}([testcond:nDimY:nDimY*nDimX],[testcond:nDimY:nDimY*nDimX]);
        else % Naive bayes classifier - just take uncertainty on each channel
            betasig_givenY(testcond,k,:,:) = diag(betas_sigma{k}(:,testcond + 1,testcond + 1)) + ...
                diag(betas_sigma{k}(:,1,1)) + 2*diag(betas_sigma{k}(:,1,testcond + 1));
        end
    end
end

% iterate through time, computing likelihood of X given Gamma:
if isfield(model,'Omega')
    if isfield(model.Omega,'pseudomean')
        CovMat = model.Omega.pseudomean(1:nDimX,1:nDimX);
    elseif strcmp(model.train.covtype,'uniquediag') || ...
           strcmp(model.train.covtype,'shareddiag')  % Naive Bayes case
        prec = diag(model.Omega.Gam_shape ./ model.Omega.Gam_rate(1:nDimX));
        CovMat = inv(prec);
    else % shared full covariance matrix
        prec = model.Omega.Gam_shape * model.Omega.Gam_irate(1:nDimX,1:nDimX);
        CovMat = inv(prec);
    end
else
    if strcmp(model.train.covtype,'diag')
        for k=1:K
            prec = diag(model.state(k).Omega.Gam_shape ./ model.state(k).Omega.Gam_rate(1:nDimX));
            CovMat(:,:,k) = inv(prec);
        end
    else % statewise full covariance
        for k=1:K
            prec = model.state(k).Omega.Gam_shape * model.state(k).Omega.Gam_irate(1:nDimX,1:nDimX);
            CovMat(:,:,k) = inv(prec);
        end
    end
end

% remove intercept from data if necessary:
if intercept
    for k=1:K
        X = X - Gamma(:,k) * betas_mu{k}(1,:);
    end
    Y_pred = zeros(T,nDimY-1);
else
    Y_pred = zeros(T,nDimY);
end
loglikelihoodX = zeros(size(Y_pred));
if T==T2
    for t=1:T
        CovMat_t = sum(CovMat .* repmat(permute(Gamma(t,:),[3,1,2]),nDimX,nDimX,1),3);
        if classification
            for testcond = 1:numconds
                mu = sum(squeeze(betamu_givenY(testcond,:,:)) .*repmat(Gamma(t,:)',1,nDimX));
                mu_rec(testcond,t,:) = mu;
                S = CovMat_t + squeeze(sum(betasig_givenY(testcond,:,:,:).*repmat(permute(Gamma(t,:)',[3,1,2]),1,1,nDimX,nDimX),2));
                loglikelihoodX(t,testcond) = -0.5*log(det(S)) -0.5 * (X(t,:) - mu) * inv(S) * (X(t,:) - mu)';
            end
            m = max(loglikelihoodX(t,:),[],2);
            Y_pred(t,:) = loglikelihoodX(t,:)==m;
            if sum(Y_pred(t,:))>1
                warning(['Equal scores achieved for multiple classes, t=',int2str(t),'\n\n']);
                a = find(Y_pred(t,:));
                Y_pred(t,a) = 0;
                Y_pred(t,a(randi(length(a)))) = 1; %randomly silence all but one of these entries
            end
        else
            betas_t = permute(sum(repmat(Gamma(t,:),[numconds,1,nDimX]).* betamu_givenY,2),[1,3,2]);

            %alternative implementation:
            if ~intercept
                sigma_Y = inv(inv(eye(nDimY)) + betas_t * inv(CovMat_t) * betas_t');
            else
                sigma_Y = inv(inv(eye(nDimY-1)) + betas_t * inv(CovMat_t) * betas_t');
            end
            betas_backwardmodel = inv(CovMat_t) * betas_t' * sigma_Y;
            Y_pred(t,:) = X(t,:)*betas_backwardmodel;
        end
    end
else
%implies Gamma invariant, so collapse over time:
    for testcond = 1:numconds
        mu = squeeze(betamu_givenY(testcond,:,:))';
        mu_rec(testcond,:) = mu;
        sig = CovMat + squeeze(betasig_givenY(testcond,:,:,:));
        loglikelihoodX(:,testcond) = -0.5 * log(det(sig)) -0.5 * sum(((X-repmat(mu,T,1))*inv(sig)).*(X-repmat(mu,T,1)),2);
    end
    m = max(loglikelihoodX(:,:),[],2);
    Y_pred = loglikelihoodX(:,:)==repmat(m,1,numconds);
end
if numconds==2 && all(Y_pred(:,1) == 1-Y_pred(:,2))
    % binary output:
    Y_pred = Y_pred(:,1);
end
end

author: OHBA-analysis

repo: HMM-MAR


text:
function [gx] = g_exp2d(x,phi,u,in)

gkaf = in.gkaf;
gkas = in.gkas;
%
% X = [gkaf,gkas,gkaf.^2,gkas.^2,gkaf.*gkas];
% X = [X,zeros(size(X,1),1)];
%
% gx = X*phi;
% dgdx = [];
% dgdphi = X';

gx = phi(1).*exp(0.5.*gkaf).*exp(-abs(gkas-0.5.*gkaf));
% gx = real(gx);

author: MBB-team

repo: VBA-toolbox


text:
function [p P_u P_pol] = omni_retro(u, pol)
%
% OMNI_RETRO  Retroproject pixel into 3D space.
%   P = OMNI_RETRO(UP, pol) retro-projects the point UP of the image
%   plane into the point P at depth S in 3D space
%
%   based on M=cam2world(m, ocam_model) from DAVIDE SCARAMUZZA
%   "Omnidirectional camera calibration toolbox" CAM2WORLD Project a give
%   pixel point onto the unit sphere M=CAM2WORLD=(m, ocam_model) returns
%   the 3D coordinates of the vector emanating from the single effective
%   viewpoint on the unit sphere
%
%
% Notes
% - pol: polynom coeff (vec)  -> "distortion" parametes in a0..an order, not polyval() order
%                                 we have 4th order polynome
%
%   Author: Davide Scaramuzza - email: davide.scaramuzza@ieee.org
%   Author: Grigory Abuladze - email: ryhor.a@google.com

%   Copyright (C) 2012 Grigory Abuladze @ ASL-vision.
%   Copyright (C) 2008 DAVIDE SCARAMUZZA, ETH Zurich

  if nargout == 1  % only points
    % Given an image point it returns the 3D coordinates of its correspondent optical ray
    r = (u(1,:).^2 + u(2,:).^2).^0.5;
    p = [u(1,:);
         u(2,:);
         polyval(pol(end:-1:1), r)];
    p = normc(p); % normalizes coordinates so that they have unit length (projection onto the unit sphere)

  else %  nargout > 2  ->  Jacobians

    if size(u,2) > 1
      error('Jacobians not available for multiple points');
    end

    r   = (u(1)^2 + u(2)^2)^0.5;
    n   = numel(pol);
    nn  = 0:n-1;     % orders vector
    r2n = r.^nn;     % powers vector
    p3  = r2n*pol;   %
    p   = [u(1);
           u(2);
           p3];

    [p,P_p] = normvec(p,1);

    a0 = pol(1); a1 = pol(2); a2 = pol(3); a3 = pol(4); a4 = pol(5);
    u1 = u(1);
    u2 = u(2);
    r2 = (u1^2 + u2^2);
    r  = r2^(1/2);

    % all before normalization
    P_u = [...
            [                                            1,                                            0]
            [                                            0,                                            1]
            [ 2*a2*u1 + 4*a4*u1*r2 + (a1*u1)/r + 3*a3*u1*r, 2*a2*u2 + 4*a4*u2*r2 + (a1*u2)/r + 3*a3*u2*r]
          ];

    P_pol = [...
              [ 0, 0,  0,        0,    0]
              [ 0, 0,  0,        0,    0]
              [ 1, r, r2, r2^(3/2), r2^2]
            ];

    % add normalization
    % f = norm(g()) -> use chain rule
    P_u = P_p * P_u;

  end

return

%% Jacobian
syms a0 a1 a2 a3 a4 real
syms u1 u2 s real

u   = [u1 u2]';
pol = [a0 a1 a2 a3 a4]';

p     = omni_retro(u, s, pol);

P_u   = simplify(jacobian(p,u))
P_s   = simplify(jacobian(p,s))
P_pol = simplify(jacobian(p,pol))

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function solplotl(sol,xy,x,y,fig)
%solplotl   plots nodal data on L-shaped domain
%   solplotl(sol,xy,x,y,fig);
%   input
%          sol        nodal solution vector
%          xy         nodal coordinate vector
%          x          vector of x-axis interpolation points
%          y          vector of y-axis interpolation points
%          fig        figure number
%
%   IFISS function: DJS; 30 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
fprintf('plotting solution... ')
% interpolate to a cartesian product mesh
[X,Y]=meshgrid(x,y);
xysol = griddata(xy(:,1),xy(:,2),sol,X,Y);
[II,JJ]=find(X<0 & Y<0); xysol(II,JJ)=nan;
figure(fig)
subplot(121),contour(X,Y,xysol,20),axis('square')
axis('off'),
if all([min(x),max(x),min(y),max(y)] == [-1,1,-1,1]), ellx, end
subplot(122),mesh(X,Y,xysol),axis('square')
view(330,30)
fprintf('done\n')
return


text:
%% FUNCTION Least_Dirty
%   Dirty Multi-Task Learning with Least Squares Loss.
%
%% INPUT
%   X: {n * d} * t - input matrix
%   Y: {n * 1} * t - output matrix
%   rho1: group sparsity regularization parameter
%   rho2: elementwise sparsity regularization parameter
%
%% OUTPUT
%   W: model: d * t
%   C: constant parameters
%   P: group sparsity structure (joint feature selection)
%   Q: elementwise sparsity component
%   funcVal: function (objective) value vector.
%   lossVal: loss value for every iteration.
%
%% LICENSE
%   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/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou, Pinghua Gong and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 3, 2012.
%
%% RELATED PAPERS
%
%   [1] Jalali, A. and Ravikumar, P. and Sanghavi, S. and Ruan, C. A dirty
%       model for multi-task learning, NIPS 2010.
%

function [W, C, P, Q, funcVal, lossVal] = Logistic_Dirty(X, Y, rho1, rho2, opts)

if nargin <4
    error('\n Inputs: X, Y, rho1, should be specified!\n');
end
X = multi_transpose(X);

if nargin <5
    opts = [];
end

task_num  = length (X);
dimension = size(X{1}, 1);
funcVal = [];
lossVal = [];

%initialize a starting point
C0_prep = zeros(1, task_num);
for t_idx = 1: task_num
    m1 = nnz(Y{t_idx} == 1);
    m2 = nnz(Y{t_idx} == -1);
    if ( m1==0 || m2==0 )
        C0_prep(t_idx) = 0;
    else
        C0_prep(t_idx) = log(m1/m2);
    end
end

if opts.init==2
    P0 = zeros(dimension, task_num);
    Q0 = zeros(dimension, task_num);
    C0 = zeros(1, task_num);
elseif opts.init== 0
    P0 = randn(dimension, task_num);
    Q0 = randn(dimension, task_num);
    C0 = C0_prep;
elseif opts.init== 1
    if isfield(opts,'P0')
        P0=opts.P0;
        if (nnz(size(P0)-[dimension, task_num]))
            error('\n Check the input .P0');
        end
    else
        error('\n check opt.init');
    end

     if isfield(opts,'Q0')
        Q0=opts.Q0;
        if (nnz(size(Q0)-[dimension, task_num]))
            error('\n Check the input .Q0');
        end
    else
        error('\n check opt.init');
    end

    if isfield(opts,'C0')
        C0=opts.C0;
    else
        error('\n check opt.init');
    end
end

Pz= P0;
Qz= Q0;
Cz= C0;
Pz_old = P0;
Qz_old = Q0;
Cz_old = C0;

t = 1;
t_old = 0;
iter = 0;
gamma = 1;
gamma_inc = 2;

while iter < opts.maxIter
    alpha = (t_old - 1) /t;

    Ps = (1 + alpha) * Pz - alpha * Pz_old;
    Qs = (1 + alpha) * Qz - alpha * Qz_old;
    Cs = (1 + alpha) * Cz - alpha * Cz_old;

    % compute function value and gradients of the search point
    [gWs, gCs, Fs ]  = gradVal_eval(Ps+Qs, Cs);

    % the Armijo Goldstein line search scheme
    while true
        Pzp = proximalL1infnorm(Ps - gWs/gamma, rho1 / gamma);
        Qzp = proximalL11norm(Qs - gWs/gamma, rho2 / gamma);

        Czp = Cs - gCs/gamma;
        Fzp = funVal_eval(Pzp+Qzp, Czp);

        delta_Pzp = Pzp - Ps;
        delta_Qzp = Qzp - Qs;
        delta_Czp = Czp - Cs;
        nm_delta_Pzp=norm(delta_Pzp, 'fro')^2;
        nm_delta_Qzp=norm(delta_Qzp, 'fro')^2;
        nm_delta_Czp=norm(delta_Czp, 'fro')^2;

        r_sum = (nm_delta_Pzp+nm_delta_Qzp+nm_delta_Czp)/2;

%         Fzp_gamma = Fs + sum(sum((delta_Pzp+delta_Qzp).* gWs))...
%             + sum(sum(delta_Czp .* gCs))...
%             + gamma/2 * norm(delta_Pzp + delta_Qzp, 'fro')^2 ...
%             + gamma/2 * nm_delta_Czp;

        Fzp_gamma = Fs + sum(sum((delta_Pzp+delta_Qzp).* gWs))...
            + sum(sum(delta_Czp .* gCs))...
            + gamma/2 * (nm_delta_Pzp +nm_delta_Qzp) ...
            + gamma/2 * nm_delta_Czp;

        if (r_sum <=1e-20)
            break;
        end

        if (Fzp <= Fzp_gamma)
            break;
        else
            gamma = gamma * gamma_inc;
        end
    end

    Pz_old = Pz;
    Qz_old = Qz;
    Cz_old = Cz;
    Pz = Pzp;
    Qz = Qzp;
    Cz = Czp;

    funcVal = cat(1, funcVal, Fzp + rho1*L1infnorm(Pzp) + rho2*L11norm(Qzp));
    lossVal = cat(1, lossVal, Fzp);

    % test stop condition.
    switch(opts.tFlag)
        case 0
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <= opts.tol)
                    break;
                end
            end
        case 1
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <=...
                        opts.tol* funcVal(end-1))
                    break;
                end
            end
        case 2
            if ( funcVal(end)<= opts.tol)
                break;
            end
        case 3
            if iter>=opts.maxIter
                break;
            end
    end

    iter = iter + 1;
    t_old = t;
    t = 0.5 * (1 + (1+ 4 * t^2)^0.5);

end

P=Pzp;
Q=Qzp;
W = Pzp+Qzp;
C = Czp;

% private functions

    function [grad_W, grad_C, funcVal] = gradVal_eval(W, C)
        grad_W = zeros(dimension, task_num);
        grad_C = zeros(1, task_num);
        lossValVect = zeros (1 , task_num);
        for i = 1:task_num
                [ grad_W(:, i), grad_C(:, i), lossValVect(:, i)] = unit_grad_eval( W(:, i), C(i), X{i}, Y{i});
        end
        % here when computing function value we do not include
        % l1 norm.
        funcVal = sum(lossValVect);
    end

    function [funcVal] = funVal_eval (W, C)
        funcVal = 0;
        for i = 1: task_num
            funcVal = funcVal + unit_funcVal_eval( W(:, i), C(i), X{i}, Y{i});
        end
        % here when computing function value we do not include
        % l1 norm.
    end

end

function [ grad_w, grad_c, funcVal ] = unit_grad_eval( w, c, x, y)
    %gradient and logistic evaluation for each task
    m = length(y);
    weight = ones(m, 1)/m;
    weighty = weight.* y;
    aa = -y.*(x'*w + c);
    bb = max( aa, 0);
    funcVal = weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb );
    pp = 1./ (1+exp(aa));
    b = -weighty.*(1-pp);
    grad_c = sum(b);
    grad_w = x * b;
end

function [ funcVal ] = unit_funcVal_eval( w, c, x, y)
    %function value evaluation for each task
    m = length(y);
    weight = ones(m, 1)/m;
    aa = -y.*(x'*w + c);
    bb = max( aa, 0);
    funcVal = weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb );
end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


function [region_meta_info, region_meta_info_true]=preprocess_cpmc_candidates(imnames, cpmcdir, cpmcgtdir, gtdir, ovoutdir, N)
if(~exist(ovoutdir, 'file'))
    mkdir(ovoutdir);
end

for i=1:numel(imnames)
    %load sprep
    x1=load(fullfile(cpmcdir, [imnames{i} '.mat']));
    sp=x1.sp;
    reg2sp=x1.sp_app;

    %compute overlaps between everything
    ov=double(reg2sp')*double(reg2sp);
    ov=ov>0;
    if(~exist(fullfile(ovoutdir, [imnames{i} '.mat']), 'file'))
    save(fullfile(ovoutdir, [imnames{i} '.mat']), 'ov');
    end
    %load gt
    if(isempty(gtdir))
        %no gt
        inst=zeros(size(sp));
        categories=[];
    else
    [cls, inst, categories]=load_gt(gtdir, imnames{i});
    end
    x2=load(fullfile(cpmcgtdir, imnames{i}));
    %compute overlap of projected ground truth with ground truth
    gt_overlap=get_gt_overlaps(x2.sp_app, x2.sp, inst);

    %for each projected gt get the instance it corresponds to
    [m, assigned_index]=max(gt_overlap,[],1);

    %compute overlaps with ground truth
    overlap=get_gt_overlaps(reg2sp, sp, inst);
    box_overlap=get_gt_overlaps_box(reg2sp, sp, inst);
    %populate meta info
    region_meta_info.num_regs(i)=size(reg2sp,2)+numel(categories(assigned_index));
    region_meta_info.overlaps{i}=overlap(assigned_index,:);
    region_meta_info.box_overlaps{i}=box_overlap(assigned_index,:);
    region_meta_info.gt{i}=[categories(assigned_index(:))' zeros(1,size(reg2sp,2))];

    region_meta_info_true.num_regs(i)=size(reg2sp,2)+numel(categories);
    region_meta_info_true.overlaps{i}=overlap;
    region_meta_info_true.box_overlaps{i}=box_overlap;
    region_meta_info_true.gt{i}=[categories(:)' zeros(1,size(reg2sp,2))];

    fprintf('Done %d\n', i);
end

function write_sprep_text(sp, filename)
fid=fopen(filename, 'w');
fprintf(fid,'%d %d\n', size(sp,1), size(sp,2));
fprintf(fid, '%d ', sp(:));
fclose(fid);

function [sp, sp2reg]=intersect_partitions(sp1, sp2, sp2reg1, sp2reg2)
%do a unique
[unique_pairs, junk, sp]=unique([sp1(:) sp2(:)], 'rows');
sp=reshape(sp, size(sp1));

%create new sp2reg
sp2reg=[sp2reg1(unique_pairs(:,1),:) sp2reg2(unique_pairs(:,2),:)];

%break disconnected sp
 spNew = zeros(size(sp));
  sp2regNew = false(size(sp2reg));
  cnt = 0;
  %% Check connectivity of superpixels
  for i = 1:max(sp(:))
    tt = bwconncomp(sp == i);
    for j = 1:tt.NumObjects,
      cnt = cnt+1;
      spNew(tt.PixelIdxList{j}) = cnt;
      sp2regNew(cnt, :) = sp2reg(i, :);
    end
  end

sp=spNew;
sp2reg=sp2regNew;

function [spNew sp2regNew] = compressSp2reg(sp, sp2reg)
  assert(islogical(sp2reg), 'sp2reg should be logical');
  [sp2reg bb spSame] = unique(sp2reg, 'rows');
  sp = spSame(sp);

  spNew = zeros(size(sp));
  sp2regNew = false(size(sp2reg));
  cnt = 0;
  %% Check connectivity of superpixels
  for i = 1:max(sp(:))
    tt = bwconncomp(sp == i);
    for j = 1:tt.NumObjects,
      cnt = cnt+1;
      spNew(tt.PixelIdxList{j}) = cnt;
      sp2regNew(cnt, :) = sp2reg(i, :);
    end
  end


text:
classdef TestExposureCompensator
    %TestExposureCompensator

    methods (Static)
        function test_1
            etypes = {'GainCompensator', 'BlocksGainCompensator'};
            for i=1:numel(etypes)
                obj = cv.ExposureCompensator(etypes{i});
                typename = obj.typeid();
                validateattributes(typename, {'char'}, {'row', 'nonempty'});

                %{
                %TODO
                obj.feed(corners, images, masks);
                img = obj.apply(index, corner, img, mask);
                if strcmp(etypes{i}, 'GainCompensator')
                    g = obj.gains();
                end
                %}
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% SetValueOfAssignment Sets the value of a variable assignment in a factor.
%
%   F = SetValueOfAssignment(F, A, v) sets the value of a variable assignment,
%   A, in factor F to v. The order of the variables in A are assumed to be the
%   same as the order in F.var.
%
%   F = SetValueOfAssignment(F, A, v, VO) sets the value of a variable
%   assignment, A, in factor F to v. The order of the variables in A are given
%   by the vector VO.
%
%   Note that SetValueOfAssignment *does not modify* the factor F that is
%   passed into the function, but instead returns a modified factor with the
%   new value(s) for the specified assignment(s). This is why we have to
%   reassign F to the result of SetValueOfAssignment in the code snippets
%   shown above.
%
%   See also GetValueOfAssignment.m and SampleFactors.m
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function F = SetValueOfAssignment(F, A, v, VO)

if (nargin == 3),
    indx = AssignmentToIndex(A, F.card);
else
    map = zeros(length(F.var), 1);
    for i = 1:length(F.var),
        map(i) = find(VO == F.var(i));
    end;
    indx = AssignmentToIndex(A(map), F.card);
end;

F.val(indx) = v;

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
classdef(Abstract) AbstractGradientOptimizer  < AbstractOptimizer

    methods
        setGradientCalculationMethod(obj, newGradCalcMethod);
    end

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function z=getUVDirection(zC,zRx,M,includeW)
%%GETUVDIRECTION Convert a position into local direction cosines [u;v] or
%         [u;v;w]. Direction cosines u and v are just the x and y
%         coordinates of a unit vector from the receiver to the target in
%         the coordinate system at the receiver. This basically assumes
%         that the boresight direction of the receiver is the z axis.
%         Assuming the target is in front of the receiver, the third unit
%         vector coordinate is not needed. However, with the includeW
%         option, it can be provided, resulting in r-u-v-w coordinates.
%
%INPUT: zC A 3XN matrix of points in global [x;y;z] Cartesian coordinates.
%      zRx The 3XN [x;y;z] location vectors of the receivers in Cartesian
%          coordinates.  If this parameter is omitted or an empty matrix is
%          passed, then the receivers are assumed to be at the origin. If
%          only a single vector is passed, then the receiver location is
%          assumed the same for all of the target states being converted.
%        M A 3X3XN hypermatrix of the rotation matrices to go from the
%          alignment of the global coordinate system to that at the
%          receiver. The z-axis of the local coordinate system of the
%          receiver is the pointing direction of the receiver. If omitted
%          or an empty matrix is passed, then it is assumed that the local
%          coordinate system is aligned with the global and M=eye(3) --the
%          identity matrix is used. If only a single 3X3 matrix is passed,
%          then is=t is assumed to be the same for all of the N
%          conversions. Typically, if includeW is false, it doesn't make
%          sense to return the u-v values that are not local, so one will
%          usually omit M in that instance.
%  includeW An optional boolean value indicating whether a third direction
%          cosine component should be included. The u and v direction
%          cosines are two parts of a 3D unit vector. Generally, one might
%          assume that the target is in front of the sensor, so the third
%          component would be positive and is not needed. However, the
%          third component can be included if ambiguity exists. The default
%          if this parameter is omitted or an empty matrix is passed is
%          false.
%
%OUTPUT: z The 2XN (or 3XN if includeW is true) matrix of direction cosines
%          of the converted points in the form [u;v] or [u;v;w].
%
%Details of the conversion are given in [1].
%
%REFERENCES:
%[1] David F. Crouse , "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%February 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

N=size(zC,2);

if(nargin<4||isempty(includeW))
    includeW=false;
end

if(nargin<3||isempty(M))
    M=repmat(eye(3),[1,1,N]);
elseif(size(M,3)==1)
    M=repmat(M,[1,1,N]);
end

if(nargin<2||isempty(zRx))
    zRx=zeros(3,N);
elseif(size(zRx,2)==1)
    zRx=repmat(zRx,[1,N]);
end

%Allocate space for the return values.
if(includeW==true)
    z=zeros(3,N);
else
    z=zeros(2,N);
end
for curPoint=1:N
    %The target location in the receiver's coordinate system.
    zCL=M(:,:,curPoint)*(zC(:,curPoint)-zRx(1:3,curPoint));

    %Perform the conversion.
    r1=norm(zCL);%Receiver to target.

    u=zCL(1)/r1;
    v=zCL(2)/r1;

    z(1:2,curPoint)=[u;v];
    if(includeW)
        z(3,curPoint)=zCL(3)/r1;
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function shortcutUtils = GetShortcutUtils()
%GETSHORTCUTUTILS Gets an instance of ShortcutUtils.
%
%  SHORTCUTUTILS = GETSHORTCUTUTILS() gets an instance of ShortcutUtils.
%
% Examples:
%
% shortcutUtils = GetShortcutUtils();
% methods(shortcutUtils)
% methodsview(shortcutUtils)

% $Author: rcotton $	$Date: 2010/08/23 14:35:09 $	$Revision: 1.1 $
% Copyright: Health and Safety Laboratory 2010

CheckForJVM();
shortcutUtils = com.mathworks.mlwidgets.shortcuts.ShortcutUtils;

end

repo: mcbench-benchmarks


text:
function obj=prepare(obj)
% prepare for collecting coverage information
%
% obj=prepare(obj)
%
% Input:
%   obj                 MOcovMFileCollection instance
%
% Output:
%   obj                 MOcovMFileCollection instance with, if
%                       method=='file', the original path restored as well
%                       as mfiles rewritten to collect coverage
%                       information.

    monitor=obj.monitor;

    fns=mocov_find_files(obj.root_dir,'*.m',monitor,obj.exclude_pat);
    n=numel(fns);

    mfiles=cell(n,1);
    for k=1:n
        fn=fns{k};
        mfiles{k}=MOcovMFile(fn);
    end

    obj.mfiles=mfiles;

    if ~ischar(obj.method)
        error('method must be char, found %s', class(obj.method));
    end

    switch obj.method
        case 'profile'
            profile on
            notify(monitor,'','Enabling profiler');

        case 'file'
            % store original path
            obj.orig_path=path();
            notify(monitor,sprintf('Preserving original path\n'));

            temp_dir=tempname();
            notify(monitor,sprintf('Rewriting m-files\n'));
            obj=rewrite_mfiles(obj,temp_dir);

            addpath(genpath(temp_dir));
            notify(monitor,'',sprintf('Path is: %s\n', path()));

        otherwise
            error('illegal method %s', obj.method);
    end

    notify(monitor,sprintf('Coverage preparation complete\n'));

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%%  RDT testing for VISTASOFT
%
% Some examples and testing of the download of vistadata from the
% RemoteDataToolbox (RDT) site.

%% Create an RDT client object

rd = RdtClient('vistasoft');

% Here is a look at the general repository
rd.openBrowser('fancy',true)

%% Download examples

% First some examples from the AFQ remote path.
rd.crp('/AFQ/templates');

% This is a list of the artifacts and we print it
a = rd.listArtifacts('print',true);

% Notice they have a variety of types.  We automatically handle matlab file
% types by loading them.  Other types may need special handling.

%% For Matlab files, the handling is simpler

% Here are just the matlab types
a = rd.listArtifacts('print',true,'type','mat');

% In this case, the returned variable (data) is a struct with each of the
% variables in the Remote Matlab file
data = rd.readArtifact(a(1))

%% Downloading a text file

% Here is a text file, the outcome is putting the text file in a default
% location with a lousy, default name.
a = rd.listArtifacts('print',true,'type','txt');

fname = rd.readArtifact(a(1),'type','txt');

% This is the file that was downloaded by default.  It is ugly and includes
% terms thare are used by the database.  Readable, but ugly.  The advantage
% of doing it this way is the server will know that you downloaded it and
% if you request it again, it will just use the cached, local copy.
disp(fname)

% An alternative is that you can control the file destination. Used in this way,
% the file name is better and you control the destination folder precisely.
dFolder = fullfile(vistaRootPath,'local');
fname = rd.readArtifact(a(1),'destinationFolder',dFolder);
dir(dFolder)

% Notice that the file name is correct here, without all the annoying
% database formatting.  It is yours.  The server doesn't know about it.
disp(fname)

%% nii.gz files are stored as type 'gz'.  Sorry about that

% Here they are.  Vistasoft has a lot of nii.gz files in general
a = rd.listArtifacts('print',true,'type','gz');

% Here is the download with the destination specified
fname = rd.readArtifact(a(1),'destinationFolder',dFolder);

% But the file that is written out is good
disp(fname)

% And niftiRead works on the downloaded fname.
ni = niftiRead(fname)

%% Recursive listing

rd.crp('/vistadata/anatomy');

a = rd.listArtifacts('print',true,'recursive',true);

%% To upload vista data files requires a credential.

% See the files rdtPublishFunctional and rdtPublishAFQ for the methods used
% to move AFQ data and vistadata up to the RemoteDataToolbox site.

% cd(vistadata)
rd = RdtClient('vistasoft');
rd.credentialsDialog;

%%

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ output_struct ] = parsesubapertureinputs( varargin )
%PARSESUBAPERTUREINPUTS Parse input arguments for subaperture processing
%
% Used in both in-memory and in-file version of subaperture processing
%
% Written by: Wade Schwartzkopf
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

p = inputParser;
p.KeepUnmatched=true;
p.addParamValue('frames', 7, @(x) isscalar(x)&&(x>=2)&&(mod(x,1)==0));
p.addParamValue('apfraction', 0.25, @(x) isscalar(x)&&(x<1)&&(x>0));
p.addParamValue('method', 'normal', @(x) any(strcmpi(x,{'minimal', 'normal','fullpixel'})));
p.addParamValue('platformdir', 'right', @(x)any(strcmpi(x,{'right','left'})));
p.addParamValue('dim', 1, @(x) isequal(x,1)||isequal(x,2));
p.addParamValue('fill', 1, @(x) isscalar(x)&&(x>=1));
p.FunctionName = 'SUBAPERTURE';
p.parse(varargin{:});

output_struct=p.Results;
output_struct.reverse_frames=(output_struct.dim==1)&&strcmp(output_struct.platformdir, 'right');

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function A=linspace(x1,x2,N)
%LINPSPACE (overloaded)

if nargin < 3
    N = 100;
elseif isa(N,'sdpvar')
    error('SDPVAR/LINSPACE does not support N as a decision variable');
end

A = x1 + linspace(0,1,N)*(x2-x1);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [StopFlag, Status] = StopCritBDCA(nfxk, Niter, Nmap, T, MaxNumIter, MaxNumMapEval, TimeLimit, epsilon, Stopping_Crit)
% Function checking that one of the stopping criteria
% holds to terminate LLM and GLM. It perepare the status determining why
% the algorithm is stopped.
%
% USAGE:
%
%    [StopFlag, Status] = StopCritBDCA(nfxk, Niter, Nmap, T, MaxNumIter, MaxNumMapEval, TimeLimit, epsilon, Stopping_Crit)
%
% INPUTS:
%    nhxk:             the norm 2 of `h(xk)`
%    Niter:            the number of iterations
%    Nmap:             the number of mapping calls
%    T:                the running time
%    MaxNumIter:       maximum number of iterations
%    MaxNumMapEval:    maximum number of function evaluations
%    TimeLimit:        maximum running time
%    epsilon:          accuracy parameter
%    Stopping_Crit:    stopping criterion:
%
%                        1. stop if :math:`||nfxk|| \leq \epsilon`
%                        2. stop if `MaxNumIter` is reached
%                        3. stop if `MaxNumMapEval` is reached
%                        4. stop if `TimeLimit` is reached
%                        5. stop if (default) :math:`||hxk|| \leq \epsilon` or `MaxNumIter` is reached
%
% OUTPUTS:
%    StopFlag:         1: if one of the stopping criteria holds, 0: if none of the stopping criteria holds
%    Status:           the reason of the scheme termination

switch Stopping_Crit

  case 1
    if nhxk <= epsilon
      StopFlag = 1;
      Status   = 'A solution of nonlinear system is found.';
    else
      StopFlag = 0;
      Status   = [];
    end

  case 2
    if Niter >= MaxNumIter
      StopFlag = 1;
      Status   = 'Maximum number of iterations is reached.';
    else
      StopFlag = 0;
      Status   = [];
    end

  case 3
    if Nmap >= MaxNumMapEval
      StopFlag = 1;
      Status   = 'Maximum number of mapping evaluations is reached.';
    else
      StopFlag = 0;
      Status   = [];
    end

  case 4
    if T >= TimeLimit
      StopFlag = 1;
      Status   = 'Time limit is reached.';
    else
      StopFlag = 0;
      Status   = [];
    end

  case 5
    if (nfxk <= epsilon || Niter >= MaxNumIter)
      StopFlag = 1;
      if Niter < MaxNumIter
          Status = 'a solution is found';
      else
          Status = 'Maximum number of iterations is reached.';
      end
    else
      StopFlag = 0;
      Status   = [];
    end

end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% End of StopCritBDCA.m %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: opencobra

repo: cobratoolbox


text:
function fused = Ying_2017_CAIP(I, mu, k, a, b) % camera a, b
%%
% @inproceedings{ying2017new,
%   title={A New Image Contrast Enhancement Algorithm Using Exposure Fusion Framework},
%   author={Ying, Zhenqiang and Li, Ge and Ren, Yurui and Wang, Ronggang and Wang, Wenmin},
%   booktitle={International Conference on Computer Analysis of Images and Patterns},
%   pages={36--46},
%   year={2017},
%   organization={Springer}
% }
%
% Please feel free to contact me (yingzhenqiang-at-gmail-dot-com) if you
% have any questions or concerns.

if  ~exist( 'mu', 'var' )
    mu = 0.5;
end

if ~exist( 'a', 'var' )
    a = -0.3293;
end

if ~exist( 'b', 'var' )
    b = 1.1258;
end

if ~isfloat(I)
    I = im2double( I );
end

lambda = 0.5;
sigma = 5;

%% t: scene illumination map
t_b = max( I, [], 3 ); % also work for single-channel image
t_our =  imresize( tsmooth( imresize( t_b, 0.5 ), lambda, sigma ), size( t_b ) );

%% k: exposure ratio
if  ~exist( 'k', 'var' ) || isempty(k)
    isBad = t_our < 0.5;
    J = maxEntropyEnhance(I, isBad);
else
    J = applyK(I, k, a, b); %k
    J = min(J, 1); % fix overflow
end

%% W: Weight Matrix
t = repmat(t_our, [1 1 size(I,3)]);
W = t.^mu;

I2 = I.*W;
J2 = J.*(1-W);

fused = I2 + J2;

    function J = maxEntropyEnhance(I, isBad)
        Y = rgb2gm(real(max(imresize(I, [50 50]), 0))); % max - avoid complex number

        if exist('isBad', 'var')
            isBad = (imresize(isBad, [50 50]));
            Y = Y(isBad);
        end

        if isempty(Y)
           J = I; % no enhancement k = 1
           return;
        end

        opt_k = fminbnd(@(k) ( -entropy(applyK(Y, k)) ),1, 7);
        J = applyK(I, opt_k, a, b) - 0.01;

    end
end

function I = rgb2gm(I)
if size(I,3) == 3
    I = im2double(max(0,I)); % negative double --> complex double
    I = ( I(:,:,1).*I(:,:,2).*I(:,:,3) ).^(1/3);
end
end

function J = applyK(I, k, a, b)

if ~exist( 'a', 'var' )
    a = -0.3293;
end

if ~exist( 'b', 'var' )
    b = 1.1258;
end

f = @(x)exp((1-x.^a)*b);
beta = f(k);
gamma = k.^a;
J = I.^gamma.*beta;
end

function S = tsmooth( I, lambda, sigma, sharpness)
if ( ~exist( 'lambda', 'var' ) )
    lambda = 0.01;
end
if ( ~exist( 'sigma', 'var' ) )
    sigma = 3.0;
end
if ( ~exist( 'sharpness', 'var' ) )
    sharpness = 0.001;
end
I = im2double( I );
x = I;
[ wx, wy ] = computeTextureWeights( x, sigma, sharpness);
S = solveLinearEquation( I, wx, wy, lambda );
end

function [ W_h, W_v ] = computeTextureWeights( fin, sigma, sharpness)

dt0_v = [diff(fin,1,1);fin(1,:)-fin(end,:)];
dt0_h = [diff(fin,1,2)';fin(:,1)'-fin(:,end)']';

gauker_h = filter2(ones(1,sigma),dt0_h);
gauker_v = filter2(ones(sigma,1),dt0_v);
W_h = 1./(abs(gauker_h).*abs(dt0_h)+sharpness);
W_v = 1./(abs(gauker_v).*abs(dt0_v)+sharpness);

end

function OUT = solveLinearEquation( IN, wx, wy, lambda )
[ r, c, ch ] = size( IN );
k = r * c;
dx =  -lambda * wx( : );
dy =  -lambda * wy( : );
tempx = [wx(:,end),wx(:,1:end-1)];
tempy = [wy(end,:);wy(1:end-1,:)];
dxa = -lambda *tempx(:);
dya = -lambda *tempy(:);
tempx = [wx(:,end),zeros(r,c-1)];
tempy = [wy(end,:);zeros(r-1,c)];
dxd1 = -lambda * tempx(:);
dyd1 = -lambda * tempy(:);
wx(:,end) = 0;
wy(end,:) = 0;
dxd2 = -lambda * wx(:);
dyd2 = -lambda * wy(:);

Ax = spdiags( [dxd1,dxd2], [-k+r,-r], k, k );
Ay = spdiags( [dyd1,dyd2], [-r+1,-1], k, k );

D = 1 - ( dx + dy + dxa + dya);
A = (Ax+Ay) + (Ax+Ay)' + spdiags( D, 0, k, k );

if exist( 'ichol', 'builtin' )
    L = ichol( A, struct( 'michol', 'on' ) );
    OUT = IN;
    for ii = 1:ch
        tin = IN( :, :, ii );
        [ tout, ~ ] = pcg( A, tin( : ), 0.1, 50, L, L' );
        OUT( :, :, ii ) = reshape( tout, r, c );
    end
else
    OUT = IN;
    for ii = 1:ch
        tin = IN( :, :, ii );
        tout = A\tin( : );
        OUT( :, :, ii ) = reshape( tout, r, c );
    end
end
end

author: baidut

repo: OpenCE

save_path: github-repos/MATLAB/baidut-OpenCE


text:
function [models,M] = esvm_subset_of_models(models,M,subset)
%Choose a subset of models and a subset of the calibration matrix
% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

if ~exist('subset','var')
  s1 = cellfun(@(x)size(x.model.w,1),models);
  s2 = cellfun(@(x)size(x.model.w,2),models);
  urows = unique([s1; s2]','rows');
  [c,d] = ismember([s1; s2]',urows,'rows');
  [aa,bb] = max(hist(d,1:size(urows,1)));
  target = urows(bb,:);
  subset = find(s1==target(1) & s2==target(2));
end
hits = repmat(1:length(models),2,1);

models = models(subset);

subset2 = find(ismember(hits(:),subset));
M.w = cellfun2(@(x)x(subset2),M.w);
M.w = M.w(subset2);
M.b = M.b(subset2);
M.C = M.C(subset2,subset2);

author: quantombone

repo: exemplarsvm


text:
function AC = simpleTRsizing(AC)
% gives us R, sigma, Omega, c, TS, and b for tail rotor

%% check TSoffset value
TSoffset = 7.5; %m/s slower than main rotor
ARgoal = 6;

AC.TRotor.TS = AC.Rotor.TS - TSoffset;

discarea = AC.Rotor.DA;
discloading = AC.W.AUM./discarea;
ratio = 1./(7.15-.0553*discloading);
AC.TRotor.R = AC.Rotor.R.*ratio;

AC.TRotor.DA = pi*AC.TRotor.R.^2;

sigma_mr = AC.Rotor.sigma;
AC.TRotor.sigma = 1.5277 * sigma_mr +.0461; % Trendline from historical aircraft

AC.TRotor.b = round(pi*AC.TRotor.sigma*ARgoal);

AC.TRotor.c = AC.TRotor.sigma*AC.TRotor.DA./(AC.TRotor.b*AC.TRotor.R);

AC.TRotor.Omega = AC.TRotor.TS/AC.TRotor.R;

AC.TRotor.X = AC.Rotor.R+1.1*AC.TRotor.R;

AC.TRotor.BA = AC.TRotor.sigma*AC.TRotor.DA;
end

author: Sable


text:
% LIONSIMBA example script
% Different heat exchange coefficients: this script provides the example simulation shown in the
% paper.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

% Clear the workspace
clear

% Define the integration times.
t0 = 0;
tf = 10^4;
% Define the parameters structure. By default the hcell parameter is set to
% 1 [W / (m^2 K)]
param{1} = Parameters_init;

% Change the hcell parameter and set it to 0.01 [W / (m^2 K)]
param{1}.hcell = 0.01;

% Start the simulation. Note that the final integration time is 10^4 and
% LIONSIMBA will stop automatically when reached the Cutoff Voltage of
% 2.5V. Note that no Jacobian among one simulation and the other. This is
% due to the fact that the model structure changes for each of the
% simulation (because hcell changes) and consequently also the Jacobian
% matrix changes and cannot be used for all the scenarios.

out1 = startSimulation(t0,tf,[],-30,param);

% Change the hcell parameter and set it to 1 [W / (m^2 K)]
param{1}.hcell = 1;

% Run the simulation
out2 = startSimulation(t0,tf,[],-30,param);

% Change the hcell parameter and set it to 100 [W / (m^2 K)]
param{1}.hcell = 100;

% Run the simulation
out3 = startSimulation(t0,tf,[],-30,param);

%% Plot the results

figure(1)
plot(out1.time{1},out1.Temperature{1}(:,end),'LineWidth',6)
hold on
plot(out2.time{1},out2.Temperature{1}(:,end),'--','LineWidth',6)
plot(out3.time{1},out3.Temperature{1}(:,end),'-.','LineWidth',6)
xlabel('Time [s]')
ylabel('Temperature [K]')
grid on
box on

figure(2)
plot(out1.time{1},out1.Voltage{1},'LineWidth',6)
hold on
plot(out2.time{1},out2.Voltage{1},'--','LineWidth',6)
plot(out3.time{1},out3.Voltage{1},'-.','LineWidth',6)
xlabel('Time [s]')
ylabel('Voltage [V]')
grid on
box on

author: lionsimbatoolbox

repo: LIONSIMBA


text:
%%***********************************************************
%% etp: Education testing problem.
%%
%% (dual problem) maximize     e'*d
%%                  subject to   B - diag(d) >= 0
%%                               d >= 0
%%
%% (primal problem) minimize     Tr B*X
%%                  subject to   X >= 0
%%                               diag(X) >= e
%%
%% Ref: M.T. Chu, J.W. Wright, IMA J. of Numerical Anal.,
%%      15 (1995), pp. 141--160.
%%-----------------------------------------------------------
%% [blk,Avec,C,b,X0,y0,Z0,objval,d] = etp(B,feas,solve);
%%
%% B = nxn positive definite.
%% feas  = 1 if want feasible starting point
%%       = 0 if otherwise.
%% solve = 0 just to initialize
%%       = 1 if want to solve the problem.
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%***********************************************************

   function [blk,Avec,C,b,X0,y0,Z0,objval,d] = etp(B,feas,solve);

   if nargin < 2; feas = 0; end;
   if nargin < 3; solve = 0; end;
   if (~isreal(B))
      error('only real B allowed');
   elseif (norm(B-B','fro') > 1e-13);
      error(' B must be symmetric');
   end;
%%
%% validate B
%%
   n = length(B);
   d = eig(B); d = real(d);
   if (min(d) < 0);
      error('B must be positive def');
   end;
%%
%%
   blk{1,1} = 's'; blk{1,2} = n;
   blk{2,1} = 'l'; blk{2,2} = n;
   b = ones(n,1);
   C{1,1} = B;
   C{2,1} = zeros(n,1);

   A = cell(2,n);
   for k = 1:n
       A{1,k} = sparse(k,k,1,n,n);
       A{2,k} = [zeros(k-1,1); -1; zeros(n-k,1)];
   end;

   Avec = svec(blk,A,ones(size(blk,1),1));
   if (feas == 1);
      y0 = 0.9*min(d)*ones(n,1);
      Z0 = ops(C,'-',Atyfun(blk,Avec,[],[],y0));
      X0{1,1} = 1.1*eye(n);
      X0{2,1} = 0.1*ones(n,1);
   elseif (feas == 0);
      [X0,y0,Z0] = infeaspt(blk,Avec,C,b);
   end;
   if (solve)
      [obj,X,y,Z] = sqlp(blk,Avec,C,b,[],X0,y0,Z0);
      objval = obj(2);
      d = y;
   else
      objval = []; d = [];
   end
%%===========================================================

repo: REST


text:
function geometry_test192 ( )

%*****************************************************************************80
%
%% TEST192 tests SPHERE_UNIT_SAMPLE_3D_2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2005
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;

  n_sample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST192\n' );
  fprintf ( 1, '  For the unit sphere in 3 dimensions:\n' );
  fprintf ( 1, '  SPHERE_UNIT_SAMPLE_3D_2 samples;\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  A few sample values: \n' );
  fprintf ( 1, '\n' );

  for i = 1 : 5
    [ x, seed ] = sphere_unit_sample_3d_2 ( seed );
    fprintf ( 1, '  %10f  %10f  %10f\n', x(1:dim_num) );
  end

  average(1:dim_num) = 0.0;

  for i = 1 : n_sample
    [ x, seed ] = sphere_unit_sample_3d_2 ( seed );
    average(1:dim_num) = average(1:dim_num) + x(1:dim_num);
  end

  average(1:dim_num) = average(1:dim_num) / n_sample;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now average the points, which should get a value\n' );
  fprintf ( 1, '  close to zero, and closer as N_SAMPLE increases.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Average:        %10f  %10f  %10f\n', average(1:dim_num) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now choose a random direction, sample the same\n' );
  fprintf ( 1, '  number of points, and compute the dot product with\n' );
  fprintf ( 1, '  the direction.\n' );
  fprintf ( 1, '  Take the absolute value of each dot product \n' );
  fprintf ( 1, '  and sum and average.\n' );

  for j = 1 : 5

    [ v, seed ] = sphere_unit_sample_3d_2 ( seed );

    dot_average = 0.0;

    for i = 1 : n_sample
      [ x, seed ] = sphere_unit_sample_3d_2 ( seed );
      dot_average = dot_average + abs ( x(1:dim_num) * v(1:dim_num)' );
    end

    dot_average = dot_average / n_sample;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  V:                %10f  %10f  %10f\n', v(1:dim_num) );
    fprintf ( 1, '  Average |(XdotV)| %10f\n', dot_average );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [W_RF,W_BB]=OMP_Combining(F_RF,F_BB)
global  Nr Ns Nrf H Vn Codebook_w
W_RF = [];
W_MMSE = ((F_BB' * F_RF' * H' * H * F_RF * F_BB + Vn * Ns * eye(Ns))^(-1) * F_BB' * F_RF' * H')';
Wres = W_MMSE;
n = 1 / Ns * H * F_RF * F_BB * F_BB' * F_RF' *H' + Vn * eye(Nr);
for i = 1 : Nrf
    y = Codebook_w' * n * Wres;
    k = find(diag(y * y')==max(diag(y * y')));
    W_RF(:,i) = Codebook_w(:,k);
    W_BB = (W_RF' * n * W_RF)^(-1) * W_RF' * n * W_MMSE;
    Wres = (W_MMSE - W_RF * W_BB) / norm(W_MMSE - W_RF * W_BB,'fro');
end

author: Zzhaoxingyu

repo: hybrid-beamforming-for-three-scenes


text:
function [x, infos] = convex_mu_nmf(V, rank, in_options)
% (Kernel) Convex multiplicative update for non-negative matrix factorization ((Kernel-)Convex-MU-NMF).
%
% The problem of interest is defined as
%
%  (standard)   min || VWH - V ||_F^2,
%               where
%               {W, H} >= 0, and
%
%  (kernel)     min || K(V,V) WH - K(V,V) ||_F^2,
%               where
%               {W, H} >= 0, and K(V,V) is a kernel matrix.
%
%
% Given a V with mixed-signs, factorized non-negative matrices {W, H} are calculated.
%
%
% Inputs:
%       matrix      V
%       rank        rank
%       in_options
%           sub_mode: 'std' (default) or 'kernel'
%           kernel  : kernel type (rbf (default), polynomial, linear, sigmoid)
%
% Output:
%       x           solution of x
%       infos       information
%
% References:
%       C. Ding, T. Li, and M.I. Jordan,
%       "Convex and semi-nonnegative matrix factorizations,"
%       IEEE Transations on Pattern Analysis and Machine Intelligence,
%       vol. 32, no. 1, pp. 45-55, 2010.
%
%       T. Li and C. Ding,
%       "The Relationships Among Various Nonnegative Matrix Factorization Methods for Clustering,"
%       International Conference on Data Mining,
%       2006.
%
%       Y. Li and A. Ngom,
%       "A New Kernel Non-Negative Matrix Factorization and Its Application in Microarray Data Analysis,"
%       CIBCB,
%       2012.
%
%
% This file is part of NMFLibrary
%
% This file has been ported from
% convexnmfrule.m and kernelconvexnmf.m at https://sites.google.com/site/nmftool/home/source-code
% by Yifeng Li.
%
%%%%
% Copyright (C) <2012>  <Yifeng Li>
%
% 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
% 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/>.
%
% Contact Information:
% Yifeng Li
% University of Windsor
% li11112c@uwindsor.ca; yifeng.li.cn@gmail.com
% May 01, 2011
%%%%
%
%
% This file was originally created by Graham Grindla.
%
% 2010-01-14 Graham Grindlay (grindlay@ee.columbia.edu)
%
% Copyright (C) 2008-2028 Graham Grindlay (grindlay@ee.columbia.edu)
%
% 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/>.
%
%
% Ported by H.Kasai on June 30, 2022
%
% Change log:
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.sub_mode = 'std';
    local_options.kernel ='rbf';
    local_options.kernel_param = [];
    local_options.special_nmf_cost = @(V, W, H, R, options) convex_mu_nmf_cost_func(V, W, H, R, options);
    local_options.special_init_factors = @(V, rank, wh_flag, r_flag, options) convex_mu_initialize(V, rank, wh_flag, r_flag, options);

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    W = init_factors.W;
    H = init_factors.H;

    % initialize
    epoch = 0;
    grad_calc_count = 0;

    % initialize for this algorithm
    if strcmp(options.sub_mode, 'std')
        Ak = V' * V;
        X = V;
        options.name = 'std';
    elseif strcmp(options.sub_mode, 'kernel')
        Ak = computeKernelMatrix(V, V, options);
        X = Ak;
        options.name = sprintf('kernel-%s', options.kernel);
    else
        error('Invalid sub_mode')
    end
    method_name = sprintf('Convex-MU (%s)', options.name);

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    Ap = (abs(Ak) + Ak) ./ 2;
    An = (abs(Ak) - Ak) ./ 2;

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(X, W, H, [], options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('Convex-MU (%s): Epoch = 0000, cost = %.16e, optgap = %.4e\n', options.name, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        ApW = Ap * W;
        AnW = An * W;
        WH  = W * H;

        % update H
        H = H .* sqrt((ApW' + AnW' * WH) ./ (AnW' + ApW' * WH));
        H = max(H, eps);
        HHt = H * H';

        % update W
        W = W .* sqrt((Ap * H' + AnW * HHt) ./ (An * H' + ApW * HHt));
        W = max(W, eps);

        % measure gradient calc count
        grad_calc_count = grad_calc_count + m*n;

        % measure elapsed time
        elapsed_time = toc(start_time);

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(X, W, H, [], options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = W;
    x.H = H;

end

function res = convex_mu_nmf_cost_func(V, W, H, R, options)

    res = norm(V - V * W * H, 'fro');

end

function [init_factors, init_factors_opts] = convex_mu_initialize(V, rank, wh_flag, r_flag, options)

    init_factors_opts = [];

    [~, n] = size(V);

    if wh_flag
        H = rand(rank, n);
        W = H' * diag(1./sum(H,2)');
    else
        W = options.x_init.W;
        H = options.x_init.H;
    end

    init_factors.W = W;
    init_factors.H = H;
    init_factors.R = zeros(size(V));

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function y = norm_pdf(x,mu,sigma)
%NORM_PDF Normal probability density function (pdf).
%
%   Y = NORMPDF(X,MU,SIGMA) Returns the normal pdf with
%   mean, MU, and standard deviation, SIGMA, at the values in X.
%
%   The size of X is the common size of the input arguments. A scalar input
%   functions as a constant matrix of the same size as the other inputs.
%
%   Default values for MU and SIGMA are 0 and 1 respectively.

% Copyright (c) 1998-2004 Aki Vehtari

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

if nargin < 3,
  sigma = 1;
end

if nargin < 2;
  mu = 0;
end

if nargin < 1,
  error('Requires at least one input argument.');
end

y = exp(-0.5 * ((x-mu)./sigma).^2 -log(sigma) -log(2*pi)/2);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function r = subsref(p,s)
%SUBSREF      Implements subscripted references for polynomials
%
%  r = p(i)
%
%For univariate polynomials p, p(i) is the coefficient of x^i, 0<=i<=degree(p).
%  For i>degree(p), p(i):=0. Similarly, p(i:j) is the vector of coefficients
%  [ p(i) p(i+1) ... p(j) ], or, p(:) is the (row) vector of all coefficients
%  of p, the same as vector(p) for univariate polynomials.
%
%For multivariate polynomials p in k variables x_1..x_k, p(i) is the
%  coefficient of x_1^i, i.e. a polynomial in k-1 variables.
%  Similarly, p(i,[],j) or p(i,:,j) is the coefficient polynomial of
%  x1^i*x3^j. Indices i,j,... must be single indices, no range.
%Note that access to coefficients refers to the current order p.v of
%  variables of p. To change this order, see permvars.
%For example, for a polynomial in three variables p.v={'x','y','z'},
%  p(1,0,3) is the coefficient of x*z^3 (a constant), where p(1,[],3)
%  if the coefficient of x*z^3, a univariate polynomial in y.
%
%Polynomial evaluation is denoted by p{x} or p{x1,...,xn}, computing the value
%  of p at x. This is the same as polyval(p,x) of polyval(p,x1,...,xn).
%For univariate polynomials, x may be a vector or matrix yielding the vector
%  or matrix of polynomial values evaluated at the corresponding coefficients.
%For multivariate polynomials, x is a vector of values of the variables.
%  For x being a matrix, the result is the (column) vector of p{x(i,:)}.
%
%Moreover, p.mid, p.rad, p.inf, p.sup give access to the midpoint, radius,
%  infimum and supremum of p, respectively.
%
%Finally, p.e, p.c and p.v give access to the arrays of exponents, coefficients
%  and variables of p, such that polynom(p.e,p.c,p.v) is again p. Single variables
%  are accessed by p.v{i}.
%
%In the univariate case,  p == polynom(p.c,p.v)  [ degree is length(p.c)+1 ].
%In the multivariate case,  p == polynom(p.e,p.c,p.v)  [ p.e is (sparse) exponent set ].
%

% written  11/20/97     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  while 1
    if ~isa(p,'polynom')
      r = subsref(p,s(1));
    elseif strcmp(s(1).type,'()')     % index reference p(i)
      if size(p.e,2)==1               % univariate polynomial
        if length(s(1).subs)>1
          error('invalid call: more than one index')
        end
        if isequal(s(1).subs{1},':')
          r = p.c;
        else
          n1 = length(p.c);
          index = n1-s(1).subs{1};    % index vector, p(0) constant term, etc.
          if ~isreal(index) | ~isequal(index,round(index))
            error('index must be integer')
          end
          if any(index>n1)            % index negative
            error('index out of range')
          end
          if any(index<1)             % index greater than degree
            indexgt0 = ( index>0 );
            r = typeadj( zeros(1,length(index)) , typeof(p.c) );
            if any(indexgt0)
              r(indexgt0) = p.c(index(indexgt0));
            end
          else
            r = p.c(index);
          end
        end
      else                            % multivariate polynomial
        k = size(p.e,2);              % number of variables
        if length(s(1).subs)>k
          error('too many indices')
        end
        index = logical(zeros(1,k));
        exponents = [];
        for i=1:length(s(1).subs)
          if isempty(s(1).subs{i}) | isequal(s(1).subs{i},':')
            index(i) = 0;
          else
            index(i) = 1;
            exponents = [ exponents s(1).subs{i} ];
          end
        end
        if isempty(exponents)
          r = p;
          return
        end
        I = all( ( p.e(:,index) == repmat(exponents,size(p.e,1),1) ) , 2 );
        if ~any(I)                    % coefficients do not occur
          r = typeadj( 0 , typeof(p.c) );
          return
        end
        r.e = p.e(I,~index);
        if isempty(r.e)               % coefficient is single constant
          r = p.c(I);
          return
        end
        r.c = p.c(I);
        r.v = p.v(~index);
        if size(r.e,2)==1             % coefficient polynomial univariate
          if iscell(r.v)              % be sure r.v is char for univ. pol.
            r.v = r.v{1};
          end
          n = max(r.e);               % r.e is vector
          c = r.c;
          r.c = typeadj( zeros(1,n+1) , typeof(c) );
          r.c(n+1-r.e) = c;
          r.e = n;
          if r.e==0
            r = r.c;
            return
          end
        else
          r = normalize(r);
        end
        r = class(r,'polynom');
      end
    elseif strcmp(s(1).type,'{}')     % polynomial evaluation p{x}
      r = polyval(p,s(1).subs{:});
    elseif strcmp(s(1).type,'.')      % polynomial access to inf, sup, ...
      if     strcmp(s(1).subs,'mid'), r = mid(p);
      elseif strcmp(s(1).subs,'rad'), r = rad(p);
      elseif strcmp(s(1).subs,'inf'), r = inf(p);
      elseif strcmp(s(1).subs,'sup'), r = sup(p);
      elseif strcmp(s(1).subs,'e')
        if iscell(p.v)                  % multivariate
          r = flipud(sortrows(p.e));
        else
          r = p.e;                      % univariate
        end
      elseif strcmp(s(1).subs,'c')
        if iscell(p.v)                  % multivariate
          [ dummy index ] = sortrows(p.e); r = flipud(p.c(index));
        else
          r = p.c;
        end
      elseif strcmp(s(1).subs,'v'), r = p.v;
      else
        error('invalid reference for polynomial')
      end
    else
      error('invalid index reference for polynomial')
    end
    if length(s)==1
      return
    end
    s = s(2:end);
    p = r;
  end

repo: OpenMAS


text:
function edgeBoxesSweeps()
% Parameter sweeps for Edges Boxes object proposals.
%
% Running the parameter sweeps requires altering internal flags.
% The sweeps are not well documented, use at your own discretion.
%
% Structured Edge Detection Toolbox      Version 3.01
% Code written by Piotr Dollar and Larry Zitnick, 2014.
% Licensed under the MSR-LA Full Rights License [see license.txt]

% define parameter sweeps
rt = 'D:\code\research\edges\';
expNms = {'alpha','beta','minScore','edgeMinMag','edgeMergeThr',...
  'clusterMinMag','maxAspectRatio','minBoxArea','gamma','kappa'};
expNms=expNms(1:end); opts=createExp(rt,expNms); maxn=inf;

% run training and testing jobs
jobs = createJobs(rt,opts,maxn);
tic, for i=1:length(jobs), edgeBoxes(jobs{i}{:}); end; toc

% plot all results
for e=1:length(expNms)
  eval={'data',boxesData('split','val'),'names',{opts{e}.name},...
    'resDir',[rt 'boxes/sweeps/'],'maxn',maxn,'fName',expNms{e}};
  boxesEval(eval{:},'thrs',.7);
  boxesEval(eval{:},'thrs',.5:.05:1,'cnts',1000);
end

end

function jobs = createJobs( rt, opts, maxn )
% create jobs
M='models/forest/modelBsds'; M=load(M); M=M.model; M.opts.nThreads=1;
data=boxesData('split','val'); fs=data.imgs(1:min(end,maxn));
opts=[opts{:}]; N=length(opts); jobs=cell(1,N); D=zeros(1,N);
for e=1:N, opts(e).name=[rt 'boxes/sweeps/' opts(e).name '-val.mat']; end
for e=1:N, D(e)=exist(opts(e).name,'file')==2; jobs{e}={fs,M,opts(e)}; end
[~,K]=unique({opts.name},'stable'); D=D(K); jobs=jobs(K); jobs=jobs(~D);
fprintf('nJobs = %i\n',length(jobs));
end

function opts = createExp( rt, expNm )

% if expNm is a cell, call recursively and return
if( iscell(expNm) )
  N=length(expNm); opts=cell(1,N);
  for e=1:N, opts{e}=createExp(rt,expNm{e}); end; return;
end

% setup opts
opts=edgeBoxes(); opts.minScore=0;
N=100; optsDefault=opts; opts=opts(ones(1,N));
switch expNm
  case 'alpha'
    vs=45:5:75; N=length(vs);
    for e=1:N, opts(e).alpha=vs(e)/100; end
  case 'beta'
    vs=60:5:90; N=length(vs);
    for e=1:N, opts(e).beta=vs(e)/100; end
  case 'minScore'
    vs=[0 5 10 25 50 100]; N=length(vs);
    for e=1:N, opts(e).minScore=vs(e)/1000; end
  case 'edgeMinMag'
    vs=[0 50 100 200 400]; N=length(vs);
    for e=1:N, opts(e).edgeMinMag=vs(e)/1000; end
  case 'edgeMergeThr'
    vs=[25 50 100 200 400]; N=length(vs);
    for e=1:N, opts(e).edgeMergeThr=vs(e)/100; end
  case 'clusterMinMag'
    vs=[0 50 100 200 400]; N=length(vs);
    for e=1:N, opts(e).clusterMinMag=vs(e)/100; end
  case 'maxAspectRatio'
    vs=1:5; N=length(vs);
    for e=1:N, opts(e).maxAspectRatio=vs(e); end
  case 'minBoxArea'
    vs=[100 250 500 1000 2500 5000]; N=length(vs);
    for e=1:N, opts(e).minBoxArea=vs(e); end
  case 'gamma'
    vs=[25 50 100 200 400 10000]; N=length(vs);
    for e=1:N, opts(e).gamma=vs(e)/100; end
  case 'kappa'
    vs=50:25:200; N=length(vs);
    for e=1:N, opts(e).kappa=vs(e)/100; end
  otherwise, error('invalid exp: %s',expNm);
end

% produce final set of opts and find default opts
O=1:N; opts=opts(O); d=0;
for e=1:N, if(isequal(optsDefault,opts(e))), d=e; break; end; end
if(d==0), disp(expNm); assert(false); end; opts(d).name='Default';
for e=1:N, if(e~=d), opts(e).name=[expNm int2str2(vs(e),5)]; end; end

end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function o2 = montage(obj, varargin)
% This function displays an atlas object on a standard slice montage
% Call with 'nosymmetric' or 'colors' ... to match colors with wedge plots
%
% :Usage:
% ::
%
%    montage(obj, [optional inputs])
%
% - takes all optional inputs to canlab_results_fmridisplay and addblobs method
%
% :Input:
%
%   **obj:**
%        a region object
%
% :Optional Inputs:  - see help canlab_results_fmridisplay
%
%   **o2***
%        An existing fmridisplay object, with no keyword strings
%
%   **symmetric** [default]
%       Mirror left/right blobs with same colors
%       See match_colors_left_right
%
%   **nosymmetric**
%       Standard color map, no L/R color-matching for symmetry
%       Call with 'nosymmetric' to match colors with wedge plots
%
%   **'noblobs':**
%        do not display blobs
%
%   **'nooutline':**
%        do not display blob outlines
%
%   **'addmontages':**
%        when entering existing fmridisplay obj, add new montages
%
%   **'noremove':**
%        do not remove current blobs when adding new ones
%
%   **'outlinecolor:**
%        followed by new outline color
%
%   **'splitcolor':**
%        followed by 4-cell new split colormap colors (help fmridisplay or edit code for defaults as example)
%
%   **'montagetype':**
%        'full' for full montages of axial and sagg slices.
%
%        'full hcp' for full montage, but with surfaces and volumes from
%        HCP data
%
%        'compact' [default] for single-figure parasagittal and axials slices.
%
%        'compact2': like 'compact', but fewer axial slices.
%
%        'multirow': followed by number of rows
%           e.g., o2 = canlab_results_fmridisplay([], 'multirow', 2);
%
%   **'noverbose':**
%        suppress verbose output, good for scripts/publish to html, etc.
%
%   **'overlay':**
%        specify anatomical image for montage (not surfaces), followed by
%        image name
%        e.g., o2 = canlab_results_fmridisplay([], 'overlay', 'icbm152_2009_symmetric_for_underlay.img')';
%
%         The default brain for overlays is based on Keuken et al. 2014
%         For legacy SPM8 single subject, enter as arguments:
%         'overlay', which('SPM8_colin27T1_seg.img')
%
% Other inputs to addblobs (fmridisplay method) are allowed, e.g., 'cmaprange', [-2 2], 'trans'
%
% See help fmridisplay and region/montage
% e.g., 'color', [1 0 0]
%
%
% Examples:
%
% o2 = canlab_results_fmridisplay([], 'noverbose');
% o2 = montage(r, o2);      % symmetric colors left/right
% o2 = removeblobs(o2);
% o2 = montage(r, o2, 'map');

% PROGRAMMERS' NOTES:
% Tor Wager, Feb 2018.
% This function simply invokes the region montage method.

r = atlas2region(obj);

o2 = montage(r, varargin{:});

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore