text:
function [image_roidb, bbox_means, bbox_stds] = rfcn_prepare_image_roidb(conf, imdbs, roidbs, bbox_means, bbox_stds)
% [image_roidb] = rfcn_prepare_image_roidb(conf, imdbs, roidbs, cache_img)
% --------------------------------------------------------
% R-FCN implementation
% Modified from MATLAB Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2016, Jifeng Dai
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

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

    if ~iscell(imdbs)
        imdbs = {imdbs};
        roidbs = {roidbs};
    end

    imdbs = imdbs(:);
    roidbs = roidbs(:);

    image_roidb = cellfun(@(x, y) ... // @(imdbs, roidbs)
                    arrayfun(@(z) ... //@([1:length(x.image_ids)])
                        struct('image_path', x.image_at(z), 'image_id', x.image_ids{z}, 'im_size', x.sizes(z, :), 'imdb_name', x.name, ...
                        'overlap', y.rois(z).overlap, 'boxes', y.rois(z).boxes, 'class', y.rois(z).class, 'image', [], 'bbox_targets', []), ...
                        [1:length(x.image_ids)]', 'UniformOutput', true), imdbs, roidbs, 'UniformOutput', false);

    image_roidb = cat(1, image_roidb{:});

    % enhance roidb to contain bounding-box regression targets
    [image_roidb, bbox_means, bbox_stds] = append_bbox_regression_targets(conf, image_roidb, bbox_means, bbox_stds);
end

function [image_roidb, means, stds] = append_bbox_regression_targets(conf, image_roidb, means, stds)
    % means and stds -- (k+1) * 4, include background class

    num_images = length(image_roidb);
    % Infer number of classes from the number of columns in gt_overlaps
    if conf.bbox_class_agnostic
        num_classes = 1;
    else
        num_classes = size(image_roidb(1).overlap, 2);
    end

    valid_imgs = true(num_images, 1);
    for i = 1:num_images
       rois = image_roidb(i).boxes;
       [image_roidb(i).bbox_targets, valid_imgs(i)] = ...
           compute_targets(conf, rois, image_roidb(i).overlap);
    end
    if ~all(valid_imgs)
        image_roidb = image_roidb(valid_imgs);
        num_images = length(image_roidb);
        fprintf('Warning: rfcn_prepare_image_roidb: filter out %d images, which contains zero valid samples\n', sum(~valid_imgs));
    end

    if ~(exist('means', 'var') && ~isempty(means) && exist('stds', 'var') && ~isempty(stds))
        % Compute values needed for means and stds
        % var(x) = E(x^2) - E(x)^2
        class_counts = zeros(num_classes, 1) + eps;
        sums = zeros(num_classes, 4);
        squared_sums = zeros(num_classes, 4);
        for i = 1:num_images
           targets = image_roidb(i).bbox_targets;
           for cls = 1:num_classes
              cls_inds = find(targets(:, 1) == cls);
              if ~isempty(cls_inds)
                 class_counts(cls) = class_counts(cls) + length(cls_inds);
                 sums(cls, :) = sums(cls, :) + sum(targets(cls_inds, 2:end), 1);
                 squared_sums(cls, :) = squared_sums(cls, :) + sum(targets(cls_inds, 2:end).^2, 1);
              end
           end
        end

        means = bsxfun(@rdivide, sums, class_counts);
        stds = (bsxfun(@minus, bsxfun(@rdivide, squared_sums, class_counts), means.^2)).^0.5;

        % add background class
        means = [0, 0, 0, 0; means];
        stds = [0, 0, 0, 0; stds];
    end

    % Normalize targets
    for i = 1:num_images
        targets = image_roidb(i).bbox_targets;
        for cls = 1:num_classes
            cls_inds = find(targets(:, 1) == cls);
            if ~isempty(cls_inds)
                image_roidb(i).bbox_targets(cls_inds, 2:end) = ...
                    bsxfun(@minus, image_roidb(i).bbox_targets(cls_inds, 2:end), means(cls+1, :));
                image_roidb(i).bbox_targets(cls_inds, 2:end) = ...
                    bsxfun(@rdivide, image_roidb(i).bbox_targets(cls_inds, 2:end), stds(cls+1, :));
            end
        end
    end
end

function [bbox_targets, is_valid] = compute_targets(conf, rois, overlap)

    overlap = full(overlap);

    [max_overlaps, max_labels] = max(overlap, [], 2);

    % ensure ROIs are floats
    rois = single(rois);

    bbox_targets = zeros(size(rois, 1), 5, 'single');

    % Indices of ground-truth ROIs
    gt_inds = find(max_overlaps == 1);

    if ~isempty(gt_inds)
        % Indices of examples for which we try to make predictions
        ex_inds = find(max_overlaps >= conf.bbox_thresh);

        % Get IoU overlap between each ex ROI and gt ROI
        ex_gt_overlaps = boxoverlap(rois(ex_inds, :), rois(gt_inds, :));

        assert(all(abs(max(ex_gt_overlaps, [], 2) - max_overlaps(ex_inds)) < 10^-4));

        % Find which gt ROI each ex ROI has max overlap with:
        % this will be the ex ROI's gt target
        [~, gt_assignment] = max(ex_gt_overlaps, [], 2);
        gt_rois = rois(gt_inds(gt_assignment), :);
        ex_rois = rois(ex_inds, :);

        [regression_label] = rfcn_bbox_transform(ex_rois, gt_rois);

        if conf.bbox_class_agnostic
            bbox_targets(ex_inds, :) = [max_labels(ex_inds)>0, regression_label];
        else
            bbox_targets(ex_inds, :) = [max_labels(ex_inds), regression_label];
        end
    end

    % Select foreground ROIs as those with >= fg_thresh overlap
    is_fg = max_overlaps >= conf.fg_thresh;
    % Select background ROIs as those within [bg_thresh_lo, bg_thresh_hi)
    is_bg = max_overlaps < conf.bg_thresh_hi & max_overlaps >= conf.bg_thresh_lo;

    % check if there is any fg or bg sample. If no, filter out this image
    is_valid = true;
    if ~any(is_fg | is_bg)
        is_valid = false;
    end
end

repo: R-FCN


text:
function ndx = merge(param1,param2)
% ndx = merge(ndx_obj_array)
% ndx = merge(ndx1,ndx2)
% Merges Ndx objects.  The function can be called in two ways: with
% two ndxs or with an object array of ndxs.
% Inputs:
%   Either:
%   ndx1: An Ndx object.
%   ndx2: Another Ndx object.
%   Or:
%   ndx_obj_array: An array of Ndx objects.
% Outputs:
%   ndx: An Ndx object that contains the information from the input ndxs.

assert(nargin<3)
if nargin==2
    ndx1 = param1;
    ndx2 = param2;
    assert(isa(ndx1,'Ndx'))
    assert(isa(ndx2,'Ndx'))
    assert(ndx1.validate())
    assert(ndx2.validate())

    ndx_obj_arr = Ndx.empty(2,0);
    ndx_obj_arr(1) = ndx1;
    ndx_obj_arr(2) = ndx2;

    ndx = Ndx.merge(ndx_obj_arr);
else
    % the output ndx must have all models and segment in the input
    % ndxs (only once).  A trial in any ndx becomes a trial in the
    % output ndx.
    assert(nargin==1)
    obj_array = param1;
    assert(isa(obj_array,'Ndx'))
    numndxs = length(obj_array);
    ndx = Ndx();
    for ii=1:numndxs
	ndx_new = Ndx();
	ndx1 = ndx;
	ndx2 = obj_array(ii);

	% create new ndx with empty masks
	ndx_new.modelset = union(ndx1.modelset, ndx2.modelset);
	ndx_new.segset = union(ndx1.segset, ndx2.segset);

	% expand ndx1 mask
	trials_1 = false(length(ndx_new.modelset), length(ndx_new.segset));
	[dummy,model_index_a,model_index_b] = intersect(ndx_new.modelset,ndx1.modelset);
	[dummy,seg_index_a,seg_index_b] = intersect(ndx_new.segset,ndx1.segset);
	trials_1(model_index_a,seg_index_a) = ndx1.trialmask(model_index_b,seg_index_b);

	% expand ndx2 mask
	trials_2 = false(length(ndx_new.modelset), length(ndx_new.segset));
	[dummy,model_index_a,model_index_b] = intersect(ndx_new.modelset,ndx2.modelset);
	[dummy,seg_index_a,seg_index_b] = intersect(ndx_new.segset,ndx2.segset);
	trials_2(model_index_a,seg_index_a) = ndx2.trialmask(model_index_b,seg_index_b);

	% merge masks
	trials = trials_1 | trials_2;

	% build new ndx
	ndx_new.trialmask = trials;
	ndx = ndx_new;
    end
end

assert(ndx.validate())

author: nesl


text:
%% Author: epokh
%% Website: www.epokh.org/drupy
%% This software is under GPL

%%This function plan a trajectory passing in N
%% points using N-1 3th order polynomial:
%% the user must choose: intial,final and intermediate points
%% timePoints is a matrix [2,N] where first row is the time,second row
%% is the position points
%% the velocities this time are computed as a linear interpolation between
%% the position points
%% Output: a matrix [N-1,4] whose columns are the polynomail coefficients

function [TJ]=trajectory3VelComputed(timePoints)
N=length(timePoints);

if(N<2)
fprintf('Error: you defined only 1 point');
return;
end
TJ=[];

for k=1:N-1
    tk=timePoints(1,k);
    tk1=timePoints(1,k+1);
    qk=timePoints(2,k);
    qk1=timePoints(2,k+1);
    if(k==1)
        velk=0;
    else
        velk=(qk-timePoints(2,k-1))/(tk-timePoints(1,k-1));
    end

    if(k==N-1)
        velk1=0;
    else
       velk1=( timePoints(2,k+1)- timePoints(2,k))/(timePoints(1,k+1)-timePoints(1,k));
    end

    if(sign(velk)~=sign(velk1))
        velk=0;
        velk1=0;
    else
        velk=0;
        velk1=(velk+velk1)/2;

    end

    coeffs=pol3interpol(tk,tk1,qk,qk1,velk,velk1);
    TJ=[TJ;coeffs'];
end

end

author: Sable

repo: mcbench-benchmarks


text:
% rsadjust() - adjust l-values (Ramberg-Schmeiser distribution)
%                with respect to signal mean and variance
%
% Usage: p = rsadjust(l3, l4, m, var, skew)
%
% Input:
%   l3   - value lambda3 for Ramberg-Schmeiser distribution
%   l4   - value lambda4 for Ramberg-Schmeiser distribution
%   m    - mean of the signal distribution
%   var  - variance of the signal distribution
%   skew - skewness of the signal distribution (only the sign of
%          this parameter is used).
%
% Output:
%   l1  - value lambda3 for Ramberg-Schmeiser distribution
%   l2  - value lambda4 for Ramberg-Schmeiser distribution
%   l3  - value lambda3 for Ramberg-Schmeiser distribution (copy
%         from input)
%   l4  - value lambda4 for Ramberg-Schmeiser distribution (copy
%         from input)
%
% Author: Arnaud Delorme, SCCN, 2003
%
% See also: rsfit(), rsget()
%
% Reference: Ramberg, J.S., Tadikamalla, P.R., Dudewicz E.J., Mykkytka, E.F.
%            A probability distribution and its uses in fitting data.
%            Technimetrics, 1979, 21: 201-214.

% Copyright (C) 2003 Arnaud Delorme, SCCN, arno@salk.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 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [l1,l2,l3,l4] = rsadjust( l3, l4, mu, sigma2, m3);

    % swap l3 and l4 for negative skewness
    % ------------------------------------
    if m3 < 0
        ltmp = l4;
        l4   = l3;
        l3   = ltmp;
    end;

    A = 1/(1 + l3) - 1/(1 + l4);
    B = 1/(1 + 2*l3) + 1/(1 + 2*l4) - 2*beta(1+l3, 1+l4);
    C = 1/(1 + 3*l3) - 1/(1 + 3*l4) ...
             - 3*beta(1+2*l3, 1+l4) + 3*beta(1+l3, 1+2*l4);

    % compute l2 (and its sign)
    % ------------------------
    l2 = sqrt( (B-A^2)/sigma2 );
    if m3 == 0, m3 = -0.000000000000001; end;
    if (m3*(C - 2*A*B + 2*A^3)) < 0, l2 = -l2; end;
    %l22 = ((C - 2*A*B + 2*A^3)/m3)^(1/3) % also equal to l2

    % compute l1
    % ----------
    l1 = mu - A/l2;

    return;

    % fitting table 1 of
    % ---------------
    [l1 l2]  = pdffitsolv2(-.0187,-.0388, 0, 1, 1)
    [l1 l2]  = pdffitsolv2(-.1359,-.1359, 0, 1, 0)

    % sign problem
    [l1 l2]  = pdffitsolv2(1.4501,1.4501, 0, 1)

    % numerical problem for l1
    [l1 l2]  = pdffitsolv2(-.00000407,-.001076, 0, 1, 2)
    [l1 l2]  = pdffitsolv2(0,-.000580, 0, 1, 2)

author: PatternRecognition

repo: OpenBMI


text:
% Script to test actuator at current parameter values using frestimate
% Copyright 2010 MathWorks, Inc.

mdl = 'actuator_freq_resp';

amplitude_100 = 0.2;
frequency_100 = 43;
amplitude_40  = amplitude_100*0.4;
frequency_40  = 55;

ios = getlinio(mdl);
in_2pts = frest.Sinestream(...
   'Frequency',[20, 20.1, 43, 57], ...
   'Amplitude', [amplitude_100, amplitude_40, amplitude_100, amplitude_40], ...
   'SimulationOrder', 'Sequential', ... %Set to OneAtATime to restart simuation for each frequency
   'FreqUnits', 'Hz');

sys = frestimate(mdl,ios,in_2pts);
R = sys.ResponseData(:);
Lag = angle(R)/pi*180;
%disp('');
%disp(['Desired Phase: -59 -50 -90 -90']);
%disp(['Actual Phase:  ' sprintf('%3.0f %3.0f %3.0f %3.0f',Lag(1),Lag(2),Lag(3),Lag(4))]);


text:
function jtwgft = gsp_jtwgft(G, w, s, param)
%GSP_JTWGFT Compute the joint vertex-time windowed graph fourier transform
%   Usage:  jtwgft = gsp_jtwgft(G, w, s)
%           jtwgft = gsp_jtwgft(G, w, s, param)
%
%   Input parameters:
%         G          : Graph structure
%         w          : windows
%         s          : vertex-time signal
%         param      : structure of optional parameters
%   Output parameters:
%         jtwgft     : Joint vertex-time windowed graph Fourier transform
%
%   Additional parameters
%   ---------------------
%   * *param.boundary*    : 'periodic' or 'reflecting' (default 'periodic')
%   * *param.a*           : hop size in time (default 1)
%   * *param.M*           : number of frequencies (default size(s,2))
%

% Author : Nathanael Perraudin
% Date   : 30 April 2016

if nargin<4
    param=struct;
end
if ~isfield(param,'boundary'), param.boundary = 'periodic'; end
if ~isfield(param,'a'), param.a = 1; end
if ~isfield(param,'M'), param.M = size(s,2); end

switch param.boundary
    case 'periodic'
        % 1) perform DGT
        S = dgt(transpose(s),w,param.a,param.M);
        % 2) Reshape
        S = reshape(S,[],G.N);
        % 3) GFT
        Shat = transpose(gsp_gft(G,transpose(S)));
        % 4) Reshape
        jtwgft = reshape(Shat,param.M,[],G.N);
    case 'reflecting'
        % 1) perform DGT
        S = wmdct(transpose(s),w,param.M,size(s,2));
        % 2) Reshape
        S = reshape(S,[],G.N);
        % 3) GFT
        Shat = transpose(gsp_gft(G,transpose(S)));
        % 4) Reshape
        jtwgft = reshape(Shat,param.M,[],G.N);
    otherwise
        error('Unknown boundary condition');
end

end

repo: gspbox


text:
function [S,a,t] = normalizeS(S)

t = mean(S,2);
S = S - t*ones(1,size(S,2));
a = mean(std(S, 1, 2));
S = S / a;

author: geopavlakos


text:
function [ t ] = net_establishConnection( ip )
%% This function is used to establish a connection with the server on KUKA

%% Syntax:
% function [ t ] = net_establishConnection( ip )
% This function is used to establish a connection with the server on KUKA.

%% Arregument:
% ip: a string variable, the IP of the robot, example '172.31.1.147',
% To figure out your robot IP, do the following:
% from the teach pendant, click on 'Station' button, then the 'information'
% button, from the information page check out your KUKA iiwa IP.

%% Return value:
% t: is a TCP/IP object, this object is empty when the connection could not be established.

% Copy right, Mohammad SAFEEA, 3rd of May 2017

t=[];
if ~exist('ip','var') || isempty(ip)
  warning('Connection could not be establised');
  return;
else
    try
            port=30001;
            % Connect to server
            t= tcpip(ip, port, 'NetworkRole', 'client');
            t.TransferDelay='off';
            t.Timeout=3.0;
            fopen(t);
            pause(1);
    catch exception
       warning('Connection could not be establised');
    end

end

end

repo: KST-Kuka-Sunrise-Toolbox


text:
classdef C10MOP5 < EvoXBenchProblem
% <multi> <real> <large/none> <expensive/none>

%------------------------------- Reference --------------------------------
% Z. Lu, R. Cheng, Y. Jin, K. C. Tan, and K. Deb, Neural architecture
% search as multiobjective optimization benchmarks: Problem formulation and
% performance assessment, IEEE Transactions on Evolutionary Computation,
% 2023.
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            config.name = 'nb201';
            config.args.fidelity = 200;
            config.args.objs = ['err&params&flops&edgegpu_latency&' ...
                'edgegpu_energy'];
            config.args.dataset = 'cifar10';
            config.args.normalized_objectives = true;
            obj.Setting@EvoXBenchProblem(config);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function determ = pascal1_determinant ( n )

%*****************************************************************************80
%
%% PASCAL1_DETERMINANT returns the determinant of the PASCAL1 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM, the determinant.
%
  determ = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function fs_retinotopicTemplate(subject, out_path, subjects_dir)
%
% function fs_retinotopicTemplate(subject, out_path, subjects_dir)
%
% Parameters
% ----------
% subject : str
%    The subject ID in the Freesurfer context
% out_path : str
%    Full path to the location where output files will be generated
% subjects_dir : str (optional)
%    Full path to the Freesurfer SUBJECTS_DIR. Per default this will be set
%    to the environment-defined $SUBJECTS_DIR variable
%
% Notes
% -----
% Based on the work of Benson et al. (2012). The requirements
% for this to work are:
% 1. The data should have been processed using the cross-sectional stream
% in FreeSurfer (i.e. recon-all -s subjid -autorecon-all). Preferably using
% Freesurfer version 5.1
% 2. The scripts (surfreg) and atlas (fsaverage_sym), which are not part of
% the standard FreeSurfer 5.1 distribution, should be installed
% 3. surfreg should be installed in $FREESURFER_HOME/bin/
% 4. fsaverage_sym should be copied to the FreeSufer subject data directory
%
% The 'fsaverage_sym' subject used to be abvailable at:
%
% ftp://surfer.nmr.mgh.harvard.edu/transfer/outgoing/flat/greve/
%
% Until it returns to this location, you can copy it from, or symlink to:
%
% /home/arokem/data/freesurfer_subjects/fsaverage_sym/
%
%
% References
% ----------
% NC Benson, OH Butt, R Datta, PD Radoeva, DH Brainard, GK Aguirre (2012)
% The retinotopic organization of striate cortex is well predicted by
% surface topology. Current Biology 22: 2081-2085.
%
% Greve, Douglas N., Lise Van der Haegen, Qing Cai, Steven Stufflebeam, Mert
% R. Sabuncu, Bruce Fischl, and Marc Bysbaert. "A surface-based analysis of
% language lateralization and cortical asymmetry." (2013). Journal of
% Cognitive Neuroscience. In press.

% To accomodate user set SUBJECTS_DIR system variable:
if notDefined('subjects_dir')
    subjects_dir = getenv('SUBJECTS_DIR');
else
    syscall(sprintf('export SUBJECTS_DIR=%s', subjects_dir));
end

maps = {'ecc', 'pol'};
hemis = {'lh', 'rh'};
for map_idx = 1:length(maps)
    template = fullfile(vistaRootPath, 'fileFilters', 'freesurfer', ...
        sprintf('mh.V1.%stmp.sym.mgh', maps{map_idx}));
    for hemi_idx = 1:length(hemis)
        file_root = fullfile(out_path, sprintf('%s_%s_%s', subject, ...
                             hemis{hemi_idx} , maps{map_idx}));
        % These commands differ slightly for the two hemispheres (xhemi):
        if strcmp(hemis{hemi_idx}, 'lh')
            % If the registration file is not there, we'll have to make it:
            if ~(exist(fullfile(subjects_dir, subject, ...
                    '/surf/lh.fsaverage_sym.sphere.reg'), 'file') == 2)
                fprintf('[%s]: Registering... This might take a while... \n',...
                    mfilename);

                cmd_str1 = sprintf('surfreg --s %s --t fsaverage_sym --lh', subject);
            else
                cmd_str1 = '';
            end
            cmd_str2 = ['mri_surf2surf --srcsubject fsaverage_sym --trgsubject ' ...
                sprintf('%s --sval %s --tval ', subject, template) ...
                sprintf('%s.mgh --hemi lh', file_root)];
        else
            if ~(exist(fullfile(subjects_dir, subject, ...
                    'xhemi/surf/lh.fsaverage_sym.sphere.reg'), 'file') == 2)
                fprintf('[%s]: Registering... This might take a while... \n',...
                    mfilename);

                % Yep - the following line is with 'lh':
                cmd_str1 = sprintf('surfreg --s %s --t fsaverage_sym --lh --xhemi', subject);
            else
                cmd_str1 = '';
            end
            cmd_str2 = ['mri_surf2surf --srcsubject fsaverage_sym --trgsubject ' ...
                sprintf('%s/xhemi --sval %s --tval ', subject, template) ...
                sprintf('%s.mgh --srcsurfreg sphere.reg --trgsurfreg ', file_root, subject) ...
                'fsaverage_sym.sphere.reg --hemi lh'];
        end
        syscall(cmd_str1);
        syscall(cmd_str2);

        cmd_str = [sprintf('mri_surf2vol --surfval %s.mgh --projfrac 1 ', file_root) ...
            sprintf('--identity %s --o %s.mgz --hemi %s ', subject, file_root, hemis{hemi_idx}) ...
            sprintf('--template %s/%s/mri/orig.mgz', subjects_dir, subject)];
        syscall(cmd_str);

        % Convert to volume:
        outfile = [file_root '.nii'];
        cmd_str = [sprintf('mri_convert  %s.mgz ', file_root), outfile];
        syscall(cmd_str);
        gzip(outfile);
        delete(outfile);

    end
end
end

% Helper function: throw an error if the system call doesn't work as
% expected:
    function [status, result] = syscall(cmd_str)
        % Allow for noops:
        if strcmp(cmd_str, '')
            return
        end
        fprintf('[%s]: Executing "%s" \n', mfilename, cmd_str);
        [status, result] = system(cmd_str);
        if status~=0
            error(result);
        end
    end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% VBRFA2010_DEMO - A simple comparison of VB PCA and two robust extensions.
%
% The algorithms do not use ARD prior for the loadings W in order to keep
% the visualizations of the subspaces more clear. In addition, it helps
% avoiding pruning out relevant components too easily.
%
% NOTE: The robust algorithms might be less sensitive to the initialization
% if they estimate many components. With only one estimated component, they
% might kill all the components.
%
% This demo doesn't always show good results for robust algorithms. You may
% want to run a few times. In addition, sometimes there are several good
% interpretations of the data and the robust algorithms pick one.

% Last modified 2010-06-10
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function Q = vbrfa2010_demo

n = 50;
m = 2;
d = 1;
dh = d + 0; % number of components to estimate

randn('state', 6);
rand('state', 6);

mu = [0;0];
W = [2; 1];
X = orth(randn(d,n)')' * sqrt(n);
s2 = 0.2;

% Generate normal data
Y = W*X + repmat(mu,1,n);

% Add noise
Yn = Y + sqrt(s2)*randn(m,n);

% Generate some outliers
Yno = Yn;
p = (rand(m,n) < 0.05);
Yno(p) = 5 + 8 * rand(sum(p(:)),1);

% Generate some missing values
Ynom = Yno;
pmv = (rand(m,n) < 0.0);
Ynom(pmv) = NaN;

% Run different algorithms
options.init.tau = 1e2*ones(m,1); % the initialization of this can be crucial..
options.init.nu = 1;
options.prior.aw = 1e10*ones(d,1); % fix w by setting strong prior
options.prior.bw = 1e15*ones(d,1); % fix w by setting strong prior
options.update_nu = 1;
options.update_w = 1;
options.rotate = 0;
options.common_nu = false;
options.common_tau = true
options.maxiter = 50;

% PCA
disp('Run PCA')
options.robustness = 'none';
results_p = vbrfa(Ynom, dh, options);
W_p = results_p.W;
X_p = results_p.X;
S_p = results_p.Sv;
mu_p = results_p.mu;
nu_p = results_p.nu;

% Robust PCA with multivariate Student-t
disp('Run robust PCA with multivariate Student-t')
options.robustness = 'multivariate-t';
results_rp = vbrfa(Ynom, dh, options);
W_rp = results_rp.W;
X_rp = results_rp.X;
S_rp = results_rp.Sv;
mu_rp = results_rp.mu;
nu_rp = results_rp.nu;

% Robust PCA with independent Student-t
disp('Run robust PCA with independent Student-t')
options.robustness = 'independent-t';
results_rvb = vbrfa(Ynom, dh, options);
W_rvb = results_rvb.W;
X_rvb = results_rvb.X;
S_rvb = results_rvb.Sv;
mu_rvb = results_rvb.mu;
nu_rvb = results_rvb.nu;

% Reconstruct
Y_p = W_p*X_p + repmat(mu_p,1,n);
Y_rp = W_rp*X_rp + repmat(mu_rp,1,n);
Y_rvb = W_rvb*X_rvb + repmat(mu_rvb,1,n);

plot_results(Ynom, []);
%print('-depsc2', '/home/jluttine/papers/vbrfa/figures_journal/fig_demo_data.eps');
plot_results(Ynom, Y_p)
%print('-depsc2', '/home/jluttine/papers/vbrfa/figures_journal/fig_demo_pca.eps');
plot_results(Ynom, Y_rp)
%print('-depsc2', '/home/jluttine/papers/vbrfa/figures_journal/fig_demo_multi-t.eps');
plot_results(Ynom, Y_rvb)
%print('-depsc2', '/home/jluttine/papers/vbrfa/figures_journal/fig_demo_indep-t.eps');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
function plot_results(Y, Yh)
%subspace2d(Yh, [], Y)
figure
plot(Y(1,:), Y(2,:), 'x', 'markersize', 5, 'color', [.2 .2 .2]);
if ~isempty(Yh)
  hold on
  plot(Yh(1,:), Yh(2,:), 'ko', 'markersize', 3);
  nans = nan*ones(1,size(Y,2));
  plot([Y(1,:);Yh(1,:);nans], [Y(2,:);Yh(2,:);nans], ':', 'color', [.0 .0 .0]);
end
set(gca, 'DataAspectRatioMode', 'manual');
set(gca, 'DataAspectRatio', [1 1 1]);
set(gca, 'XTick', [], 'YTick', []);
set(gcf, 'units', 'centimeters');
pos = get(gcf, 'position');
set(gcf, 'position', [pos(1), pos(2), 5 5])
set(gcf, 'PaperPositionMode', 'auto', 'paperunits', 'centimeters', 'PaperSize', [5 5]);

xmg = 0.05 * (max(Y(1,:)) - min(Y(1,:)));
ymg = 0.05 * (max(Y(2,:)) - min(Y(2,:)));
xl = [min(Y(1,:))-xmg, max(Y(1,:))+xmg];
yl = [min(Y(2,:))-ymg, max(Y(2,:))+ymg];
set(gca, 'xlim', xl, 'ylim', yl);

% $$$ % Mark outliers
% $$$ hold on
% $$$ indeces = sum(p,1)>0;
% $$$ plot(VYno(1,indeces), VYno(2, indeces), 'go', 'MarkerSize', 8);

% $$$ % Mark observations with missing values
% $$$ indeces = sum(pmv,1)>0;
% $$$ scatter(VYno(1,indeces), VYno(2, indeces), 'yo');

return

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function frame = FrameFromEvents(polarity, contrast, dims)

%{
Polarity is a structure containing at least the fields x, y, and polarity,
column vectors of equal lengths where polarity is boolean
and x and y are zero-based coordinates.
returns 'frame'(uint8) a matrix of plus/minus plots of polarity data a la jAER.
Encoding is from 0 (full scale OFF events)
to contrast * 2 (full scale ON events, with contrast representing the
neutral value.
dims is the expected dimension of the output - a tuple of [X Y]

Note that the frame returned accords to the maximum size of the x and y
events.
%}

% Output class is double, because the numerical operation on polarity
% results in a double
if exist('dims', 'var')
    frame = accumarray([polarity.y polarity.x] + 1, polarity.polarity * 2 - 1, dims(end : -1 : 1));
else
    frame = accumarray([polarity.y polarity.x] + 1, polarity.polarity * 2 - 1);
end

% Clip the values according to the contrast
frame(frame > contrast) = contrast;
frame(frame < - contrast) = -contrast;
frame = uint8(frame + contrast + 1);

author: panpanfei


text:
function transfo = createBasisTransform(source, target)
%CREATEBASISTRANSFORM Compute matrix for transforming a basis into another basis.
%
%   TRANSFO = createBasisTransform(SOURCE, TARGET)
%   Both SOURCE and TARGET represent basis, in the following form:
%   [x0 y0  ex1 ey1  ex2 ey2]
%   [y0 y0] is the origin of the basis, [ex1 ey1] is the first direction
%   vector, and [ex2 ey2] is the second direction vector.
%
%   The result TRANSFO is a 3-by-3 matrix such that a point expressed with
%   coordinates of the first basis will be represented by new coordinates
%   P2 = transformPoint(P1, TRANSFO) in the target basis.
%
%   TRANSFO = createBasisTransform(TARGET)
%   Assumes the source is the standard (Oij) basis, with origin at (0,0),
%   first direction vector equal to (1,0) and second direction  vector
%   equal to (0,1).
%
%
%   Example
%     % define source and target bases
%     src = [ 0 0   1  0    0  1];
%     tgt = [20 0  .5 .5  -.5 .5];
%     trans = createBasisTransform(src, tgt);
%     % create a polygon in source basis
%     poly = [10 10;30 10; 30 20; 20 20;20 40; 10 40];
%     figure;
%     subplot(121); drawPolygon(poly, 'b'); axis equal; axis([-10 50 -10 50]);
%     hold on; drawLine([0 0 1 0], 'k'); drawLine([0 0 0 1], 'k');
%     drawLine([20 0 1 1], 'r'); drawLine([20 0 -1 1], 'r');
%     t = -1:5; plot(t*5+20, t*5, 'r.'); plot(-t*5+20, t*5, 'r.');
%     % transform the polygon in target basis
%     poly2 = transformPoint(poly, trans);
%     subplot(122); drawPolygon(poly2, 'b'); axis equal; axis([-10 50 -10 50]);
%     hold on; drawLine([0 0 1 0], 'r'); drawLine([0 0 0 1], 'r');
%     t = -1:5; plot(t*10, zeros(size(t)), 'r.'); plot(zeros(size(t)), t*10, 'r.');
%
%   See also
%   transforms2d
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2010-12-03, using Matlab 7.9.0.529 (R2009b)
% Copyright 2010-2022 INRA - Cepia Software Platform

% init basis transform to identity
t1 = eye(3);
t2 = eye(3);

if nargin == 2
    % from source to reference basis
    t1(1:2, 1) = source(3:4);
    t1(1:2, 2) = source(5:6);
    t1(1:2, 3) = source(1:2);
else
    % if only one input, use first input as target basis, and leave the
    % first matrix to identity
    target = source;
end

% from reference to target basis
t2(1:2, 1) = target(3:4);
t2(1:2, 2) = target(5:6);
t2(1:2, 3) = target(1:2);

% compute transform matrix
transfo = zeros(3, 3);
maxSz = 1;
for i = 1:maxSz
    % coordinate of three reference points in source basis
    po = t1(1:2, 3, i)';
    px = po + t1(1:2, 1, i)';
    py = po + t1(1:2, 2, i)';

    % express coordinates of reference points in the new basis
    t2i = inv(t2(:,:,i));
    pot = transformPoint(po, t2i);
    pxt = transformPoint(px, t2i);
    pyt = transformPoint(py, t2i);

    % compute direction vectors in new basis
    vx = pxt - pot;
    vy = pyt - pot;

    % concatenate result in a 3-by-3 affine transform matrix
    transfo(:,:,i) = [vx' vy' pot' ; 0 0 1];
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function plotranking(scores, X, dims)
% PLOTRANKING(DAT)
% PLOTRANKING(SCORES, X)
%
% DAT is a data object returned by the MRANK method TEST.
%
% SCORES may be one column [scores], two columns [indices, scores] or three
% columns [indices, scores, ranks] (the latter is the format returned by
% READRANKING or in the result field of an MRANK algorithm).
%
% SCORES may also be a string, in which case it is automatically passed
% through READRANKING (to plot the contents of an output file).
%
% Use NaN in the scores column to mark labelled points.
%
% X may be a string: if so, a variable X is loaded from the named file.
% If X is omitted, 'twomoons' is assumed. The first two columns of X
% provide the coordinates for the plotted data points.

if nargin < 3, dims = [1 2]; end
if nargin < 2, X = ''; end
if nargin < 1, scores = ''; end

if nargin == 1 & (isa(scores, 'data') | isa(scores, 'data_global'))
	X = get_x(scores);
	scores = get_y(scores);
end

if isempty(X), X = 'twomoons'; end
if isstr(X), load(X, 'X'); end
if isa(X, 'data') | isa(X, 'data_global'), X = get_x(X); end
X = X(:, dims);

if isstr(scores), scores = readranking(scores); end

m = size(X, 1);

if isempty(scores), scores = [1:m]'; end

if size(scores, 2) == 1
	if size(scores, 1) ~= m, error('SCORES must have one element per data point unless indices are supplied in the first column'), end
	scores = [[1:m]' scores];
end
if size(scores, 2) == 2
	labelled = isnan(scores(:, 2));
	scores(labelled, 3) = nan;
	t = compute_ranks(-scores(~labelled, 2));
	scores(~labelled, 3) = t(:);
end

X = X(scores(:, 1), :);
ranks = scores(:, 3);
scores = scores(:, 2);

labelled = isnan(scores);

plotpoints(X(~labelled, :), ranks(~labelled))
plot(X(labelled, 1), X(labelled, 2), 'marker', 'o', 'color', [1 0 0], 'markerfacecolor', [0 1 0], 'markersize', 10, 'linestyle', 'none')
axis square
figure(gcf)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plotpoints(X, ranks)

fullcolour = 0;
m = size(X, 1);
cmap = jet(m);
[ranks i] = sort(max(ranks)-ranks);
X = X(i, :);
s = 10 + 90*ranks/m;
c = [1:m]';
if ~fullcolour, c = round(m/4)+round(m/2)*(c<m/2); end
c = cmap(c, :);
washeld = ishold;
plot(X(:, 1), X(:, 2), 'color', [1, 0,0])
hold on
if fullcolour
	scatter(X(:, 1), X(:, 2), 20, c, 'filled')
else
	scatter(X(:, 1), X(:, 2), s, c)
end
if ~washeld, hold off, end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rm = compute_ranks(am)

rm = zeros(size(am));
for j = 1:size(am, 2)
    a = am(:, j);
    [sorted r] = sort(a);
    r(r) = 1:length(r);
    while 1
        tied = sorted(min(find(diff(sorted) == 0)));
        if isempty(tied), break, end
        sorted(sorted==tied) = nan;
        r(a==tied) = mean(r(a==tied));
    end
    rm(:, j) = r;
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function f = filterBinomial1d( r, show )
% 1D binomial filter (approximation to Gaussian filter)
%
% Creates a binomial filter of size 2*r+1 x 1.  This can be used to
% approximate the Gaussian distribution with sigma=sqrt((2*r+1)/4). For
% large r, should give same output as:
%  g = fspecial( 'Gaussian', [2*r+1,1],sqrt((2*r+1)/4) );
% Given sigma, use r ~= 2*sigma^2.
%
% Use F = f*f' to get the equivalent 2d filter.
%
% USAGE
%  f = filterBinomial1d( r, [show] )
%
% INPUTS
%  r       - mask will have length 2r+1 and var=(2*r+1)/4
%  show    - [0] figure to use for optional display
%
% OUTPUTS
%  f       - 1D binomial filter
%
% EXAMPLE
%  r = 10;
%  fbinom = filterBinomial1d( r, 1 );
%  fgauss = filterGauss( 2*r+1, [], (2*r+1)/4, 2);
%
% See also FILTERGAUSS
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<2 ); show=0; end;
if( mod(r,1)~=0 ); error( 'r must be an integer'); end;

f = diag(fliplr(pascal(2*r+1))) / 4^r;

% display
if(show); filterVisualize( f, show ); end;

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function e = ExplorationNoise(t)
%% LocalExploration
% Generate exploratoion noise
  e = (0.01*sin(10*t)+0.01*sin(3*t)+0.01*sin(100*t));
end

author: yu-jiang

repo: radpbook

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


text:
function Show_EWT2D_Curvelet(ewtc)

%==========================================================================
% function Show_EWT2D_Curvelet(ewtc)
%
% This function displays the curvelet coefficient obtained by the Empirical
% Curvelet Transform (one figure per scale)
%
% Inputs:
%   -ewtc: output of the empirical curvelet transform
%
% Author: Jerome Gilles
% Institution: UCLA - Department of Mathematics
% Year: 2013
% Version: 1.0
%==========================================================================

Ns=length(ewtc);

% Show the low pass filter
figure;
imshow(ewtc{1},[]);

for s=2:Ns
   Nt=length(ewtc{s});
   figure;
   for t=1:Nt
      subplot(1,Nt,t);
      imshow(ewtc{s}{t},[]);
   end
end


text:
function mesh = spm_eeg_inv_spatnorm(mesh)
% Spatial Normalisation (using Unified Segmentation)
% Transforms individual sMRI into MNI space and save the [inverse]
% deformations that will be needed for computing the individual mesh
%
% FORMAT mesh = spm_eeg_inv_spatnorm(mesh)
%
% mesh        - input data struct
%
% mesh        - same data struct including the inverse deformation .mat file
%               and filename of normalised (bias corrected) sMRI
%__________________________________________________________________________
% Copyright (C) 2009-2014 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_eeg_inv_spatnorm.m 6217 2014-09-29 17:54:33Z guillaume $

spm('Pointer','Watch');

%-Initialise
%--------------------------------------------------------------------------
sMRI = mesh.sMRI;

%-Spatial Transformation into MNI space
%--------------------------------------------------------------------------
[pth,nam] = spm_fileparts(sMRI);
def       = fullfile(pth,['y_' nam '.nii']);
mat       = fullfile(pth,[nam '_seg8.mat']);

if ~(exist(def, 'file') && exist(mat, 'file'))
    % Assume it is an image, so derive deformation field.
    T = fullfile(spm('dir'),'tpm','TPM.nii,');
    p = struct('channel',struct('vols',     {{sMRI}},...
        'biasreg',  0.0001,...
        'biasfwhm', 60,...
        'write',    [0 0]),...
        'tissue', struct('tpm',   {{[T,'1']},{[T,'2']},{[T,'3']},{[T,'4']},{[T,'5']},{[T,'6']}},...
        'ngaus', {2,2,2,3,4,2},...
        'native',{[0 0],[0 0],[0 0],[0 0],[0 0],[0 0]},...
        'warped',{[0 0],[0 0],[0 0],[0 0],[0 0],[0 0]}),...
        'warp',   struct('reg',[0 0.001 0.5 0.05 0.2], 'affreg', 'mni', 'samp', 3, 'write', [0 1],'mrf',0));

    spm_preproc_run(p);
end

mesh.def    = def;
mesh.Affine = getfield(load(mat, 'Affine'), 'Affine');

%-Clean up
%--------------------------------------------------------------------------
spm('Pointer','Arrow');

author: spm


text:
function res = nsamples(obj)
% Method for getting the number of samples per trial
% FORMAT res = nsamples(obj)
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: nsamples.m 1373 2008-04-11 14:24:03Z spm $

res = obj.Nsamples;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function alpha = angle3d(varargin)
%ANGLE3D compute angle between 2 3D points
%
%   ALPHA = anglePoints3d(P1, P2)
%   compute angle (P1, O, P2), in radians, between 0 and PI.
%
%   ALPHA = anglePoints3d(P1, P2, P3)
%   compute angle (P1, P2, P3), in radians, between 0 and PI.
%
%   ALPHA = anglePoints3d(PTS)
%   PTS is a 3x3 or 2x3 array containing coordinate of points.
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 21/02/2005.
%

%   HISTORY
%   20/09/2005 : add case of single argument for all points

p2 = [0 0 0];
if length(varargin)==1
    pts = varargin{1};
    if size(pts, 1)==2
        p1 = pts(1,:);
        p0 = [0 0 0];
        p2 = pts(2,:);
    else
        p1 = pts(1,:);
        p0 = pts(2,:);
        p2 = pts(3,:);
    end
elseif length(varargin)==2
    p1 = varargin{1};
    p0 = [0 0 0];
    p2 = varargin{2};
elseif length(varargin)==3
    p1 = varargin{1};
    p0 = varargin{2};
    p2 = varargin{3};
end

% normalized points
p1 = normalize3d(p1-p0);
p2 = normalize3d(p2-p0);
alpha = acos(dot(p1, p2, 2));

function vn = normalize3d(v)
%NORMALIZE3D normalize a 3D vector, such that its norm equals 1.

n = sqrt(v(:,1).*v(:,1) + v(:,2).*v(:,2) + v(:,3).*v(:,3));
vn = v./[n n n];

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
% Change to your downloaded location
clear
addpath('C:\liblinear\matlab')
addpath('../training_code')
addpath('../utilities')

%% load shared definitions and AU data
shared_defs;

% Set up the hyperparameters to be validated
hyperparams.c = 10.^(-7:1:1);
hyperparams.e = 10.^(-3);

hyperparams.validate_params = {'c', 'e'};

% Set the training function
svr_train = @svm_train_linear;

% Set the test function (the first output will be used for validation)
svr_test = @svm_test_linear;

%%
for a=1:numel(aus)

    au = aus(a);

    rest_aus = setdiff(all_aus, au);

    % make sure validation data's labels are balanced
    [users_train, users_valid] = get_balanced_fold(DISFA_dir, users, au, 1/3, 1);

    % need to split the rest
    [train_samples, train_labels, valid_samples, valid_labels, ~, PC, means, scaling, valid_ids, valid_success] = Prepare_HOG_AU_data_generic(users_train, users_valid, au, rest_aus, DISFA_dir, hog_data_dir);

    train_labels(train_labels > 1) = 1;
    valid_labels(valid_labels > 1) = 1;

    train_samples = sparse(train_samples);
    valid_samples = sparse(valid_samples);

    %% Validate here
    hyperparams.success = valid_success;

    [ best_params, ~ ] = validate_grid_search_no_par(svr_train, svr_test, false, train_samples, train_labels, valid_samples, valid_labels, hyperparams);

    model = svr_train(train_labels, train_samples, best_params);

    [~, prediction] = svr_test(valid_labels, valid_samples, model);

    name = sprintf('classifiers/AU_%d_stat.mat', au);

    [ accuracies, F1s, corrs, ccc, rms, classes ] = evaluate_regression_results( prediction, valid_labels );

    save(name, 'model', 'accuracies', 'F1s', 'corrs', 'rms', 'ccc', 'prediction', 'valid_labels');

    name = sprintf('classifiers/AU_%d_stat.dat', au);

    pos_lbl = model.Label(1);
    neg_lbl = model.Label(2);

    w = model.w(1:end-1)';
    b = model.w(end);

    svs = bsxfun(@times, PC, 1./scaling') * w;
    write_lin_svm(name, means, svs, b, pos_lbl, neg_lbl);

end

author: TadasBaltrusaitis

repo: OpenFace


text:
% PROBABILITY DISTRIBUTION FUNCTIONS (in the dist-folder):
%
% Priors
%   PRIOR_FIXED       Fix parameter to its current value
%   PRIOR_GAMMA       Gamma prior structure
%   PRIOR_GAUSSIAN    Gaussian prior structure
%   PRIOR_INVGAMMA    Inverse-gamma prior structure
%   PRIOR_INVT        Inverse Student-t prior structure
%   PRIOR_INVUNIF     Inverse uniform prior structure
%   PRIOR_LAPLACE     Laplace (double exponential) prior structure
%   PRIOR_LOGGAUSSIAN Log-Gaussian prior structure
%   PRIOR_LOGLOGUNIF  Uniform prior structure for the log(log(parameter))
%   PRIOR_LOGT        Student-t prior structure for the logarithm of the parameter
%   PRIOR_LOGUNIF     Uniform prior structure for the logarithm of the parameter
%   PRIOR_SINVCHI2    Scaled inverse-chi-square prior structure
%   PRIOR_SQINVGAMMA  Gamma prior structure for square inverse of the parameter
%   PRIOR_SQINVLOGUNIF Uniform prior structure for the log of the square
%                      inverse of parameter
%   PRIOR_SQINVSINVCHI2 Scaled-Inv-Chi^2 prior structure for square inverse
%                       of the parameter
%   PRIOR_SQINVUNIF   Uniform prior structure for the square inverse of the
%                     parameter
%   PRIOR_SQRTINVT    Student-t prior structure for the square root of
%                     inverse of the parameter
%   PRIOR_INVSQRTUNIF Uniform prior structure for the square root of
%                     inverse of the parameter
%   PRIOR_SQRTT       Student-t prior structure for the square root of the
%                     parameter
%   PRIOR_SQRTUNIF    Uniform prior structure for the square root of the
%                     parameter
%   PRIOR_T           Student-t prior structure
%   PRIOR_UNIF        Uniform prior structure
%
% Probability (log/cumulative) density functions
%   BETA_CDF      - Beta cumulative distribution function
%   BETA_INV      - Inverse of the beta cumulative distribution function (cdf)
%   BETA_LPDF     - Beta log-probability density function (lpdf)
%   BETA_PDF      - Beta probability density function (pdf)
%   DIR_LPDF      - Log probability density function of uniform Dirichlet
%                   distribution
%   DIR_PDF       - Probability density function of uniform Dirichlet
%                   distribution
%   GAM_CDF       - Cumulative of Gamma probability density function (cdf)
%   GAM_LPDF      - Log of Gamma probability density function (lpdf)
%   GAM_PDF       - Gamma probability density function (pdf)
%   GEO_LPDF      - Geometric log probability density function (lpdf)
%   INVGAM_LPDF   - Inverse-Gamma log probability density function
%   INVGAM_PDF    - Inverse-Gamma probability density function
%   LAPLACE_LPDF  - Laplace log-probability density function (lpdf)
%   LAPLACE_PDF   - Laplace probability density function (pdf)
%   LOGN_LPDF     - Log normal log-probability density function (lpdf)
%   LOGT_LPDF     - Log probability density function (lpdf) for log Student's T
%   MNORM_LPDF    - Multivariate-Normal log-probability density function (lpdf)
%   MNORM_PDF     - Multivariate-Normal log-probability density function (lpdf)
%   NBIN_CDF      - Negative binomial cumulative distribution function (cdf)
%   NBIN_INV      - Inverse of Negative binomial cumulative distribution function (inv)
%   NBIN_PDF      - Negative binomial probability density function (pdf)
%   NEGBIN_LPDF   - Negative binomial log probability density function
%   NEGBIN_PDF    - Negative binomial probability density function
%   NEGBINZTR_LPDF - Zero trunc. negative binomial log probability density function
%   NEGBINZTR_PDF - Zero trunc. negative binomial log probability density function
%   NORM_CDF      - Normal cumulative probability density function (cdf)
%   NORM_INV      - Inverse of the normal cumulative distribution function (cdf)
%   NORM_LPDF     - Normal log-probability density function (lpdf)
%   NORM_PDF      - Normal probability density function (pdf)
%   POISS_LPDF    - Poisson log-probability density function
%   POISS_PDF     - Poisson probability density function
%   SINVCHI2_LPDF - Scaled inverse-chi log-probability density function
%   SINVCHI2_PDF  - Scaled inverse-chi probability density function
%   T_CDF         - Student's t cumulative distribution function (cdf)
%   T_INV         - Inverse of Student's T cumulative distribution function (cdf)
%   T_LPDF        - Student's T log-probability density function (lpdf)
%   T_PDF         - Student's T probability density function (pdf)
%
% Random numbers
%   BETARAND      - Random matrices from beta distribution
%   CATRAND       - Random matrices from categorical distribution
%   DIRRAND       - Uniform dirichlet random vectors
%   EXPRAND       - Random matrices from exponential distribution
%   GAMRAND       - Random matrices from gamma distribution
%   GAMRAND1      - Random matrices from gamma distribution (mex)
%   INTRAND       - Random matrices from uniform integer distribution
%   INVGAMRAND    - Random matrices from inverse gamma distribution
%   INVGAMRAND1   - Random matrices from inverse gamma distribution (mex)
%   INVWISHRND    - Random matrices from inverse Wishart distribution
%   NORMLTRAND    - Random draws from a left-truncated normal
%                   distribution, with mean = mu, variance = sigma2
%   NORMRTRAND    - Random draws from a right-truncated normal
%                   distribution, with mean = mu, variance = sigma2
%   NORMTRAND     - Random draws from a normal truncated to interval
%   NORMTZRAND    - Random draws from a normal distribution truncated by zero
%   SINVCHI2RAND  - Random matrices from scaled inverse-chi distribution
%   TRAND         - Random numbers from Student's t-distribution
%   UNIFRAND      - Generate unifrom random numberm from interval [A,B]
%   WISHRND       - Random matrices from Wishart distribution
%
% Others
%   HAMMERSLEY    - Hammersley quasi-random sequence
%

author: gpstuff-dev

repo: gpstuff

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


text:
function create_tsne_Callback(hObject, eventdata, handles)
% Creates a 2-dimensional t-sne image of calls.

padding = 1000; % Pad the temp image by this amount, so that calls near the border still fit.
blackLevel = 10; % Subtract this value from each call image to make a nicer picture.

% Select embedding type
embeddingType = questdlg('Embed with UMAP or t-SNE?', 'Embedding Method', 't-SNE' , 'UMAP', 't-SNE');
if isempty(embeddingType); return; end
if strcmp(embeddingType, 'UMAP')
    if ~exist('run_umap.m', 'file')
        msgbox('Please download UMAP and add it to MATLAB''s path and try again')
        web('   https://www.mathworks.com/matlabcentral/fileexchange/71902-uniform-manifold-approximation-and-projection-umap');
        return
    end
end

%% Choose clustering dimensions
inputParameters= questdlg('Select the input dimensions', 'Input Dimensions', 'Variational autoencoder embeddings', 'Contour shape, frequency, and duration', 'Contour shape, frequency, and duration');
if isempty(inputParameters); return; end

% Choose to assign colors by the call classification or by pitch.
colorType = questdlg({'Color the calls by frequecy (pitch), or by cluster identity?','If coloring by cluster, you may not use pre-extracted contours'}, 'Choose Color', 'Frequency' , 'Cluster', 'Frequency');
if isempty(colorType); return; end

switch inputParameters
    case 'Contour shape, frequency, and duration'
        % Get the clustering parameters, prepare data as if performing k-means
        clusterParameters = inputdlg({'Images height (pixels)','Image width (pixels)','Perplexity (if using t-SNE, ignore for UMAP)','Max number of calls to plot (set to 0 to plot everything)', 'Shape weight','Frequency weight','Duration weight'},...
            'Choose cluster parameters:',1,...
            {'6000','6000','30','2000','1','1','1'});
        if isempty(clusterParameters); return; end

        slope_weight = str2double(clusterParameters{5});
        freq_weight = str2double(clusterParameters{6});
        duration_weight = str2double(clusterParameters{7});

        % Get the data
        [ClusteringData, clustAssign] = CreateClusteringData(handles, 'forClustering', true);

        %% Extract features
        ReshapedX   = cell2mat(cellfun(@(x) imresize(x',[1 9]) ,ClusteringData.xFreq,'UniformOutput',0));
        slope       = diff(ReshapedX,1,2);
        slope       = zscore(slope);
        freq        = cell2mat(cellfun(@(x) imresize(x',[1 8]) ,ClusteringData.xFreq,'UniformOutput',0));
        freq        = zscore(freq);
        duration    = repmat(ClusteringData.Duration,[1 8]);
        duration    = zscore(duration);

        data = [
            freq     .*  freq_weight,...
            slope    .*  slope_weight,...
            duration .*  duration_weight,...
            ];

    case 'Variational autoencoder embeddings'
        clusterParameters = inputdlg({'Images height (pixels)','Image width (pixels)','Perplexity (if using t-SNE, ignore for UMAP)','Max number of calls to plot (set to 0 to plot everything)'},...
            'Choose cluster parameters:',1,...
            {'6000','6000','30','2000'});
        if isempty(clusterParameters); return; end

        FromExisting = questdlg('From existing model?','Cluster','Yes','No','No');
        if isempty(FromExisting); return; end
        switch FromExisting % Load Model
            case 'No'
                [encoderNet, decoderNet, options, ClusteringData] = create_VAE_model(handles);
                [FileName, PathName] = uiputfile(fullfile(handles.data.squeakfolder, 'Clustering Models', 'Variational Autoencoder Model.mat'), 'Save clustering model (optional)');
                if ~isnumeric(FileName) % Save the new model
                    save(fullfile(PathName, FileName), 'encoderNet', 'decoderNet', 'options');
                end
            case 'Yes' % Load the VAE model
                [FileName, PathName] = uigetfile(fullfile(handles.data.squeakfolder, 'Clustering Models', '*.mat'), 'Select VAE model');
                if isnumeric(FileName);return;end
                load(fullfile(PathName,FileName),'encoderNet','options');
                [ClusteringData, clustAssign] = CreateClusteringData(handles, 'spectrogramOptions', options.spectrogram, 'scale_duration', options.maxDuration, 'freqRange', options.freqRange, 'forClustering', true);
                if isempty(ClusteringData); return; end
        end
        data = extract_VAE_embeddings(encoderNet, options, ClusteringData);

end

imsize = str2double(clusterParameters(1:2))';
perplexity = str2double(clusterParameters{3});
NumberOfCalls = str2double(clusterParameters{4});

% Plot all calls if number of calls is set to 0
if NumberOfCalls == 0
    NumberOfCalls = size(ClusteringData,1);
end
NumberOfCalls = min(size(ClusteringData,1), NumberOfCalls);

%% Get parameters

% Calculate embeddings
rng default;

% Run embedding
switch embeddingType
    case 't-SNE'
        embed = tsne(data,'Verbose',1,'Perplexity',perplexity);
    case 'UMAP'
        embed = run_umap(data);
end
% Rescale values between 0 and 1
embed = (embed - min(embed)) ./ (max(embed)-min(embed));

ClusteringData.embedY = 1-embed(:,2); % flip Y coordinates so the images looks like the UMAP figure
ClusteringData.embedX = embed(:,1);

switch colorType
    case 'Frequency'
        minfreq = prctile(ClusteringData.MinFreq, 1);
        maxfreq = prctile(ClusteringData.MinFreq + ClusteringData.Bandwidth, 99);
        ColorData = HSLuv_to_RGB(256, 'H',  [0 270], 'S', 100, 'L', 75); % Make a color map for each category
        ColorData = reshape(ColorData,size(ColorData,1),1,size(ColorData,2));
    case 'Cluster'
        [clustAssignID, cName] = findgroups(clustAssign); % Convert categories into numbers
        ClusteringData.Cluster = clustAssignID; % Append the category number to clustering data

        % make it so that adjacent clusters are generally different colors.
        % it turns out that this isn't trivial, so try 200 different color
        % orders and use the best one.
        embedings=[ClusteringData.embedY,  ClusteringData.embedX];
        z=unique(clustAssignID);
        for c=1:height(z)
            idx=clustAssignID==z(c);
            clusterCentroids(c,:)=mean(embedings(idx,:));
        end
        % = splitapply(@mean,[ClusteringData.embedY,  ClusteringData.embedX], clustAssignID);
        clusterCentroids = pdist2(clusterCentroids, clusterCentroids);
        hueAngle = [
            sin(linspace(0,2*pi, length(cName)))
            cos(linspace(0,2*pi, length(cName)))]';
        hueDistance = pdist2(hueAngle, hueAngle, 'cosine');
        colorSeperation = zeros(200,1);
        colorOrders = zeros(200,length(cName));
        for i = 1:200
            colorOrders(i,:) = randperm(length(cName)); % Randomize the order of the colors
            colorSeperation(i) = sum((1-clusterCentroids).* hueDistance(colorOrders(i,:),colorOrders(i,:)), 'all');
        end
        [~, idx] = max(colorSeperation)
        colorOrder = colorOrders(idx,:)

        hueRange = [0 360-360/length(cName)];
        cMap = HSLuv_to_RGB(length(cName), 'H', hueRange, 'S', 100, 'L', 75); % Make a color map for each category
        cMap = cMap(colorOrder,:);
        figure('Color','w') % Display the colors
        h = image(reshape(cMap,[],1,3));
        yticklabels(h.Parent, cellstr(cName));
        yticks(h.Parent,1:length(cName));
end

%% Create the image
im = zeros([imsize+padding*2,3],'uint8');

% Only plot the X number of calls
calls2plot = randsample(size(ClusteringData, 1), NumberOfCalls, false);

for i = calls2plot'
    call = ClusteringData(i, :);
    % Get x and y coordinates to place with image
    iy = imsize(1) * call.embedY;
    iy = iy:iy + size(call.Spectrogram{:},1) - 1;
    iy = round(iy - mean(size(call.Spectrogram{:},1))) + padding;
    ix = imsize(2) * call.embedX;
    ix = ix:ix + size(call.Spectrogram{:},2) - 1;
    ix = round(ix - mean(size(call.Spectrogram{:},2))) + padding;

    % Either use the call pitch or the cluster id to apply a color mask
    switch colorType
        case 'Frequency'
            % Interpolate the color data
            freqdata = linspace(call.MinFreq, call.MinFreq + call.Bandwidth, size(call.Spectrogram{:},1));
            colorMask = interp1(linspace(minfreq, maxfreq, size(ColorData,1)), ColorData, freqdata, 'nearest', 'extrap');
        case 'Cluster'
            colorMask = reshape(cMap(call.Cluster,:),1,1,3);
    end
    im(iy,ix,:) = max(im(iy,ix,:),uint8(single(call.Spectrogram{:}) .* colorMask - blackLevel));
end

% Crop the image at the first and last non-empty pixel
[y1,~] = find(max(max(im,[],3),[],2),1,'first');
[y2,~] = find(max(max(im,[],3),[],2),1,'last');

[~,x1] = find(max(max(im,[],3),[],1),1,'first');
[~,x2] = find(max(max(im,[],3),[],1),1,'last');

[fname,fpath] = uiputfile({'*.jpg'; '*.png'},'Save image', 'embeddings');
imwrite(im2uint8(im(y1:y2,x1:x2,:)),fullfile(fpath,fname))

% Open the image in a file manager
if ispc % Open the file in windows explorer
    system(['explorer.exe /select,"' fullfile(fpath,fname) '"']);
    % winopen(fullfile(fpath,fname));
elseif ismac % Open the file in finder - UNTESTED (I don't have a mac)
    system(['open -a Finder "' fullfile(fpath,fname) '"']);
elseif isunix % Open the file in linux file manager - UNTESTED
    system(['xdg-open "' fullfile(fpath) '"']); % open folder in file manager
    % system(['xdg-open "' fullfile(fpath, fname) '"']); % open image in default viewer
end

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
classdef MatNatModality
    % MatNatModality Enumeration representing a modality for imaging data
    %
    %     Licence
    %     -------
    %     Part of MatNat. https://github.com/tomdoel/matnat
    %     Author: Tom Doel, 2015.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        Name
        XnatType
    end

    methods
        function obj = MatNatModality(name, xnatType)
            obj.Name = name;
            obj.XnatType = xnatType;
        end
    end

    enumeration
        MR ('MR', 'xnat:mrScanData')
    end

    properties
    end

    methods (Static)
        function modality = getModalityFromXnatString(xnatType)
            % Get the modality type from the Xnat type string

            allEnums = enumeration('MatNatModality');
            for enum = allEnums
                if strcmp(xnatType, enum.XnatType)
                    modality = enum;
                    return;
                end
            end
            modality = [];
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function Y = au_qmr_ilu_s(tr,X,i)
%
%  Solves shifted linear systems with the matrix A or its transpose A' by
%  the QMR iteration with ILU preconditioning:
%
%  for tr = 'N':
%
%    Y = inv(A+p(i)*I)*X;
%
%  for tr = 'T':
%
%    Y = inv(A.'+p(i)*I)*X.
%
%  The incomplete LU factors of the ILU preconditioner for A+p(i)*I are
%  provided as global data if mc has been set to 'C' when calling
%  'au_qmr_ilu_m_i'. In this case, this data must be generated by calling
%  'au_qmr_ilu_s_i' before calling this routine! If mc has been set to 'M',
%  the incomplete factorizations are computed within this routine to save
%  memory. However, the computational cost is increased in this case.
%
%  Calling sequence:
%
%    Y = au_qmr_ilu_s(tr,X,i)
%
%  Input:
%
%    tr        (= 'N' or 'T') determines whether shifted systems with
%              A or A' should be solved;
%    X         a matrix of proper size;
%    i         the index of the shift parameter.
%
%  Output:
%
%    Y         the resulting solution matrix.
%
%
%   LYAPACK 1.0 (Thilo Penzl, August 1999)

if nargin~=3
  error('Wrong number of input arguments.');
end

global LP_A LP_P LP_MC LP_MAX_IT_QMR LP_TOL_QMR LP_TOL_ILU LP_INFO_QMR

if ~length(LP_A) | ~length(LP_P) | ~length(LP_MC) | ~length(LP_MAX_IT_QMR) ...
  | ~length(LP_TOL_QMR) | ~length(LP_TOL_ILU) | ~length(LP_INFO_QMR)
  error('This routine needs global data which must be generated by calling ''au_qmr_ilu_m_i'' first.');
end

[n,m] = size(X);

if LP_MC=='C'

  eval(lp_e( 'global LP_L',i,' LP_U',i ));

  is_init1 = eval(lp_e( 'length(LP_L',i,')' ));
  is_init2 = eval(lp_e( 'length(LP_U',i,')' ));
  if ~is_init1 | ~is_init2
    error('This routine needs global data which must be generated by calling ''au_qmr_ilu_s_i'' first.');
  end

else

  eval(lp_e( '[LP_L',i,',LP_U',i,'] = luinc(LP_A+LP_P(i)*speye(n),LP_TOL_ILU);' ));

%
%  This shows the amount of fill-in created by the ILU preconditioner:
%
  if LP_INFO_QMR >= 5
    nnA = nnz(LP_A);
    eval(lp_e( 'nnLU = nnz(LP_L',i,')+nnz(LP_U',i,');'  ));
    disp(sprintf('ILU-Ratio (nnz(L)+nnz(U))/nnz(A) = %6g',nnLU/nnA));
  end

end

Y = zeros(n,m);

if tr=='N'
  for j = 1:m
    eval(lp_e(  '[Y(:,j),fl,rr] = qmr(LP_A+LP_P(i)*speye(n),X(:,j),LP_TOL_QMR,LP_MAX_IT_QMR,LP_L',i,',LP_U',i,');'  ));
    if LP_INFO_QMR & fl
      disp('WARNING in ''au_qmr_ilu_s'': QMR terminated not properly!');
      disp(sprintf('  ( normalized precond. QMR residual = %6g )',rr));
      rr0 = norm(LP_A*Y(:,j)+LP_P(i)*Y(:,j)-X(:,j))/norm(X(:,j));
      disp(sprintf('  ( normalized UNprecond. QMR residual = %6g )',rr0));
    end
  end
elseif tr=='T'
  for j = 1:m
    eval(lp_e(  '[Y(:,j),fl,rr] = qmr((LP_A.''+LP_P(i)*speye(n)),X(:,j),LP_TOL_QMR,LP_MAX_IT_QMR,LP_U',i,'.'',LP_L',i,'.'');'  ));
    if LP_INFO_QMR & fl
      disp('WARNING in ''au_qmr_ilu_s'': QMR terminated not properly!');
      disp(sprintf('  ( normalized QMR residual = %6g )',rr));
      rr0 = norm(LP_A.'*Y(:,j)+LP_P(i)*Y(:,j)-X(:,j))/norm(X(:,j));
      disp(sprintf('  ( normalized UNprecond. QMR residual = %6g )',rr0));
    end
  end
else
  error('tp must be either ''N'' or ''T''.');
end

author: Sable

repo: mcbench-benchmarks


text:
function r = sin(a)
%SIN          Hessian (elementwise) sine
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  K = prod(size(a.x));
  if K==1                   % scalar hessian

    r.x = sin(a.x);
    cosax = cos(a.x);
    r.dx = cosax * a.dx;
    r.hx = cosax * a.hx - reshape( ((0.5*r.x)*a.dx) * a.dx.' , size(a.hx) );

  else                      % matrix hessian

    N = getappdata(0,'INTLAB_HESSIAN_NUMVAR');
    N2 = N^2;

    r.x = sin(a.x);
    if issparse(a.hx)               % input sparse

      ax = cos(full(a.x(:)));
      sizeax = length(ax);
      [ia,ja,sa] = find(a.dx);
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if isempty(ia)
        r.dx = sparse([],[],[],N,sizeax);
        r.hx = sparse([],[],[],N2,sizeax);
      else
        adx1 = (-0.5)*r.x(:);
        if isa(a.x,'intval')          % sparse intval
          rdx = times(ax(ja),sa(:),0);
          adx1 = times(adx1(ja),sa(:),0);
          if rdx.complex
            r.dx = intval( sparse(ia,ja,rdx.mid,N,sizeax) , sparse(ia,ja,rdx.rad,N,sizeax) , 'midrad' );
          else
            r.dx = intval( sparse(ia,ja,rdx.inf,N,sizeax) , sparse(ia,ja,rdx.sup,N,sizeax) , 'infsup' );
          end
          if adx1.complex
            adx1 = intval( sparse(ia,ja,adx1.mid,N,sizeax) , sparse(ia,ja,adx1.rad,N,sizeax) , 'midrad' );
          else
            adx1 = intval( sparse(ia,ja,adx1.inf,N,sizeax) , sparse(ia,ja,adx1.sup,N,sizeax) , 'infsup' );
          end
        else                          % sparse point
          r.dx = sparse(ia,ja,ax(ja).*sa(:),N,sizeax);
          adx1 = sparse(ia,ja,adx1(ja).*sa(:),N,sizeax);
        end
        r.hx = adx2rhx(N,sizeax,adx1,a.dx);
      end
      [ia,ja,sa] = find(a.hx);        % sparse point or intval
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if ~isempty(ia)
        if isa(a.x,'intval')
          rhx = times(ax(ja),sa(:),0);
          if rhx.complex
            r.hx = r.hx + intval( sparse(ia,ja,rhx.mid,N2,sizeax) , sparse(ia,ja,rhx.rad,N2,sizeax) , 'midrad' );
          else
            r.hx = r.hx + intval( sparse(ia,ja,rhx.inf,N2,sizeax) , sparse(ia,ja,rhx.sup,N2,sizeax) , 'infsup' );
          end
        else
          r.hx = r.hx + sparse(ia,ja,ax(ja).*sa(:),N2,sizeax);
        end
      end

    else                            % input full

      r.x = sin(a.x);
      rx = r.x(:).';
      cosax = cos(a.x(:).');
      cosax = cosax(ones(N*N,1),:);
      r.dx = a.dx .* cosax(1:N,:);
      adx = repmat(0.5*rx,N,1) .* a.dx;
      r.hx = a.hx .* cosax - adx(repmat(1:N,N,1),:) .* a.dx(repmat(1:N,1,N),:);

    end

  end

  r = class(r,'hessian');

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
%% boxconstrain
% Below is a demonstration of the features of the |boxconstrain| function

%% Syntax
% |[xb]=boxconstrain(x,lb,ub,m,s);|

%% Description
% The |boxconstrain| function can be used to constrain parameters from
% [-inf,inf] to the range [lb ub]. The sigmoidal tanh function is used to
% do this mapping. The middle of the sigmoid, denoted by the parameter m
% need not be (lb+ub)/2. If x=m xb=m.
% The optional parameter s sets the slope of the sigmoid at
% x=m, the default slope is 1.

%% Examples

clear; close all; clc;

%% Example: Constraining parameters

m=5; %The middle
lb=0; %The lower bound
ub=10; %The upper bound
x=linspace(m-20,m+20,1000); % x paremeters for plotting
[y]=boxconstrain(x,lb,ub,m);

%%

cFigure; hold on;
xlabel('x'); ylabel('Constrained x');
h1=plot([min(x) max(x)],[ub ub],'b--','LineWidth',3);
h2=plot([min(x) max(x)],[lb lb],'r--','LineWidth',3);
h3=plot(m,m,'k.','MarkerSize',50);
h4=plot(x,y,'k-','MarkerSize',25,'LineWidth',3);

legend([h1;h2;h3;h4],{'Upper limit','Lower limit','Middle','Constrained'},'Location','Northwest');

axis tight;
set(gca,'FontSize',25);
grid on; box on;
drawnow;

%% Example: Assymetric mappings

m=1; %The middle
lb=0; %The lower bound
ub=10; %The upper bound

x=linspace(m-20,m+20,1000); % x paremeters for plotting

[y]=boxconstrain(x,lb,ub,m);

%%

cFigure; hold on;
xlabel('x'); ylabel('Constrained x');
h1=plot([min(x) max(x)],[ub ub],'b--','LineWidth',3);
h2=plot([min(x) max(x)],[lb lb],'r--','LineWidth',3);
h3=plot(m,m,'k.','MarkerSize',50);
h4=plot(x,y,'k-','MarkerSize',25,'LineWidth',3);

legend([h1;h2;h3;h4],{'Upper limit','Lower limit','Middle','Constrained'},'Location','Northwest');

axis tight;
set(gca,'FontSize',25);
grid on; box on;
drawnow;

%% Example: Custom slop

m=5; %The middle
lb=0; %The lower bound
ub=10; %The upper bound
s=0.5; %Slope at the middle

x=linspace(m-20,m+20,1000); % x paremeters for plotting

[y]=boxconstrain(x,lb,ub,m,s);

%%

cFigure; hold on;
xlabel('x'); ylabel('Constrained x');
h1=plot([min(x) max(x)],[ub ub],'b--','LineWidth',3);
h2=plot([min(x) max(x)],[lb lb],'r--','LineWidth',3);
h3=plot(m,m,'k.','MarkerSize',50);
h4=plot(x,y,'k-','MarkerSize',25,'LineWidth',3);

legend([h1;h2;h3;h4],{'Upper limit','Lower limit','Middle','Constrained'},'Location','Northwest');

axis tight;
set(gca,'FontSize',25);
grid on; box on;
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%  INTERNAL FUCNTION
%
%

author: jmaih


text:
function planC = CERRImport(optName, archiveDir, saveFile)
%function planC = CERRImport(optName,archiveDir,saveDir)
%
%Description: Wrapper file to import AAPM/RTOG directories.  The main work is done
%by importRTOGDir.  If called without arguments, the user selects an archive file and
%an options file via GUI dialogs, and the name of the file to save the
%CERR plan to is input via the file save GUI.  Can be called with 1, 2, or 3
%arguments.
%
%Inputs:
%optName -- is the name of an options file (an example file is initCERRImport).
%archiveDir -- is the path of the archive directory, with trailing slash.
%saveFile -- the name, including path, of a file to which to save the CERR plan.
%
%Output: None, except the saved file.
%
%Globals:  None.
%
%Author: J. O. Deasy, deasy@radonc.wustl.edu
%
%Reference:
%W. Harms, Specifications for Tape/Network Format for Exchange of Treatment Planning Information,
%version 3.22., RTOG 3D QA Center, (http://rtog3dqa.wustl.edu), 1997.
%
%Last Modified:  14 Feb 02, by JOD.
%
%Version:  1.x. (Incremented by 0.001 for every modification of this file or a called file.)
%
%Copyright:
%This software is copyright J. O. Deasy and Washington University in St Louis.
%A free license is granted to use or modify but only for non-commercial non-clinical use.
%Any user-modified software must retain the original copyright and notice of changes if redistributed.
%Any user-contributed software retains the copyright terms of the contributor.

% copyright (c) 2001-2006, Washington University in St. Louis.
% Permission is granted to use or modify only for non-commercial,
% non-treatment-decision applications, and further only if this header is
% not removed from any file. No warranty is expressed or implied for any
% use whatever: use at your own risk.  Users can request use of CERR for
% institutional review board-approved protocols.  Commercial users can
% request a license.  Contact Joe Deasy for more information
% (radonc.wustl.edu@jdeasy, reversed).

fname1 = '';
%Interactive or not:
if nargin == 0
    [fname1, pathname] = uigetfile('*.*','Select 0000 file to begin importing.');
    archiveDir = pathname;

    [fname, pathname] = uigetfile('*.m','Select CERR options .m file');
    optName = [pathname fname];

    %   % Would the user like to store the plan as a compressed file?
    %   % The following loop investigates this.
    %   % Compression (gz) is made through GNUZIP (freeware)
    %   % Default compression options are 'verbose' and 'best'
    %   gzip 1.2.4 Win32 (02 Dec 97)
    %   usage: gzip [-acdfhlLnNrtvV19] [-S suffix] [file ...]
    %       -a --ascii       ascii text; convert end-of-lines using local conventions
    %       -c --stdout      write on standard output, keep original files unchanged
    %       -d --decompress  decompress
    %       -f --force       force overwrite of output file and compress links
    %       -h --help        give this help
    %       -l --list        list compressed file contents
    %       -L --license     display software license
    %       -n --no-name     do not save or restore the original name and time stamp
    %       -N --name        save or restore the original name and time stamp
    %       -q --quiet       suppress all warnings
    %       -r --recursive   operate recursively on directories
    %       -S .suf  --suffix .suf     use suffix .suf on compressed files
    %       -t --test        test compressed file integrity
    %       -v --verbose     verbose mode
    %       -V --version     display version number
    %       -1 --fast        compress faster
    %       -9 --best        compress better
    %       file...          files to (de)compress. If none given, use standard input.
    % Logic for compressing files:
    %   Matlab binaries will be created first in all cases; subsequently,
    %       they will be replaced by the corresponding compressed binary IF
    %       the user has selected this option in 'CERROptions.m'

    %fid = fopen(optName,'r');
    %while(1)
    %    tline = fgetl(fid);
    %    if findstr((tline), 'optS.saveZip'), break, end
    %end
    %fclose(fid);
    %index = findstr(tline, '=');
    %toZipOrNot = lower(midstring(tline, (index + 1), length(tline) - 1));
    %if findstr(toZipOrNot, 'yes')
    %    toZipOrNot = 'yes';
    %elseif findstr(toZipOrNot, 'no')
    %    toZipOrNot = 'no';
    %else
    warndlg('Output file compression not set: defaulting to full.', ...
        'CERR Configuration');
end
optS = opts4Exe(optName);

if strcmp(optS.saveZip,'yes')

    %end        % APA commented this end (may be placed incorrectly, gave
    %an error during crosslisting)

    saveFlag = 'ui';

elseif nargin == 1 %Assume options given
    [fname1, pathname] = uigetfile('*.*','Select aapm0000 file to begin importing.');
    archiveDir = pathname;
    fid = fopen(optName,'r');
    while(1)
        tline = fgetl(fid);
        if strfind((tline), 'optS.saveZip'), break, end
    end
    fclose(fid);
    index = strfind(tline, '=');
    toZipOrNot = lower(midstring(tline, (index + 1), length(tline) - 1));
    if strfind(toZipOrNot, 'yes')
        toZipOrNot = 'yes';
    elseif strfind(toZipOrNot, 'no')
        toZipOrNot = 'no';
    else
        warndlg('Output file compression not set: defaulting to full.', ...
            'CERR Configuration');
    end
    optS = opts4Exe(optName);

    saveFlag = 'ui';

elseif nargin == 2
    fid = fopen(optName,'r');
    while(1)
        tline = fgetl(fid);
        if strfind((tline), 'optS.saveZip'), break, end
    end
    fclose(fid);
    index = strfind(tline, '=');
    toZipOrNot = lower(midstring(tline, (index + 1), length(tline) - 1));
    if strfind(toZipOrNot, 'yes')
        toZipOrNot = 'yes';
    elseif strfind(toZipOrNot, 'no')
        toZipOrNot = 'no';
    else
        warndlg('Output file compression not set: defaulting to full.', ...
            'CERR Configuration');
    end
    optS = opts4Exe(optName);

    saveFlag = 'ui';

elseif nargin == 3
    fid = fopen(optName,'r');
    while(1)
        tline = fgetl(fid);
        if strfind((tline), 'optS.saveZip'), break, end
    end
    fclose(fid);
    index = strfind(tline, '=');
    toZipOrNot = lower(midstring(tline, (index + 1), length(tline) - 1));
    if strfind(toZipOrNot, 'yes')
        toZipOrNot = 'yes';
    elseif strfind(toZipOrNot, 'no')
        toZipOrNot = 'no';
    else
        warndlg('Output file compression not set: defaulting to full.', ...
            'CERR Configuration');
    end
    optS = opts4Exe(optName);
    if ~strcmp(saveFile,'no')
        saveFlag = 'spec';
    else
        saveFlag = 'no';
    end
    fname1 = 'aapm0000';

end

%Call the import tool:
planC = importRTOGDir(optS, archiveDir, fname1)

switch saveFlag

    case 'ui'

        [fname pname] = uiputfile( '*.mat', ...
            'Save the plan data as:');
        saveFile = [pname fname];
        if ~strcmp(midstring(saveFile, ...
                length(saveFile)-3, length(saveFile)), '.mat')
            % '*.mat' file type not typed by the user
            % this situation can create a file error when compressing
            %       will add the '*.mat' suffix prior to archive being written
            saveFile = strcat(saveFile, '.mat');
        end
        save(saveFile,'planC')
        if strcmp(toZipOrNot, 'yes') % User wants to compress output file
            outstr = gnuCERRCompression(saveFile,'compress');
            h = dos(outstr);
            if (h) % Function returned as zero when compression complete
                errordlg('File not found by compression routine.', 'CERR compression');
            end
        end

    case 'spec'

        if isempty(saveFile) %save within the subdirectory
            saveFile = [archiveDir 'planC.mat'];
            if strcmp(toZipOrNot, 'yes') % User wants to compress output file
                outstr = gnuCERRCompression(saveFile,'compress');
                h = dos(outstr);
                if (h) % Function returned as zero when compression complete
                    errordlg('File not found by compression routine.', 'CERR compression');
                end
            end
        end

        try
            save(saveFile,'planC')
            if strcmp(toZipOrNot, 'yes') % User wants to compress output file
                outstr = gnuCERRCompression(saveFile,'compress');
                h = dos(outstr);
                if (h) % Function returned as zero when compression complete
                    errordlg('File not found by compression routine.', 'CERR compression');
                end
            end

        catch
            warndlg('That file name was not valid.', 'CERR file save')
            [fname pname] = uiputfile('*.mat','Save the plan data as:');
            saveFile = [pname fname];
            save(saveFile,'planC')
            if strcmp(toZipOrNot, 'yes') % User wants to compress output file
                outstr = gnuCERRCompression(saveFile,'compress');
                h = dos(outstr);
                if (h) % Function returned as zero when compression complete
                    errordlg('File not found by compression routine.', 'CERR compression');
                end
            end

        end
    case 'no'

end

% Unzip routines
% outstr = ['gzip.exe -vd ', fmat]
% dos(outstr)

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function angle_axis = RotationMatrix2AngleAxis(R)

% The conversion of a rotation matrix to the angle-axis form is
% numerically problematic when then rotation angle is close to zero
% or to Pi. The following implementation detects when these two cases
% occurs and deals with them by taking code paths that are guaranteed
% to not perform division by a small number.

% x = k * 2 * sin(theta), where k is the axis of rotation.
angle_axis(1) = R(1+2, 1+1) - R(1+1, 1+2);
angle_axis(2) = R(1+0, 1+2) - R(1+2, 1+0);
angle_axis(3) = R(1+1, 1+0) - R(1+0, 1+1);

% Since the right hand side may give numbers just above 1.0 or
% below -1.0 leading to atan misbehaving, we threshold.
costheta = min(max((R(1+0, 1+0) + R(1+1, 1+1) + R(1+2, 1+2) - 1.0) / 2.0,  -1.0), 1.0);

% sqrt is guaranteed to give non-negative results, so we only
% threshold above.
sintheta = min(sqrt(angle_axis(1) * angle_axis(1) + angle_axis(2) * angle_axis(2) + angle_axis(3) * angle_axis(3)) / 2.0, 1.0);

% Use the arctan2 to get the right sign on theta
theta = atan2(sintheta, costheta);

% Case 1: sin(theta) is large enough, so dividing by it is not a
% problem. We do not use abs here, because while jets.h imports
% std::abs into the namespace, here in this file, abs resolves to
% the int version of the function, which returns zero always.
%
% We use a threshold much larger then the machine epsilon, because
% if sin(theta) is small, not only do we risk overflow but even if
% that does not occur, just dividing by a small number will result
% in numerical garbage. So we play it safe.
kThreshold = 1e-12;
if ((sintheta > kThreshold) || (sintheta < -kThreshold))
    r = theta / (2.0 * sintheta);
    angle_axis = angle_axis * r;
    return;
end

% Case 2: theta ~ 0, means sin(theta) ~ theta to a good
% approximation.
if (costheta > 0.0)
    angle_axis = angle_axis * 0.5;
    return;
end

% Case 3: theta ~ pi, this is the hard case. Since theta is large,
% and sin(theta) is small. Dividing by theta by sin(theta) will
% either give an overflow or worse still numerically meaningless
% results. Thus we use an alternate more complicated formula
% here.

% Since cos(theta) is negative, division by (1-cos(theta)) cannot
% overflow.
inv_one_minus_costheta = 1.0 / (1.0 - costheta);

% We now compute the absolute value of coordinates of the axis
% vector using the diagonal entries of R. To resolve the sign of
% these entries, we compare the sign of angle_axis[i]*sin(theta)
% with the sign of sin(theta). If they are the same, then
% angle_axis[i] should be positive, otherwise negative.

for i=1:3
    angle_axis(i) = theta * sqrt((R(i, i) - costheta) * inv_one_minus_costheta);
    if (((sintheta < 0.0) && (angle_axis(i) > 0.0)) || ((sintheta > 0.0) && (angle_axis(i) < 0.0)))
        angle_axis(i) = -angle_axis(i);
    end
end

author: jianxiongxiao

repo: ProfXkit


text:
function s = rdivide(f, g)
%./   Divide DELTAFUNS with DELTAFUNS
%   RDIVIDE(F, G) computes the pointwise division F./G. The operation is only
%   defined if G does not have any delta functions and G has no roots.
%
% See also LDIVIDE, TIMES.
%
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Case of empty arguments:
if ( isempty(f) || isempty(g) )
    % Return an empty DELTAFUN:
    s = deltafun();
    return
end

% Check if inputs are other than DELTAFUNs, CLASSICFUNs or doubles.
if ( (~isa(f, 'deltafun') && ~isa(f, 'classicfun') && ~isa(f, 'double')) || ...
     (~isa(g, 'deltafun') && ~isa(g, 'classicfun') && ~isa(g, 'double')) )
    error('CHEBFUN:DELTAFUN:rdivide:rdivide', ...
        'Input can only be a DELTAFUN, a CLASSICFUN or a double')
end

%% Reciprocal: ( CLASSICFUN or DOUBLE ). / DELTAFUN
if ( isa(f, 'double') || isa( f, 'classicfun') && isa(g, 'deltafun') )
    if ( anyDelta(g) )
        error('CHEBFUN:DELTAFUN:rdivide:rdivide', ...
            'Division by delta functions is not defined.');
    end
    % A smooth function is returned in this case:
    s = f ./ g.funPart;
    return
end

%% DELTAFUN./DOUBLE
if ( isa(f, 'deltafun') && isa(g,'double') )
    % Copy the DELTAFUN in the output:
    s = f;
    % Divide the funPart with the double g:
    s.funPart = 1/g * f.funPart;
    % Update the delta functions:
    s.deltaMag = 1/g * s.deltaMag;
    if ( ~anyDelta(s) )
        s = s.funPart;
    end
    return
end

%% DELTAFUN./CLASSICFUN
if ( isa(f, 'deltafun') && isa(g, 'classicfun') )
    % Take reciprocal and make a DELTAFUN:
    g = deltafun( 1./g, [], [] );
    % Now, multiply:
    s = f .* g;
    return
end

%% DELTAFUN./DELTAFUN
if ( isa(f, 'deltafun') && isa(g, 'deltafun') )
    if ( anyDelta(g) )
        error('CHEBFUN:DELTAFUN:rdivide:rdivide', ...
            'Division by delta functions is not defined.');
    end
    % Take reciprocal of the funPart only and make a DELTAFUN:
    g = deltafun( 1 ./ g.funPart, [], []);
    % Now, multiply:
    s = f .* g;
    return
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function sp=findConmatFromEdges(nodeList,edgeList)
% function connectionMatrix=findConmatFromEdges(edgeList)
% Generates a connection matrix from a list of edges
% The returned connection matrix is max(edgeList(:)) square
%
% nodeList is an 8xn array
% nodes:  8xN array of
%    Nx(x,y,z,num_edges,edge_offset,layer,dist,pqindex).
% edges:  1xM array of node indices.  The edge_offset of
%    each node points into the starting location of its set
%    of edges.

nNodes=size(nodeList,2);

mxNode=max(edgeList(:));
sp=spalloc(nNodes,nNodes,fix(nNodes/40));

for t=1:nNodes

	nEdges=nodeList(4,t);
	offset=nodeList(5,t);

	theseNodes=edgeList(offset:(offset+nEdges-1));
	sp(t,theseNodes)=1;
end

sp=((sp+sp')~=0); % Matrix is symmetrical

author: vistalab

repo: vistasoft


text:
function triplot
    % triplot plots a triangular diagram
    % points on it are plotted by tripts
    % written by Gerry Middleton, 1995
    xscale = 2/sqrt(3);
    x = [0 xscale/2 xscale 0];		 % define coordinates of apices
    y = [-1 0 -1 -1];
    plot(x,y,'k');
    axis('equal');
    axis off;						 % turn off rectangular axes
    hold on;						 % hold, to permit adding to plot later
    for i = 1:9						 % plot dotted lines parallel to sides,
        x1 = i*xscale/20;			 % at 10 percent intervals
        x2 = (20-i)*xscale/20;
        x3 = i*xscale/10;
        y1 = -1 + i/10;
        y2 = -i/10;
        xx = [x1 x2];
        yy = [y1 y1];
        xx2 = [x1 x3];
        yy2 = [y1 -1];
        xx3 = [x3 (xscale/2 + x1)];
        yy3 = [-1 y2];
        plot(xx,yy,'k:',xx2,yy2,'k:',xx3,yy3,'k:');
    end
    set(gcf,'color','w')

author: CelsoReyes


text:
function varargout = erfcinv(varargin)

switch class(varargin{1})

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

    case 'char'

        operator = CreateBasicOperator('decreasing','callback');
        operator.inverse = @(x)(erfc(x));
        operator.derivative = @(x)-1./(exp(-((erfcinv(x))).^2)*2/sqrt(pi));
        operator.inflection = [1 -1];
        operator.domain = [0 2];

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

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


function [dy edy dx edx y ey x ex tmpw tmph] = pad(total_boxes,w,h)
	%compute the padding coordinates (pad the bounding boxes to square)
	tmpw=total_boxes(:,3)-total_boxes(:,1)+1;
	tmph=total_boxes(:,4)-total_boxes(:,2)+1;
	numbox=size(total_boxes,1);

    dx=ones(numbox,1);dy=ones(numbox,1);
	edx=tmpw;edy=tmph;

	x=total_boxes(:,1);y=total_boxes(:,2);
	ex=total_boxes(:,3);ey=total_boxes(:,4);

	tmp=find(ex>w);
	edx(tmp)=-ex(tmp)+w+tmpw(tmp);ex(tmp)=w;

	tmp=find(ey>h);
	edy(tmp)=-ey(tmp)+h+tmph(tmp);ey(tmp)=h;

	tmp=find(x<1);
	dx(tmp)=2-x(tmp);x(tmp)=1;

	tmp=find(y<1);
	dy(tmp)=2-y(tmp);y(tmp)=1;
end


text:
function C = glm_determine_noise(Y,E,tr)
% Estimate the noise covariance for each voxel using its time series
%
%    C = glm_determine_noise(Y,E,tr):
%
% This routine is called when we try to whiten the noise in the
% event-related GLM tools.
%
% Given a matrix of time series, Y, and an estimated error for each voxel,
% E, and a frame period, tr, estimate the temporal noise covariance matrix
% for each voxel, C.
%
% Size of each argument:
% Input: Y: # Time Points x # Voxels x # Runs
%        E:
%        tr: single integer
% Output: C: # Predictors x # Predictors x # Runs
%         The # of predictors is the total # of
%         different predictors in a GLM, including
%         DC baseline estimate predictors and, if
%         a selective-averaging GLM is being applied,
%         all nh hemodynamic predictors for each condition.
%         (The # predictors is the # of columns in the
%         design matrix produced by er_createDesMtx.)
%
%
% original code by gb, 11/04:
% updated by ras, 05/05

if nargin<3, help(mfilename); error('Insufficient args.'); end

% params
delayRange = 20; % # of seconds over which to estimate noise
kmax = floor(delayRange / tr);
nFrames = size(E,1);
nSlices = size(E,2);

% use a heuristic to select voxels within the skull:
withinSkullVoxels = glm_choose_voxels(Y);

% estimate the autocorrelation within these
% voxels over the selected delay range (e.g., 20 secs):
B = E(:,withinSkullVoxels); % voxels in brain
for k = 1:kmax
    corr(k,:) = sum(B(1:nFrames-k,:).*B(k+1:nFrames,:))/sum(B.^2);
end
corr = mean(corr,2);

% Determine the parameters alpha and rho using in estimating the
% autocorrelation of residual errors (Re). These parameters will be used in
% the estimation function:
% Re(k) = 1,  k=0;
%         (1-alpha) * rho^k, 0 < |k| <= kmax
%         0, |k| > kmax
% This is equation (5) in the Burock and Dale 2000 paper, and eq. (13) in
% the Greve FS-FAST theory paper.
kmax = length(corr);

[x eval] = fminsearch(@glm_fit_model, [0, 1/2], optimset('TolX',1e-8),kmax,corr);
alpha = x(1);
rho   = x(2);

% finally, compute the noise covariance matrix:
tmp = [1 ((1 - alpha)*(rho.^(1:kmax))) zeros(1,nFrames - kmax - 1)];
C = toeplitz(tmp);

return
% /-----------------------------------------------------/ %

% sum(E(1:nFrames-k,:).*e(k+1:nFrames,:))/sum(E.^2);


text:
function [connectivity] = channelconnectivity(cfg, data)

% CHANNELCONNECTIVIY creates a NxN matrix that describes whether channels
% are connected as neighbours
%
% See also FT_PREPARE_NEIGHBOURS, TRIANGLE2CONNECTIVITY

if (isfield(cfg, 'avgoverchan') && strcmp(cfg.avgoverchan, 'yes')) || isempty(cfg.neighbours)
  if nargin < 2
    nchan = numel(cfg.channel);
  else
    nchan = numel(data.label);
  end
  connectivity = false(nchan,nchan);

else
  % check whether to use channels in the cfg or data
  if nargin < 2 && isfield(cfg, 'channel')
    chans = cfg.channel;
  elseif nargin == 2
    chans = data.label;
  else
    ft_error('either the cfg needs to have both cfg.channel and cfg.neighbours, or a second (data) input argument needs to be specified');
  end

  nchan        = length(chans);
  connectivity = false(nchan,nchan);
  neighbours   = struct(cfg.neighbours);

  try
    % the original loop was the following, which is very slow:
    %   for chan=1:length(cfg.neighbours)
    %     [seld] = match_str(chans, cfg.neighbours(chan).label);
    %     [seln] = match_str(chans, cfg.neighbours(chan).neighblabel);
    %     if isempty(seld)
    %       % this channel was not present in the data
    %       continue;
    %     else
    %       % add the neighbours of this channel to the matrix
    %       connectivity(seld, seln) = true;
    %     end
    %   end
    % the above is still used in case the fast version fails (catch below)
    % the following code should be equivalent:

    if numel(neighbours) < nchan/4
      ft_error('this only works when at least 1/4th of all channels has neighbours defined, or with neighbours = []');
      % FIXME the above error is true (or maybe with <1/5th) but I don't
      % understand why, ES 25-nov-2013
      % the above error will cause the code to fall through to the catch
      % below, so no big worries
    end

    % avoid using match_str inside the loop for performance reasons
    [sel1,sel2] = match_str(chans, {neighbours.label}');

    % make sure we only have the neighbours present in the data
    neighbours = neighbours(sel2);

    % make big list of all neighbour labels...
    allneighb = {neighbours.neighblabel};

    % determine whether neighbours are stored as row or column vector
    siz1 = cellfun(@size, allneighb, repmat({1},size(allneighb)));
    siz2 = cellfun(@size, allneighb, repmat({2},size(allneighb)));
    if sum(siz2) > sum(siz1)
      catdim = 2;
    else
      catdim = 1;
    end

    % ...store dimensionality...
    numAllNeighb = cellfun(@numel, allneighb);
    allneighb = cat(catdim, allneighb{:});
    % ...so that we can do with one call of match_str (with fullout=1)...
    selneighb = match_str(chans, allneighb, 1);
    % ...and then we can reshape the indices again to get the neighbours
    selneighb = mat2cell(selneighb, numAllNeighb);
    for k = 1:numel(sel1)
      % remove zeroes from selneighb (these correspond to channels absent)
      seln = selneighb{k};
      seln = seln(seln > 0);
      % add the neighbours of this channel to the matrix
      connectivity(sel1(k), seln) = true;
    end

  catch
    % the fast version failed, use slow instead

    for chan=1:length(neighbours)
      [seld] = match_str(chans, neighbours(chan).label);
      [seln] = match_str(chans, neighbours(chan).neighblabel);
      if isempty(seld)
        % this channel was not present in the data
        continue;
      else
        % add the neighbours of this channel to the matrix
        connectivity(seld, seln) = true;
      end
    end

  end
end


text:
function s = sem(x,varargin)

%sem - Compute standard error of the mean (SEM).
%
%  USAGE
%
%    s = sem(x)
%
%    x              vector or matrix over which the SEM should be computed
%    dim            dimension to run SEM across (default = 1)

% Copyright (C) 2008-2013 by Michaël Zugaro
%
% 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.

if nargin < 1,
  error('Incorrect number of parameters (type ''help sem'' for details).');
end
if nargin == 1
    dim = 1;
else
    dim = varargin{1};
end
if ~isnumeric(x) & ~isdvector(x),
  error('Incorrect input - use vector or matrix (type ''help sem'' for details).');
end

% if any(size(x)==1), x = x(:); end

% s = std(x,[],dim)./sqrt(size(x,1));
s = nanstd(x,[],dim) ./ sqrt(sum(~isnan(x),dim));

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [y] = spm_gx_state_fmri(x,u,P,M)
% Simulated BOLD response and copied state vector
% FORMAT [y] = spm_gx_state_fmri(x,u,P,M)
% y          - BOLD response and copied state vector
%
% x          - state vector     (see spm_fx_fmri)
% P          - Parameter vector (see spm_fx_fmri)
% M          - model specification structure (see spm_nlsi)
%
% The `copied state vector' passes the first hidden variable in each region
% to the output variable y, so that 'neural activities' can be plotted
% by spm_dcm_generate.m
%
% See spm_fx_fmri.m and spm_dcm_generate.m
%__________________________________________________________________________
% Copyright (C) 2011 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_gx_state_fmri.m 6262 2014-11-17 13:47:56Z karl $

y = spm_gx_fmri(x,u,P,M);

% Copy first hidden state (neural activity) from each region
nregions=size(x,1);
for i=1:nregions,
    y=[y;x(i,1)];
end
y=full(y);


text:
function pstruct = tapas_hgf_ar1_binary_mab_namep(pvec)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2017 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/>.

pstruct = struct;

l = (length(pvec)+1)/6;

if l ~= floor(l)
    error('tapas:hgf:UndetNumLevels', 'Cannot determine number of levels');
end

pstruct.mu_0      = pvec(1:l);
pstruct.sa_0      = pvec(l+1:2*l);
pstruct.phi       = pvec(2*l+1:3*l);
pstruct.m         = pvec(3*l+1:4*l);
pstruct.ka        = pvec(4*l+1:5*l-1);
pstruct.om        = pvec(5*l:6*l-1);

return;

author: translationalneuromodeling

repo: tapas